document zou

# Email: akke.viitanen@helsinki.fi

# Date: 2024-06-26 18:14:36

"""

Implements the equations in Zou+2024

r"""

Implements Zou+2024 $p(\lambda_\mathrm{SAR} | M_\mathrm{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_\mathrm{SAR} \equiv

L_\mathrm{X} / M_\mathrm{star}$.

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_\mathrm{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_\mathrm{star} > 10^{9.5}\,M_\odot$

References

----------

https://ui.adsabs.harvard.edu/abs/2024ApJ...964..183Z/abstract

"""

import matplotlib as mpl

from lsst_inaf_agile import ueda2014

# NOTE: refer to Zou+2024 dataset readme on the grid parameters

Z_BOTTOM = 0.05

Z_TOP = 4.0

GRIDSIZE_Z = 51

LOG1PZGRID_BOUND = np.linspace(np.log10(1.0 + Z_BOTTOM), np.log10(1.0 + Z_TOP), GRIDSIZE_Z)

LOG1PZGRID = (LOG1PZGRID_BOUND[:-1] + LOG1PZGRID_BOUND[1:]) / 2.0

FUN_ZOU2024 = {}

LN10 = np.log(10)

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)

def _get_log_plambda(loglambda, logA, loglambdac, gamma1, gamma2):

    r"""Return log10 of $p(\lambda) \equiv A \times {(lambda / lambdac)}^{-gamma}."""

    gamma = _get_gamma(loglambda, loglambdac, gamma1, gamma2)

    return logA - gamma * (loglambda - loglambdac)

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)

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

    return loglambdac + np.where(log_Plambda > log_Plambda_c, rhs1, rhs2)

def _get_parameters(

    log_mstar,

    z,

    t,

    log_mstar_lim=(-np.inf, np.inf),

    z_lim=(-np.inf, np.inf),

):

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

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)

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)

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)

LN10 = np.log(10)

def get_schechter(logM, logMstar, alpha1, phi1, alpha2=np.nan, phi2=np.nan, factor=LN10):

    """Eq. 8 from Weaver+2020"""

    """

    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)

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],

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(

        [

    *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)

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

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)

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 by

    combining the accretion rate distribution of Zou+2024 (CTN AGN) and the

    CTK AGN fraction of Ueda+2014. The p_CTK is defined as:

    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)

    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))

def get_log_plambda_total(loglambda, mstar, z, t, test_integral=True, *args, **kwargs):

    """

    Get the total log plambda

    """

    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)

    *args,

    **kwargs,

):

    """

    Return log_lambda_sar sampled from the p(lambda)

    """

    """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)

# @cache

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:

# @cache

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

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

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

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")

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