Implement Zou+2024 Mbh-Mstar scaling relation (48b657df) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

Makefile

+24 −1

Original line number	Diff line number	Diff line
		@@ -219,9 +219,17 @@ fig/plambda_qu_sf_ctn_ctk.pdf: python/plot_plambda_qu_sf_ctn_ctk.py
		fig/plambda_qu_sf_ctn_ctk_20251006.pdf: src/scripts/plots/plot_plambda_qu_sf_ctn_ctk.py
		python3 $< $@

		fig/mbh_mstar.pdf: src/python/plot_mbh_mstar.py
		fig/plambda_qu_sf_ctn_ctk_20251025.pdf: src/scripts/plots/plot_plambda_qu_sf_ctn_ctk.py
		python3 $< $@


		#fig/mbh_mstar.pdf: src/python/plot_mbh_mstar.py
		fig/mbh_mstar.pdf: src/scripts/plots/plot_mbh_mstar.py
		python3 $< $@

		fig/mbh_mstar_with_zou2024.pdf: src/scripts/plots/plot_mbh_mstar.py
		python3 $< $@ --plot_zou2024

		fig/xray_luminosity_function.pdf: python/plot_xray_luminosity_function.py
		python3 $< $@

		@@ -279,5 +287,20 @@ figures_20250930:
		# gband number counts
		python3 src/scripts/plots/plot_number_counts_gband.py

		figures_20251024:
		# plambda
		python3 src/scripts/plots/plot_plambda_qu_sf_ctn_ctk.py fig/plambda_qu_sf_ctn_ctk_20251024.pdf
		# xlf
		python3 src/scripts/plots/plot_xray_luminosity_function.py data/catalog/dr1_new_new/ fig/xray_luminosity_function_20251024.pdf 0
		python3 src/scripts/plots/plot_xray_luminosity_function.py data/catalog/dr1_new_new/ fig/xray_luminosity_function_focc_20251024.pdf 1
		# qlf
		python3 src/scripts/plots/plot_quasar_luminosity_function.py fig/quasar_luminosity_function_20251024.pdf 0
		python3 src/scripts/plots/plot_quasar_luminosity_function.py fig/quasar_luminosity_function_focc_20251024.pdf 1
		# bhmf
		python3 src/scripts/plots/plot_black_hole_mass_function.py
		# gband number counts
		python3 src/scripts/plots/plot_number_counts_gband.py fig/number_counts_lsst_gband_20251024.pdf


		coverage:
		pytest --cov=src --cov-report=html tests/

data/mbh_mstar_zou2024/table1.txt

0 → 100644

+21 −0

Original line number	Diff line number	Diff line
		# logMstar_min logMstar_max logMBH(z=0.0) logMBH(z=0.5) logMBH(z=1.0) logMBH(z=1.5) logMBH(z=2.0) logMBH(z=2.5) logMBH(z=3.0) logMBH(z=3.5) logMBH(z=4.0)
		9.50 9.75 5.99 6.00 6.02 6.00 5.96 5.96 5.99 5.95 6.00
		9.50 9.75 5.99 5.99 6.02 5.99 5.96 5.95 5.96 5.96 5.98
		9.75 10.00 6.17 6.28 6.34 6.29 6.17 6.13 6.08 6.12 6.24
		9.75 10.00 6.17 6.26 6.33 6.29 6.16 6.11 6.09 6.11 6.25
		10.00 10.25 6.44 6.60 6.71 6.76 6.62 6.46 6.41 6.44 6.54
		10.00 10.25 6.43 6.55 6.69 6.79 6.58 6.43 6.38 6.42 6.52
		10.25 10.50 6.84 6.96 7.14 7.21 7.06 6.90 6.84 6.89 6.80
		10.25 10.50 6.75 6.88 7.10 7.18 7.04 6.85 6.82 6.84 6.76
		10.50 10.75 7.22 7.34 7.46 7.54 7.48 7.35 7.35 7.26 7.14
		10.50 10.75 7.21 7.32 7.45 7.55 7.46 7.34 7.37 7.27 7.01
		10.75 11.00 7.52 7.63 7.75 7.80 7.79 7.75 7.78 7.62 nan
		10.75 11.00 7.56 7.65 7.79 7.82 7.78 7.71 7.74 7.57 7.27
		11.00 11.25 7.86 7.96 8.07 8.05 8.19 8.12 8.04 nan nan
		11.00 11.25 7.91 7.99 8.10 8.16 8.11 8.09 8.05 7.88 7.52
		11.25 11.50 8.10 8.23 8.29 8.41 nan nan nan nan nan
		11.25 11.50 8.18 8.29 8.43 8.54 8.54 8.43 8.36 8.15 7.81
		11.50 11.75 8.36 8.55 nan nan nan nan nan nan nan
		11.50 11.75 8.49 8.59 8.79 8.93 8.94 8.83 8.65 nan nan
		11.75 12.00 nan nan nan nan nan nan nan nan nan
		11.75 12.00 8.82 9.00 9.16 9.37 nan nan nan nan nan

src/lsst_inaf_agile/mbh.py

+44 −0

Original line number	Diff line number	Diff line
		@@ -63,6 +63,49 @@ def get_log_mbh_const(log_mstar, A=500):
		return log_mstar - np.log10(A)


		def get_log_mbh_zou2024(log_mstar, z=0.0, is_tng300=False):
		"""
		Return Zou+2024 MBH-Mstar(z).

		The MBH-Mstar is based on a hybrid approach using the observed BHAR from
		earlier Zou work combined with mergers from IllustrisTNG. The Mbh-Mstar(z)
		is then built by following seeded black hole masses starting from z=4.0.

		Refer to Fig. 4 and Table 1 of Zou+2024:
		https://ui.adsabs.harvard.edu/abs/2024ApJ...976....6Z
		"""

		if not np.all(z == np.clip(z, 0, 4)):
		raise ValueError("Redshift must be in [0, 4]")

		assert np.all(z[0] == z)
		z = z[0]

		# Zou defines zbins with deltaz = 0.5
		zs = np.arange(0.0, 4.0 + 0.1, 0.5)
		print(zs)

		zclosest = zs[np.argmin(np.abs(zs - z))]
		print(f"Found {zclosest=} corresponding to {z}")

		def get_mstar_mbh(table, z):
		x1 = table["logMstar_min"][is_tng300::2]
		x2 = table["logMstar_max"][is_tng300::2]
		x = (x1 + x2) / 2
		y = [table[i + is_tng300][f"logMBH(z={z:.1f})"] for i in range(0, len(table), 2)]
		return np.array([x, y])

		table = Table.read("data/mbh_mstar_zou2024/table1.txt", format="ascii")
		x, y = get_mstar_mbh(table, z)

		# remove masked values...
		select = np.isfinite(y)
		x = x[select]
		y = y[select]

		return np.interp(log_mstar, x, y, left=np.nan, right=np.nan)


		def get_delta_log_mbh_shankar2019(log_mstar, low=0.0, high=np.inf):
		"""Return equation 6 from Shankar+2019."""
		delta_mbh = 0.32 - 0.1 * (log_mstar - 12)
		@@ -171,6 +214,7 @@ REFERENCES = (
		"continuity",
		"continuity_new",
		"continuity_new2",
		"Zou+24",
		)

src/scripts/plots/plot_mbh_mstar.py

0 → 100644

+104 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python3

		from argparse import ArgumentParser

		import matplotlib as mpl
		import matplotlib.pyplot as plt
		import numpy as np
		from lsst_inaf_agile import mbh

		###############################################################################
		# Command-line interface
		parser = ArgumentParser()
		parser.add_argument("savefig")
		parser.add_argument("--plot_zou2024", action="store_true")
		args = parser.parse_args()
		###############################################################################

		###############################################################################
		# Settings for the plot
		mpl.style.use("etc/agile.mplstyle")
		log_mstar = np.linspace(8.5, 12.5, 401)
		dlog_mstar = 0.10
		dz = 1.0
		lw1 = 2.5
		lw2 = 2.0
		alpha = 0.80
		cmap = mpl.colormaps["gist_heat"]
		cmap2 = mpl.colormaps["cool"]
		norm = mpl.colors.Normalize(vmin=0.0, vmax=5.5)
		###############################################################################


		def main():
		# fig, axes = plt.subplots(1, 2, figsize=(1.3 * 6.4, 1.0 * 4.8), sharex=True, sharey=True, dpi=300)
		fig, axes = plt.subplots(1, 2, figsize=(8.3, 5.85 / 1.5), sharex=True, sharey=True, dpi=300)

		# Left panel: Plots the redshift evolution of the Con Eq
		for z in [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]:
		# log_mbh = mbh.get_log_mbh_continuity(log_mstar, z)
		# log_mbh = mbh.get_log_mbh_continuity_new(log_mstar, z)
		ls = "solid" if z < 4 else "dotted"
		log_mbh = mbh.get_log_mbh_continuity_new2(log_mstar, np.full_like(log_mstar, z))
		axes[0].plot(log_mstar, log_mbh, color=cmap(norm(z)), lw=lw1, label=f"$z = {z:d}$", ls=ls)

		# Right panel: Plots the local Con Eq against literature
		# log_mbh = mbh.get_log_mbh_continuity(log_mstar, 0.0)
		# log_mbh = mbh.get_log_mbh_continuity_new(log_mstar, 0.0)
		log_mbh = mbh.get_log_mbh_continuity_new2(log_mstar, np.zeros_like(log_mstar))
		axes[1].plot(
		log_mstar,
		log_mbh,
		lw=lw1,
		color="k",
		# label="Continuity equation, $z=0$"
		)
		axes[1].fill_between(log_mstar, log_mbh - 0.50, log_mbh + 0.50, color="k", alpha=0.20)
		axes[1].plot(
		log_mstar,
		mbh.get_log_mbh_haring_rix2004(log_mstar),
		lw=lw2,
		# label="Häring & Rix 2004",
		label="H&R04",
		alpha=alpha,
		ls="dashdot",
		)
		axes[1].plot(
		log_mstar,
		mbh.get_log_mbh_kormendy_ho2013(log_mstar),
		lw=lw2,
		label="K&H13",
		alpha=alpha,
		ls="dashed",
		)
		axes[1].plot(
		log_mstar,
		mbh.get_log_mbh_reines_volonteri2015(log_mstar),
		lw=lw2,
		label="R&V15",
		alpha=alpha,
		ls="dotted",
		)

		for ax in axes:
		ax.set_xlim(9.5, 11.9)
		ax.set_ylim(6.1, 8.9)
		ax.set_xlabel(r"$\log_{10} \left( M_{\rm star} \,/\, M_\odot \right)$")
		# ax.set_aspect(1.0)
		ax.legend()

		axes[0].set_ylabel(r"$\log_{10} \left( M_{\rm BH} \,/\, M_\odot \right)$")

		if args.plot_zou2024:
		print("Plotting Zou+2024")
		for z in [0.0, 2.0, 4.0]:
		log_mbh = mbh.get_log_mbh_zou2024(log_mstar, np.full_like(log_mstar, z), is_tng300=True)
		axes[0].plot(log_mstar, log_mbh, color=cmap2(norm(z)), lw=lw1, label=f"$z = {z:d}$ (Z24)" % z)
		axes[0].legend()

		fig.tight_layout()
		fig.savefig(args.savefig)


		if __name__ == "__main__":
		main()