#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2024-06-26 18:14:36
r"""
Implements Zou+2024 p(lambda_SAR | M_star, z).
The specific accretion rate is defined as the ratio of the 2-10 keV intrinsic
X-ray luminosity and the host galaxy stellar mass lambda_SAR = LX / Mstar.
Zou+2024 model the p(lambda_SAR | Mstar, z) as a piecewise double power law
function, where the power law index depends on the critical lambda_SAR,c.
Their observed sample is based on 8000 X-ray AGN and 1.3M normal galaxies
compiled from the CANDELS, LSST DDFs, and eFEDS. Nominally the sample is valid
for z < 4 and M_star > 1e9.5 Msun.
References
----------
https://ui.adsabs.harvard.edu/abs/2024ApJ...964..183Z/abstract
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from scipy import interpolate
from scipy.optimize import minimize_scalar
from scipy.special import gamma, gammaincc
from lsst_inaf_agile import ueda2014
# NOTE: refer to Zou+2024 dataset readme on the grid parameters
[docs]
LOGMGRID_BOUND = np.linspace(LOGMSTAR_BOTTOM, LOGMSTAR_TOP, GRIDSIZE_MSTAR)
[docs]
LOGMGRID = (LOGMGRID_BOUND[:-1] + LOGMGRID_BOUND[1:]) / 2.0
[docs]
LOG1PZGRID_BOUND = np.linspace(np.log10(1.0 + Z_BOTTOM), np.log10(1.0 + Z_TOP), GRIDSIZE_Z)
[docs]
LOG1PZGRID = (LOG1PZGRID_BOUND[:-1] + LOG1PZGRID_BOUND[1:]) / 2.0
[docs]
def _get_gamma(loglambda, loglambdac, gamma1, gamma2):
r"""
Return power law index gamma.
Parameters
----------
loglambda: float or array_like
log10 of the specific accretion rate.
loglambdac: float
log10 of the critical specific accretion rate.
gamma1:
gamma, where loglambda < loglambdac
gamma2:
gamma, where loglambda < loglambdac
Returns
-------
gamma: float
power law index defined as gamma1 if loglambda < loglambdac else gamma2.
"""
return np.where(loglambda < loglambdac, gamma1, gamma2)
[docs]
def _get_log_plambda(loglambda, logA, loglambdac, gamma1, gamma2):
r"""Return log10 of p(lambda) = A * (lambda / lambda_c)^gamma."""
gamma = _get_gamma(loglambda, loglambdac, gamma1, gamma2)
return logA - gamma * (loglambda - loglambdac)
[docs]
def _get_log_Plambda(loglambda, logA, loglambdac, gamma1, gamma2):
r"""Return analytic integral of $p(\lambda)$ up to specific accretoin rate $\lambda$."""
def fun(x):
return 10 ** _get_log_plambda(x, logA, loglambdac, gamma1, gamma2)
select = loglambda < loglambdac
I1 = np.where(select, (fun(loglambda) - fun(loglambdac)) / gamma1, 0.0)
I2 = np.where(select, fun(loglambdac), fun(loglambda)) / gamma2
return np.log10((I1 + I2) / np.log(10))
[docs]
def _get_inv_log_Plambda(log_Plambda, logA, loglambdac, gamma1, gamma2):
r"""Return inverse of the analytic integral of $p(\lambda)$ i.e. $P^{-1}$."""
assert np.all(log_Plambda <= 0)
part = np.log(10) * 10 ** (log_Plambda - logA)
rhs1 = -np.ma.log10(gamma1 * part - gamma1 / gamma2 + 1) / gamma1
rhs2 = -np.ma.log10(gamma2 * part) / gamma2
log_Plambda_c = _get_log_Plambda(loglambdac, logA, loglambdac, gamma1, gamma2)
return loglambdac + np.where(log_Plambda > log_Plambda_c, rhs1, rhs2)
[docs]
def _get_parameters(log_mstar, z, t, log_mstar_lim=(-np.inf, np.inf), z_lim=(-np.inf, np.inf)):
"""
Return Zou+2024 parameters at the given physical parameters.
Parameters
----------
log_mstar: float
log10 of the galaxy stellar mass
z: float
redshift
t: str
galaxy type, one of "main", "starforming", or "quiescent"
log_mstar_lim: tuple[float, float]
stellar mass limits for controlling boundary extrapolation
z_lim: tuple[float, float]
redshift limits for controlling boundary extrapolation
Returns
-------
(A, loglambdac, gamma1, gamma2): tuple[float, float, float, float]
Zou+2024 parameters.
"""
dirname = {"all": "main", "star-forming": "starforming", "quiescent": "quiescent"}.get(t)
assert dirname
# Clip stellar mass and redshift to some range to prevent extrapolation
log_mstar = np.clip(log_mstar, *log_mstar_lim)
z = np.clip(z, *z_lim)
def _get_parameter(parameter):
filename = f"data/plambda/zou2024/{dirname}/maps_{parameter}.fits"
if filename not in FUN_ZOU2024:
m = fits.open(filename)[0].data
FUN_ZOU2024[filename] = interpolate.RegularGridInterpolator(
(LOGMGRID, LOG1PZGRID), m.T, bounds_error=False, fill_value=None
)
return FUN_ZOU2024[filename]((log_mstar, np.log10(1 + z)))
parameters = "logA", "loglambdac", "gamma1", "gamma2"
return [_get_parameter(p) for p in parameters]
[docs]
def get_log_plambda(loglambda, log_mstar, z, t, *args, **kwargs):
r"""Return log10 of $p(\lambda)$ for the given physical parameters."""
parameters = _get_parameters(log_mstar, z, t, *args, **kwargs)
loglambda_min = _get_inv_log_Plambda(0.0, *parameters)
return np.where(loglambda > loglambda_min, _get_log_plambda(loglambda, *parameters), -np.inf)
[docs]
def get_log_Plambda(loglambda, log_mstar, z, t):
r"""Return log10 of $P(\lambda)$ for the given physical parameters."""
parameters = _get_parameters(log_mstar, z, t)
return _get_log_Plambda(loglambda, *parameters)
[docs]
def get_inv_log_Plambda(log_Plambda, log_mstar, z, t):
r"""Return log10 of $P^{-1}(\lambda)$ for the given physical parameters."""
if np.any(log_Plambda > 0):
raise ValueError
parameters = _get_parameters(log_mstar, z, t)
return _get_inv_log_Plambda(log_Plambda, *parameters)
[docs]
def get_schechter(logM, logMstar, alpha1, phi1, alpha2=np.nan, phi2=np.nan, factor=LN10):
"""
Return a (double) Schechter function.
Implements Eq. 8 from Weaver+2020.
"""
def fun(alpha, phi):
return np.where(np.isfinite(alpha * phi), phi * (10 ** (logM - logMstar)) ** (alpha + 1), 0.0)
return factor * np.exp(-(10 ** (logM - logMstar))) * (fun(alpha1, phi1) + fun(alpha2, phi2))
[docs]
def get_xray_luminosity_function_analytic(log_lx, z, log_mstar_min, log_mstar_max, dlog_mstar):
"""Return the analytic form of the X-ray luminosity function."""
p = np.array(
[
[10.831, 0.153, -0.033],
[-0.579, 0.048, 0.022],
[-1.489, -0.087, 0.016],
[-2.312, -0.658, 0.016],
[-3.326, -0.158, -0.002],
]
)
log_mstar_star = p[0, 0] + p[0, 1] * z + p[0, 2] * z**2
alpha1 = p[1, 0] + p[1, 1] * z + p[1, 2] * z**2
alpha2 = p[2, 0] + p[2, 1] * z + p[2, 2] * z**2
log_phi1 = p[3, 0] + p[3, 1] * z + p[3, 2] * z**2
log_phi2 = p[4, 0] + p[4, 1] * z + p[4, 2] * z**2
log_mstar = np.arange(log_mstar_min, log_mstar_max, dlog_mstar)
logA, loglambdac, gamma1, gamma2 = _get_parameters(log_mstar + dlog_mstar / 2, z, "all")
loglambda = log_lx - log_mstar - dlog_mstar / 2
my_gamma = np.where(loglambda < loglambdac, gamma1, gamma2)
log_plambda_star = get_log_plambda(log_lx - log_mstar_star, log_mstar + dlog_mstar / 2, z, "all")
def get(log_mstar, log_phi, alpha):
z = alpha + my_gamma + 1
assert np.all(z > 0)
g = gamma(z) * gammaincc(z, 10 ** (log_mstar - log_mstar_star))
return 10 ** (log_phi + log_plambda_star) * g
fun1 = get(log_mstar, log_phi1, alpha1) + get(log_mstar, log_phi2, alpha2)
fun2 = get(log_mstar + dlog_mstar, log_phi1, alpha1) + get(log_mstar + dlog_mstar, log_phi2, alpha2)
return np.sum(fun1 - fun2)
[docs]
def get_stellar_mass_function_wright2018(log_mstar, z, t="all"):
"""Return Wright+2018 stellar mass function."""
# NOTE: A(z) = A0 + A1 * z + A2 * z ** 2, see Wright+2018 table A3
p = np.array(
[
[10.831, 0.153, -0.033],
[-0.579, 0.048, 0.022],
[-1.489, -0.087, 0.016],
[-2.312, -0.658, 0.016],
[-3.326, -0.158, -0.002],
]
)
mstar = p[0, 0] + p[0, 1] * z + p[0, 2] * z**2
alpha1 = p[1, 0] + p[1, 1] * z + p[1, 2] * z**2
alpha2 = p[2, 0] + p[2, 1] * z + p[2, 2] * z**2
phi1 = p[3, 0] + p[3, 1] * z + p[3, 2] * z**2
phi2 = p[4, 0] + p[4, 1] * z + p[4, 2] * z**2
parameters = mstar, alpha1, 10**phi1, alpha2, 10**phi2
return get_schechter(log_mstar, *parameters)
[docs]
def get_xray_luminosity_function(
log_lx,
z,
fun_phi_star=get_stellar_mass_function_wright2018,
t="all",
log_mstar_min=8.0,
log_mstar_max=13.0,
dlog_mstar=0.01,
log_lambda_sar_min=-np.inf,
log_lambda_sar_max=+np.inf,
check_plambda_ctk=False,
*args,
**kwargs,
):
r"""
Return the X-ray luminosity function implied by $p(\lambda)$.
The X-ray luminosity function may be derived as the product of the galaxy
stellar mass function and $p(\lambda)$.
"""
log_lx = np.atleast_2d(log_lx)
z = np.atleast_2d(z)
log_mstar = np.arange(log_mstar_min, log_mstar_max + 1e-9, dlog_mstar)
# Estimate loglambda
loglambda = log_lx[:, :, None] - log_mstar[None, None, :]
# Estimate plambda
y1 = 10 ** get_log_plambda(loglambda, log_mstar[None, None, :], z[:, :, None], t, *args, **kwargs)
if check_plambda_ctk:
get_log_plambda_ctk(np.unique(loglambda), log_mstar, z, t)
# Assert min/max values of loglambda in p(lambda)
select = (log_lambda_sar_min <= loglambda) * (loglambda < log_lambda_sar_max)
y1 = np.where(select, y1, 0.0)
# Estimate phi_star
_m, _z = np.meshgrid(log_mstar, z[0, :, None])
y2 = fun_phi_star(_m, _z, t)
return np.squeeze(np.trapezoid(y1 * y2, log_mstar, axis=-1))
[docs]
def get_specific_accretion_rate_distribution_function(loglambda, log_mstar, z, t="all", dlog_mstar=0.10):
"""
Return the specific accretion rate distribution function.
The specific accretion rate distribution is the product of the galaxy
stellar mass function and the accretion rate probability.
"""
plambda = 10 ** get_log_plambda(loglambda, log_mstar, z, t)
phi = get_stellar_mass_function_wright2018(log_mstar, z, t)
return plambda * phi * dlog_mstar
[docs]
def get_frac_ctk_agn(log_lx, z):
"""Return the fraction of CTK AGN from Ueda+2014."""
frac = [ueda2014.get_f(log_lx, z, n) for n in [20, 21, 22, 23, 24, 25]]
frac_ctk_agn = np.sum(frac[-2:], axis=0) / np.sum(frac, axis=0)
assert np.all(frac_ctk_agn <= 1.0)
return frac_ctk_agn
[docs]
def get_log_plambda_ctk(loglambda, m, z, t, test_integral=True, *args, **kwargs):
r"""
Return the accretion rate distribution of the CTK AGN population.
Combines the accretion rate distribution of Zou+2024 (CTN AGN) and the CTK
AGN fraction of Ueda+2014. The p_CTK is defined as:
(1) p_tot = p_ctn + p_ctk
(2) p_tot = p_ctn / (1 - f_CTK_AGN)
so that
(3) p_CTK = f_CTK_AGN / (1 - f_CTK_AGN) * p_CTN
where the CTK AGN fraction is defined as
(4) f_CTK_AGN \equiv N_CTK / (N_CTN + N_CTK)
and p_CTN is the accretion rate distribution of CTN AGN.
"""
log_p_ctn = get_log_plambda(loglambda, m, z, t, *args, **kwargs)
frac_ctk_agn = get_frac_ctk_agn(loglambda + m, z)
log_p_ctk = log_p_ctn + np.log10(frac_ctk_agn / (1 - frac_ctk_agn))
# NOTE: safety assertion that p_agn is still a probability measure
# within l \in loglambda
if test_integral:
test = np.trapezoid(10**log_p_ctn + 10**log_p_ctk, loglambda)
assert np.all(test <= 1.00), (test, m, z, t)
return log_p_ctk
[docs]
def get_log_plambda_total(loglambda, mstar, z, t, test_integral=True, *args, **kwargs):
r"""Return the total (CTN + CTK) $\log p(\lambda)$."""
log_plambda_ctn = get_log_plambda(loglambda, mstar, z, t, *args, **kwargs)
log_plambda_ctk = get_log_plambda_ctk(loglambda, mstar, z, t, test_integral, *args, **kwargs)
log_plambda_tot = np.log10(10**log_plambda_ctn + 10**log_plambda_ctk)
return log_plambda_tot
# @cache
[docs]
def _get_log_lambda_SAR_min_Y(m, z, t, add_ctk, dlog_lambda, *args, **kwargs):
log_lambda_min = get_log_lambda_min(m, z, t, add_ctk, *args, **kwargs)
xbins = np.arange(log_lambda_min, 36, dlog_lambda)
xcens = xbins[:-1] + dlog_lambda / 2.0
y_ctn = 10 ** get_log_plambda(xcens, m, z, t, *args, **kwargs)
Y_ctn = np.append([0], np.cumsum(y_ctn))
y_ctk = 10 ** get_log_plambda_ctk(xcens, m, z, t, False, *args, **kwargs)
Y_ctk = np.append([0], np.cumsum(y_ctk))
norm = Y_ctn[-1] + Y_ctk[-1]
Y_ctn /= norm
Y_ctk /= norm
assert np.allclose(Y_ctn[-1] + Y_ctk[-1], 1.0)
return xbins, Y_ctn, Y_ctk
[docs]
def get_log_lambda_SAR(
m,
z,
t,
add_ctk=True,
dlog_lambda=0.001,
Nsample=1,
*args,
**kwargs,
):
"""Return log_lambda_sar sampled from the p(lambda)."""
xbins, Y_ctn, Y_ctk = _get_log_lambda_SAR_min_Y(m, z, t, add_ctk, dlog_lambda, *args, **kwargs)
U = np.random.rand(Nsample)
ret = np.where(
Y_ctn.max() > U,
+np.interp(U, Y_ctn, xbins),
-np.interp(U, Y_ctn.max() + Y_ctk, xbins),
)
assert np.all(np.isfinite(ret))
assert np.all(ret > -36.0)
assert np.all(ret < +36.0)
return np.squeeze(ret)
[docs]
X = np.linspace(20, 40, 2001)
# @cache
[docs]
def get_fraction_agn(xmin, mstar, z, t, add_ctk=True, *args, **kwargs):
r"""
Return the AGN fraction given the physical parameters.
The AGN fraction is defined as the integral of $p(\lambda)$ from $\lambda =
10^{32}$ (erg/s/Msun).
"""
def get_p_tot(x):
p_tot = 10 ** get_log_plambda(x, mstar, z, t, *args, **kwargs)
if add_ctk:
p_tot += 10 ** get_log_plambda_ctk(x, mstar, z, t, False, *args, **kwargs)
return p_tot
x = X[xmin <= X]
# NOTE: attempt to fix a possible shape mismatch between loglambda and mstar
try:
y = get_p_tot(x)
ret = np.trapezoid(y, x)
except ValueError:
y = get_p_tot(x[:, None])
ret = np.sum(y * dX, axis=0)
return ret
# @cache
[docs]
def get_log_lambda_min(mstar, z, t, *args, **kwargs):
r"""Return $\lambda_\mathrm{min}$ at which $p(\lambda)$ integrates to unity."""
def fun(xmin):
return (get_fraction_agn(xmin + 28, mstar, z, t, *args, **kwargs) - 1.0) ** 2
# NOTE: minimization starts at 0 by default. The offset speeds up things
return minimize_scalar(fun).x + 28
[docs]
def plot_fig10():
"""Plot Fig. 10 from Zou+2024."""
# Produce p(lambda)-like fig. 10 of Zou+2024
zs = 0.1, 0.3, 0.5, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0
ms = 9.75, 10.25, 10.75, 11.25, 11.75
fig, axes = plt.subplots(3, 3, figsize=(3 * 6.4, 3 * 4.8))
lvec = np.linspace(26.0, 36.0, 101)
for z, ax in zip(zs, axes.flatten(), strict=False):
for i, m in enumerate(ms):
pq = get_log_plambda(lvec, m, z, "quiescent")
ps = get_log_plambda(lvec, m, z, "star-forming")
ax.plot(lvec, pq, color=f"C{i}", linestyle="dotted")
ax.plot(lvec, ps, color=f"C{i}", linestyle="solid")
ax.plot([], [], color=f"C{i}", label=f"logM = {m}")
ax.plot([], [], color="k", linestyle="solid", label="star-forming")
ax.plot([], [], color="k", linestyle="dotted", label="quiescent")
ax.set_title(f"z = {z:.1f}")
ax.set_xlabel(r"$\log \lambda$ [erg/s/Msun]")
ax.set_ylabel(r"$p(\lambda)$ [1/dex]")
ax.legend()
plt.savefig("fig/zou2024_plambda_quiescent_starforming.pdf")
M = 9.5
z = 0.50
t = "all"
U = np.random.rand(10**6)
sample = get_inv_log_Plambda(np.log10(U), M, z, t)
plt.hist(sample, bins=np.linspace(28, 36, 801), density=True)
loglambda = np.linspace(28, 36, 801)
plt.plot(loglambda, 10 ** get_log_plambda(loglambda, M, z, t))
plt.semilogy()
plt.savefig("fig/zou2024_plambda_fig10.pdf")
[docs]
def plot_z_mstar_lambda_min():
r"""
Plot the $\lambda_\mathrm{min}$ in terms of $z$ and $M_\mathrm{star}$.
The minimum accretion rate is defined as $\lambda_\mathrm{min}$ at which
$p(\lambda)$ integrates to unity. In other words, at this accretion rate
all galaxies are considered to host accretion events.
"""
fig, axes = plt.subplots(2, 2, figsize=(2 * 6.4, 2 * 4.8))
log_z1, m = np.meshgrid(np.log10(1 + np.arange(0.20, 5.50, 0.01)), np.arange(8.50, 12.0, 0.01))
z = 10**log_z1 - 1
log_lambda_min_sf = get_inv_log_Plambda(np.log10(0.50), m, z, t="starforming")
log_lambda_min_qu = get_inv_log_Plambda(np.log10(0.50), m, z, t="quiescent")
P_sf = 10 ** get_log_Plambda(31.0, m, z, t="starforming")
P_qu = 10 ** get_log_Plambda(31.0, m, z, t="quiescent")
norm1 = mpl.colors.Normalize(vmin=31.0 - 0.20, vmax=31.0 + 0.20)
norm2 = mpl.colors.Normalize(vmin=0.50 - 0.20, vmax=0.50 + 0.20)
cmap = mpl.cm.coolwarm
kwargs = dict(extent=(np.log10(1 + 0.20), np.log10(1 + 5.50), 8.50, 12.0), cmap=cmap, aspect="auto")
axes[0, 0].imshow(log_lambda_min_sf, norm=norm1, **kwargs)
im2 = axes[0, 1].imshow(log_lambda_min_qu, norm=norm1, **kwargs)
axes[1, 0].imshow(P_sf, norm=norm2, **kwargs)
im4 = axes[1, 1].imshow(P_qu, norm=norm2, **kwargs)
axes[0, 0].set_title("star-forming")
axes[0, 1].set_title("quiescent")
cbar = plt.colorbar(im2, ax=axes[0, 1])
cbar.set_label("log_lambda_min", fontsize="large")
cbar = plt.colorbar(im4, ax=axes[1, 1])
cbar.set_label("frac_ctn_agn(loglambda=31.0)", fontsize="large")
for ax in axes.flatten():
ax.set_xlabel("log(1 + z)")
ax.set_ylabel("logMstar [Msun]")
plt.savefig("fig/zou2024_z_mstar_lambda_min.pdf")
[docs]
def plot_plambda_ctn_ctk():
r"""Plot $p(\lambda)$ separately for CTN and CTK AGN."""
fig, ax = plt.subplots(1, 1)
ax.set_title(r"logMstar = 10.50, z = 1.0", fontsize="large")
ax.set_xlabel(r"$\log \lambda_{\rm SAR}$ [erg/s/Msun]", fontsize="large")
ax.set_ylabel(r"$\log p(\lambda_{\rm SAR})$ [1/dex]", fontsize="large")
log_lambda = np.arange(31.5, 36, 0.01)
log_plambda = get_log_plambda(log_lambda, 10.50, 1.0, "star-forming")
plt.plot(log_lambda, log_plambda, color="purple", label="star-forming CTN", lw=2)
log_plambda = get_log_plambda(log_lambda, 10.50, 1.0, "quiescent")
plt.plot(log_lambda, log_plambda, color="orange", label="quiescent CTN", lw=2)
log_plambda_ctk = get_log_plambda_ctk(log_lambda, 10.50, 1.0, "star-forming")
plt.plot(log_lambda, log_plambda_ctk, linestyle="dotted", color="purple", label="star-forming CTK", lw=2)
log_plambda_ctk = get_log_plambda_ctk(log_lambda, 10.50, 1.0, "quiescent")
plt.plot(log_lambda, log_plambda_ctk, linestyle="dotted", color="orange", label="quiescent CTK", lw=2)
plt.legend(fontsize="large")
plt.savefig("fig/plambda_qu_sf_ctn_ctk.pdf")
[docs]
def plot_xray_luminosity_function():
"""Plot the X-ray luminosity function."""
fig, axes = plt.subplots(3, 4, figsize=(4 * 6.4, 3 * 4.8))
zbins = 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.4, 1.8, 2.2, 2.7, 4.5, 5.0
log_lx = np.arange(41.5, 46.5, 0.01)
for z, ax in zip(zbins, axes.flatten(), strict=False):
# Plot Zou XLF
log_xlf = np.log10(get_xray_luminosity_function(log_lx, z, log_mstar_min=8.5))
ax.plot(log_lx, log_xlf, label="Zou+ 2024", linestyle="solid")
# Plot Zou XLF with cut at logM>9.5
log_xlf = np.log10(get_xray_luminosity_function(log_lx, z, log_mstar_min=9.5))
ax.plot(log_lx, log_xlf, label="Zou+ 2024 cut logM>9.5", linestyle="dashed")
# Plot Zou XLF with cut clipping Mstar logM>9.5
log_xlf = np.log10(get_xray_luminosity_function(log_lx, z, log_mstar_lim=(9.5, 12.0)))
ax.plot(log_lx, log_xlf, label="Zou+ 2024 clip logM=[9.5, 12]", linestyle="dotted")
if z > 4:
# Plot Zou XLF with cut clipping Mstar and z
log_xlf = np.log10(
get_xray_luminosity_function(log_lx, z, log_mstar_lim=(9.5, 12.0), z_lim=(0, 4.0))
)
ax.plot(log_lx, log_xlf, label="Zou+ 2024 clip logM=[9.5, 12], z=[0.0, 4.0]", linestyle="dotted")
import ueda2014
x, y, y_lo, y_hi = ueda2014.get_xray_luminosity_function_data(z, True)
print(x, y, y_lo, y_hi)
ax.plot(x, np.log10(y), ".", label="Ueda+ 2014")
ax.set_title(f"z = {z}")
ax.set_xlabel(r"$\log L_{\rm X}$ [erg/s]")
ax.set_ylabel(r"$\log \Phi_{\rm X}$ [1/Mpc3/dex]")
ax.legend()
fig.savefig("fig/xlf_zou2024.pdf")
if __name__ == "__main__":
test_extrapolation_accuracy()