"""Ilbert et al. (2013) galaxy stellar mass function: model and digitized data.
Implements the double Schechter stellar mass function (Eq. 2) and provides
the best-fit parameters digitized from Table 2 of Ilbert et al. (2013),
covering the full, quiescent, and star-forming galaxy populations over
0.2 < z < 4.
Reference
---------
Ilbert, O. et al. 2013, A&A, 556, A55.
DOI: 10.1051/0004-6361/201321100
Cosmology used in the paper
----------------------------
Ωm = 0.3, ΩΛ = 0.7, H0 = 70 km/s/Mpc, Chabrier (2003) IMF.
"""
from __future__ import annotations
import warnings
warnings.warn(
"sum_stat.lf_smf.ilbert2013 is deprecated and will be removed in a future release. "
"Schechter models and reference data have moved to the gga_model package.",
DeprecationWarning,
stacklevel=2,
)
import jax
import jax.numpy as jnp
import numpy as np
# ---------------------------------------------------------------------------
# Model functions
# ---------------------------------------------------------------------------
[docs]
@jax.jit
def double_schechter(
log10_mstar: jnp.ndarray,
log10_mstar_char: float,
phi1_star: float,
alpha1: float,
phi2_star: float,
alpha2: float,
) -> jnp.ndarray:
"""Double Schechter stellar mass function in log10-mass units.
Evaluates the double Schechter function of Ilbert et al. (2013) Eq. (2),
converting to per-dex units::
φ(log10 M) = ln(10) exp(−x) [φ1* x^(α1+1) + φ2* x^(α2+1)]
where x = M / M* = 10^(log10_M − log10_M*). For a single Schechter
function set ``phi2_star = 0.0``.
Parameters
----------
log10_mstar : jnp.ndarray
log10(M / M_sun) at which to evaluate φ.
log10_mstar_char : float
log10(M* / M_sun) — characteristic stellar mass.
phi1_star : float
φ1* [Mpc^-3] — normalisation of the first Schechter component.
alpha1 : float
α1 — slope of the first Schechter component.
phi2_star : float
φ2* [Mpc^-3] — normalisation of the second Schechter component.
Set to 0.0 for a single Schechter function.
alpha2 : float
α2 — slope of the second Schechter component.
Ignored when ``phi2_star = 0.0``.
Returns
-------
phi : jnp.ndarray
φ(log10 M) [Mpc^-3 dex^-1], same shape as ``log10_mstar``.
References
----------
Ilbert et al. (2013), A&A 556, A55, Eq. (2).
"""
x = 10.0 ** (log10_mstar - log10_mstar_char)
term1 = phi1_star * x ** (alpha1 + 1.0)
term2 = phi2_star * x ** (alpha2 + 1.0)
return jnp.log(10.0) * jnp.exp(-x) * (term1 + term2)
# ---------------------------------------------------------------------------
# Digitized Table 2 — best-fit parameters
# ---------------------------------------------------------------------------
# All φ* values are in Mpc^-3 (table values × 10^-3 converted here).
# Errors are stored as (err_hi, err_lo) pairs (positive values).
# NaN is used where a parameter is undefined (single Schechter bins).
#
# Column layout of each _DATA array:
# 0 z_lo
# 1 z_hi
# 2 log_m_complete [dex]
# 3 log_m_star [dex]
# 4 log_m_star_err_hi [dex]
# 5 log_m_star_err_lo [dex]
# 6 phi1_star [Mpc^-3]
# 7 phi1_star_err_hi [Mpc^-3]
# 8 phi1_star_err_lo [Mpc^-3]
# 9 alpha1
# 10 alpha1_err_hi
# 11 alpha1_err_lo
# 12 phi2_star [Mpc^-3] (NaN → single Schechter)
# 13 phi2_star_err_hi [Mpc^-3]
# 14 phi2_star_err_lo [Mpc^-3]
# 15 alpha2 (NaN → single Schechter)
# 16 alpha2_err_hi
# 17 alpha2_err_lo
# 18 log_rho_star [log10(M_sun Mpc^-3)]
# 19 log_rho_star_err_hi [dex]
# 20 log_rho_star_err_lo [dex]
nan = np.nan
# fmt: off
_FULL_DATA = np.array([
# z_lo z_hi logMc logM* +e -e phi1* +e -e a1 +e -e phi2* +e -e a2 +e -e logρ* +e -e
[0.2, 0.5, 7.93, 10.88, 0.10, 0.10, 1.68e-3, 0.61e-3, 0.61e-3, -0.69, 0.40, 0.36, 0.77e-3, 0.40e-3, 0.53e-3, -1.42, 0.07, 0.14, 8.308, 0.080, 0.095],
[0.5, 0.8, 8.70, 11.03, 0.08, 0.10, 1.22e-3, 0.31e-3, 0.39e-3, -1.00, 0.31, 0.31, 0.16e-3, 0.32e-3, 0.32e-3, -1.64, 0.20, 0.50, 8.226, 0.065, 0.073],
[0.8, 1.1, 9.13, 10.87, 0.06, 0.06, 2.03e-3, 0.27e-3, 0.32e-3, -0.52, 0.35, 0.27, 0.29e-3, 0.30e-3, 0.30e-3, -1.62, 0.19, 0.32, 8.251, 0.069, 0.073],
[1.1, 1.5, 9.42, 10.71, 0.08, 0.08, 1.35e-3, 0.34e-3, 0.35e-3, -0.08, 0.55, 0.52, 0.67e-3, 0.41e-3, 0.44e-3, -1.46, 0.16, 0.29, 8.086, 0.069, 0.071],
[1.5, 2.0, 9.67, 10.74, 0.07, 0.06, 0.88e-3, 0.10e-3, 0.12e-3, -0.24, 0.27, 0.28, 0.33e-3, 0.06e-3, 0.07e-3, -1.60, nan, nan, 7.909, 0.090, 0.072],
[2.0, 2.5,10.04, 10.74, 0.07, 0.07, 0.62e-3, 0.07e-3, 0.07e-3, -0.22, 0.29, 0.29, 0.15e-3, 0.04e-3, 0.04e-3, -1.60, nan, nan, 7.682, 0.110, 0.081],
[2.5, 3.0,10.24, 10.76, 0.16, 0.15, 0.26e-3, 0.05e-3, 0.08e-3, -0.15, 0.86, 0.68, 0.14e-3, 0.11e-3, 0.06e-3, -1.60, nan, nan, 7.489, 0.230, 0.123],
[3.0, 4.0,10.27, 10.74, 0.44, 0.20, 0.03e-3, 0.02e-3, 0.02e-3, 0.95, 1.05, 1.21, 0.09e-3, 0.07e-3, 0.07e-3, -1.60, nan, nan, 7.120, 0.234, 0.168],
])
# Quiescent: single Schechter at z > 0.5 (phi2_star = alpha2 = NaN)
_QUIESCENT_DATA = np.array([
# z_lo z_hi logMc logM* +e -e phi1* +e -e a1 +e -e phi2* +e -e a2 +e -e logρ* +e -e
[0.2, 0.5, 8.24, 10.91, 0.07, 0.08, 1.27e-3, 0.21e-3, 0.19e-3, -0.68, 0.22, 0.13, 0.03e-3, 0.07e-3, 0.07e-3, -1.52, 0.27, 0.44, 7.986, 0.087, 0.102],
[0.5, 0.8, 8.96, 10.93, 0.04, 0.04, 1.11e-3, 0.10e-3, 0.09e-3, -0.46, 0.05, 0.05, nan, nan, nan, nan, nan, nan, 7.920, 0.054, 0.058],
[0.8, 1.1, 9.37, 10.81, 0.03, 0.03, 1.57e-3, 0.09e-3, 0.09e-3, -0.11, 0.05, 0.05, nan, nan, nan, nan, nan, nan, 7.985, 0.044, 0.049],
[1.1, 1.5, 9.60, 10.72, 0.03, 0.03, 0.70e-3, 0.03e-3, 0.03e-3, 0.04, 0.06, 0.06, nan, nan, nan, nan, nan, nan, 7.576, 0.041, 0.046],
[1.5, 2.0, 9.87, 10.73, 0.03, 0.04, 0.22e-3, 0.01e-3, 0.01e-3, 0.10, 0.09, 0.09, nan, nan, nan, nan, nan, nan, 7.093, 0.049, 0.053],
[2.0, 2.5,10.11, 10.59, 0.06, 0.06, 0.10e-3, 0.01e-3, 0.01e-3, 0.88, 0.23, 0.21, nan, nan, nan, nan, nan, nan, 6.834, 0.076, 0.084],
[2.5, 3.0,10.39, 10.27, 0.10, 0.08, 0.003e-3,0.006e-3,0.002e-3, 3.26, 0.93, 0.93, nan, nan, nan, nan, nan, nan, 6.340, 0.079, 0.121],
])
# Star-forming: alpha2 fixed at -1.6 for z > 1.5 (errors set to NaN)
_SF_DATA = np.array([
# z_lo z_hi logMc logM* +e -e phi1* +e -e a1 +e -e phi2* +e -e a2 +e -e logρ* +e -e
[0.2, 0.5, 7.86, 10.60, 0.16, 0.11, 1.16e-3, 0.31e-3, 0.41e-3, 0.17, 0.57, 0.65, 1.08e-3, 0.29e-3, 0.31e-3, -1.40, 0.04, 0.04, 8.051, 0.091, 0.101],
[0.5, 0.8, 8.64, 10.62, 0.17, 0.10, 0.77e-3, 0.22e-3, 0.30e-3, 0.03, 0.58, 0.79, 0.84e-3, 0.28e-3, 0.35e-3, -1.43, 0.06, 0.09, 7.933, 0.069, 0.073],
[0.8, 1.1, 9.04, 10.80, 0.11, 0.12, 0.50e-3, 0.33e-3, 0.31e-3, -0.67, 0.75, 0.67, 0.48e-3, 0.41e-3, 0.41e-3, -1.51, 0.11, 0.67, 7.908, 0.065, 0.067],
[1.1, 1.5, 9.29, 10.67, 0.11, 0.09, 0.53e-3, 0.23e-3, 0.18e-3, 0.11, 0.61, 0.78, 0.87e-3, 0.29e-3, 0.40e-3, -1.37, 0.08, 0.15, 7.916, 0.059, 0.061],
[1.5, 2.0, 9.65, 10.66, 0.07, 0.06, 0.75e-3, 0.08e-3, 0.10e-3, -0.08, 0.28, 0.31, 0.39e-3, 0.07e-3, 0.07e-3, -1.60, nan, nan, 7.841, 0.097, 0.071],
[2.0, 2.5,10.01, 10.73, 0.08, 0.08, 0.50e-3, 0.07e-3, 0.07e-3, -0.33, 0.33, 0.33, 0.15e-3, 0.05e-3, 0.05e-3, -1.60, nan, nan, 7.614, 0.123, 0.084],
[2.5, 3.0,10.20, 10.90, 0.20, 0.21, 0.15e-3, 0.08e-3, 0.08e-3, -0.62, 1.04, 0.84, 0.11e-3, 0.07e-3, 0.07e-3, -1.60, nan, nan, 7.453, 0.193, 0.128],
[3.0, 4.0,10.26, 10.74, 0.29, 0.17, 0.02e-3, 0.01e-3, 0.01e-3, 1.31, 0.87, 0.87, 0.10e-3, 0.06e-3, 0.04e-3, -1.60, nan, nan, 7.105, 0.245, 0.170],
])
# fmt: on
# Column name → column index mapping for the data arrays
_COL = {
"z_lo": 0,
"z_hi": 1,
"log_m_complete": 2,
"log_m_star": 3,
"log_m_star_err_hi": 4,
"log_m_star_err_lo": 5,
"phi1_star": 6,
"phi1_star_err_hi": 7,
"phi1_star_err_lo": 8,
"alpha1": 9,
"alpha1_err_hi": 10,
"alpha1_err_lo": 11,
"phi2_star": 12,
"phi2_star_err_hi": 13,
"phi2_star_err_lo": 14,
"alpha2": 15,
"alpha2_err_hi": 16,
"alpha2_err_lo": 17,
"log_rho_star": 18,
"log_rho_star_err_hi": 19,
"log_rho_star_err_lo": 20,
}
#: Digitized Table 2 of Ilbert et al. (2013).
#:
#: Keys are ``'full'``, ``'quiescent'``, and ``'star_forming'``.
#: Each value is a numpy array of shape (n_zbins, 21).
#: Use :data:`ILBERT2013_COLUMNS` to map column names to indices.
ILBERT2013: dict[str, np.ndarray] = {
"full": _FULL_DATA,
"quiescent": _QUIESCENT_DATA,
"star_forming": _SF_DATA,
}
#: Column name → column index in the :data:`ILBERT2013` arrays.
ILBERT2013_COLUMNS: dict[str, int] = _COL
[docs]
def ilbert2013_smf(
log10_mstar: np.ndarray | jnp.ndarray,
z_bin: int,
population: str = "full",
) -> jnp.ndarray:
"""Evaluate the Ilbert et al. (2013) SMF for a given redshift bin.
Looks up the best-fit double Schechter parameters from Table 2 and
evaluates the model via :func:`double_schechter`.
Parameters
----------
log10_mstar : array-like
log10(M / M_sun) values at which to evaluate φ.
z_bin : int
Zero-based index into the redshift bins for the chosen population.
See :data:`ILBERT2013` for the bin ordering.
population : {'full', 'quiescent', 'star_forming'}
Galaxy population.
Returns
-------
phi : jnp.ndarray
φ(log10 M) [Mpc^-3 dex^-1].
Examples
--------
Evaluate the full SMF in the 0.2 < z < 0.5 bin:
>>> import numpy as np
>>> log_m = np.linspace(9.0, 12.5, 200)
>>> phi = ilbert2013_smf(log_m, z_bin=0, population='full')
"""
if population not in ILBERT2013:
raise ValueError(f"population must be one of {list(ILBERT2013)!r}, got {population!r}")
data = ILBERT2013[population]
if not (0 <= z_bin < len(data)):
raise ValueError(
f"z_bin must be in [0, {len(data) - 1}] for population '{population}', " f"got {z_bin}"
)
row = data[z_bin]
log10_mstar_char = row[_COL["log_m_star"]]
phi1_star = row[_COL["phi1_star"]]
alpha1 = row[_COL["alpha1"]]
phi2_star = row[_COL["phi2_star"]]
alpha2 = row[_COL["alpha2"]]
# NaN phi2_star → single Schechter (second term vanishes)
phi2_star = 0.0 if np.isnan(phi2_star) else phi2_star
alpha2 = 0.0 if np.isnan(alpha2) else alpha2 # value irrelevant when phi2_star=0
return double_schechter(
jnp.asarray(log10_mstar, dtype=float),
log10_mstar_char,
phi1_star,
alpha1,
phi2_star,
alpha2,
)