"""Zalesky et al. (2025) galaxy stellar mass function: models and digitized data.
Implements the single and double Schechter SMF in log10-mass units (Eqs. 9–10)
with the Eddington bias convolution kernel (Eq. 11) from the *Euclid* Cosmic
Dawn Survey stellar mass function paper. Digitized Tables 1, A.1–A.3 are
provided.
Reference
---------
Zalesky, L. et al. (Euclid Collaboration) 2025, A&A (submitted),
arXiv:2504.17867.
Cosmology used in the paper
---------------------------
H0 = 70 km/s/Mpc, Ωm = 0.3, ΩΛ = 0.7, Chabrier (2003) IMF.
"""
from __future__ import annotations
import warnings
warnings.warn(
"sum_stat.lf_smf.zalesky25 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 schechter_mass(
log10_mstar: jnp.ndarray,
log10_phi_star: float,
log10_mstar_char: float,
alpha: float,
) -> jnp.ndarray:
"""Single Schechter SMF in log10-mass (per-dex) units.
Evaluates (Zalesky et al. 2025 Eq. 9)::
φ(m) = ln(10) Φ* 10^{(α+1)(m−m*)} exp(−10^{m−m*})
Parameters
----------
log10_mstar : jnp.ndarray
log10(M / M_sun) at which to evaluate φ.
log10_phi_star : float
log10(Φ* / [Mpc^-3 dex^-1]).
log10_mstar_char : float
log10(M* / M_sun) — characteristic stellar mass.
alpha : float
Faint-end slope.
Returns
-------
phi : jnp.ndarray
φ(log10 M) [Mpc^-3 dex^-1], same shape as ``log10_mstar``.
References
----------
Schechter (1976), ApJ 203, 297.
Zalesky et al. (2025), arXiv:2504.17867, Eq. (9).
"""
phi_star = 10.0**log10_phi_star
x = 10.0 ** (log10_mstar - log10_mstar_char)
return jnp.log(10.0) * phi_star * x ** (alpha + 1.0) * jnp.exp(-x)
[docs]
@jax.jit
def double_schechter_mass(
log10_mstar: jnp.ndarray,
log10_mstar_char: float,
log10_phi1: float,
alpha1: float,
log10_phi2: float,
alpha2: float,
) -> jnp.ndarray:
"""Double Schechter SMF with shared characteristic mass.
Evaluates (Zalesky et al. 2025 Eq. 10)::
φ(m) = ln(10) exp(−x) [Φ1* x^{α1+1} + Φ2* x^{α2+1}]
where x = 10^{m − m*}.
Parameters
----------
log10_mstar : jnp.ndarray
log10(M / M_sun) at which to evaluate φ.
log10_mstar_char : float
log10(M* / M_sun) — shared characteristic stellar mass.
log10_phi1 : float
log10(Φ1* / [Mpc^-3 dex^-1]).
alpha1 : float
Slope of the first component.
log10_phi2 : float
log10(Φ2* / [Mpc^-3 dex^-1]).
alpha2 : float
Slope of the second component.
Returns
-------
phi : jnp.ndarray
φ(log10 M) [Mpc^-3 dex^-1], same shape as ``log10_mstar``.
References
----------
Zalesky et al. (2025), arXiv:2504.17867, Eq. (10).
"""
phi1 = 10.0**log10_phi1
phi2 = 10.0**log10_phi2
x = 10.0 ** (log10_mstar - log10_mstar_char)
return jnp.log(10.0) * jnp.exp(-x) * (phi1 * x ** (alpha1 + 1.0) + phi2 * x ** (alpha2 + 1.0))
[docs]
@jax.jit
def eddington_kernel(
delta_m: jnp.ndarray,
z: float,
sigma_edd: float = 0.6,
tau_c: float = 0.05,
) -> jnp.ndarray:
"""Eddington bias convolution kernel.
Product of a Gaussian (photometric scatter) and a Lorentzian (redshift
uncertainty) as defined in Zalesky et al. (2025) Eq. (11)::
D(δm, z) = exp(−δm² / 2σ²_Edd) × [τ_Edd / (2π)] / [(τ_Edd/2)² + δm²]
where τ_Edd = τ_c (1 + z). The kernel is unnormalised; the overall
normalisation is absorbed by Φ* during fitting.
Parameters
----------
delta_m : jnp.ndarray
Mass offset δm = m_obs − m_true [dex].
z : float
Redshift used to compute τ_Edd = τ_c (1 + z).
sigma_edd : float
Gaussian standard deviation [dex]. Default 0.6.
tau_c : float
Lorentzian width parameter [dex]. Default 0.05.
Returns
-------
kernel : jnp.ndarray
Unnormalised kernel values, same shape as ``delta_m``.
References
----------
Zalesky et al. (2025), arXiv:2504.17867, Eq. (11).
"""
tau_edd = tau_c * (1.0 + z)
gauss = jnp.exp(-0.5 * delta_m**2 / sigma_edd**2)
lorentz = (tau_edd / (2.0 * jnp.pi)) / ((0.5 * tau_edd) ** 2 + delta_m**2)
return gauss * lorentz
[docs]
def convolve_smf_eddington(
phi_fn,
log10_mstar: jnp.ndarray,
z: float,
sigma_edd: float = 0.6,
tau_c: float = 0.05,
n_kernel: int = 201,
kernel_range: float = 3.0,
) -> jnp.ndarray:
"""Forward-model: convolve intrinsic SMF with the Eddington bias kernel.
Computes the observed (Eddington-broadened) SMF via::
φ_obs(m) = ∫ φ_true(m − δm) D(δm, z) dδm
using a uniform grid over δm ∈ [−kernel_range, +kernel_range].
Parameters
----------
phi_fn : callable
Intrinsic SMF: ``phi_fn(log10_mstar) → jnp.ndarray``.
Must accept a 1-D JAX array and return φ [Mpc^-3 dex^-1].
log10_mstar : jnp.ndarray, shape (N,)
log10(M / M_sun) grid at which to evaluate φ_obs.
z : float
Effective redshift (for Lorentzian width τ_Edd = τ_c (1+z)).
sigma_edd : float
Eddington kernel Gaussian width [dex]. Default 0.6.
tau_c : float
Eddington kernel Lorentzian width parameter [dex]. Default 0.05.
n_kernel : int
Number of quadrature points in δm. Default 201.
kernel_range : float
Half-width of δm integration range [dex]. Default 3.0 (±5σ_Edd).
Returns
-------
phi_obs : jnp.ndarray, shape (N,)
Observed SMF [Mpc^-3 dex^-1].
"""
dm_arr = jnp.linspace(-kernel_range, kernel_range, n_kernel)
d_dm = 2.0 * kernel_range / (n_kernel - 1)
kernel_vals = eddington_kernel(dm_arr, z, sigma_edd, tau_c)
# m_prime[i, k] = mstar[i] − dm[k] — shape (N, n_kernel)
m_prime = log10_mstar[:, None] - dm_arr[None, :]
phi_prime = phi_fn(m_prime.reshape(-1)).reshape(m_prime.shape)
return jnp.sum(phi_prime * kernel_vals[None, :], axis=1) * d_dm
# ---------------------------------------------------------------------------
# Digitized Table 1 — survey volume, mass limits, galaxy counts
# ---------------------------------------------------------------------------
# Columns:
# 0 z_lo
# 1 z_hi
# 2 volume_1e6_Mpc3 [10^6 Mpc^3]
# 3 log10_mstar_lim [log10(M_sun)]
# 4 n_gal (M > M_lim)
nan = np.nan
# fmt: off
ZALESKY2025_TABLE1 = np.array([
# z_lo z_hi V[1e6 Mpc3] logMlim N_gal
[0.2, 0.5, 6.37, 8.94, 81111],
[0.5, 0.8, 15.35, 9.08, 164318],
[0.8, 1.1, 23.35, 9.20, 140565],
[1.1, 1.5, 39.84, 9.33, 140558],
[1.5, 2.0, 57.48, 9.47, 122640],
[2.0, 2.5, 60.61, 9.58, 89622],
[2.5, 3.0, 60.51, 9.68, 68940],
[3.0, 3.5, 58.80, 9.77, 40216],
[3.5, 4.5, 109.89, 9.93, 22477],
[4.5, 5.5, 98.76, 10.05, 2182],
[5.5, 6.5, 88.36, 10.16, 219],
])
# fmt: on
#: Column name → column index in :data:`ZALESKY2025_TABLE1`.
ZALESKY2025_TABLE1_COLUMNS: dict[str, int] = {
"z_lo": 0,
"z_hi": 1,
"volume_1e6_Mpc3": 2,
"log10_mstar_lim": 3,
"n_gal": 4,
}
# ---------------------------------------------------------------------------
# Digitized Tables A.1–A.3 — Schechter parameters from MCMC
# ---------------------------------------------------------------------------
# Column layout (22 columns):
# 0 z_lo
# 1 z_hi
# 2 log10_mstar_med median posterior [dex]
# 3 log10_mstar_ehi 84th − median [dex]
# 4 log10_mstar_elo median − 16th [dex]
# 5 log10_mstar_map MAP [dex]
# 6 alpha1_med
# 7 alpha1_ehi
# 8 alpha1_elo
# 9 alpha1_map
# 10 log10_phi1_med [Mpc^-3 dex^-1]
# 11 log10_phi1_ehi
# 12 log10_phi1_elo
# 13 log10_phi1_map
# 14 alpha2_med (NaN → single Schechter)
# 15 alpha2_ehi
# 16 alpha2_elo
# 17 alpha2_map
# 18 log10_phi2_med (NaN → single Schechter)
# 19 log10_phi2_ehi
# 20 log10_phi2_elo
# 21 log10_phi2_map
#
# For quiescent at z > 0.8: alpha1/log10_phi1 are NaN; alpha2/log10_phi2
# carry the single Schechter parameters.
# fmt: off
_TOTAL_DATA = np.array([
# zlo zhi logM*med +e -e MAP a1 +e -e MAP logP1 +e -e MAP a2 +e -e MAP logP2 +e -e MAP
[0.2, 0.5, 10.92, 0.09, 0.09, 10.94, -1.42, 0.07, 0.08, -1.55, -3.10, 0.19, 0.16, -3.38, -0.61, 0.37, 0.32, -0.84, -2.88, 0.17, 0.17, -2.79],
[0.5, 0.8, 10.88, 0.04, 0.04, 10.87, -1.37, 0.05, 0.05, -1.40, -3.03, 0.12, 0.12, -3.09, -0.54, 0.27, 0.22, -0.56, -2.78, 0.06, 0.04, -2.74],
[0.8, 1.1, 10.86, 0.04, 0.05, 10.89, -1.46, 0.07, 0.08, -1.54, -3.30, 0.14, 0.10, -3.46, -0.35, 0.26, 0.20, -0.55, -2.98, 0.07, 0.07, -2.93],
[1.1, 1.5, 10.98, 0.04, 0.04, 10.99, -1.44, 0.09, 0.08, -1.55, -3.62, 0.22, 0.14, -3.84, -0.77, 0.19, 0.16, -0.86, -3.18, 0.08, 0.09, -3.12],
[1.5, 2.0, 11.15, 0.04, 0.04, 11.14, -1.39, 0.06, 0.10, -1.55, -3.70, 0.21, 0.21, -4.05, -1.01, 0.41, 0.25, -1.10, -3.77, 0.24, 0.06, -3.56],
[2.0, 2.5, 11.15, 0.02, 0.02, 11.15, -1.55, 0.00, 0.01, -1.55, -3.83, 0.03, 0.02, -3.84, nan, nan, nan, nan, nan, nan, nan, nan],
[2.5, 3.0, 11.04, 0.02, 0.02, 11.05, -1.70, 0.00, 0.01, -1.70, -3.93, 0.03, 0.02, -3.94, nan, nan, nan, nan, nan, nan, nan, nan],
[3.0, 3.5, 11.04, 0.02, 0.02, 11.05, -1.69, 0.01, 0.01, -1.70, -4.07, 0.03, 0.03, -4.09, nan, nan, nan, nan, nan, nan, nan, nan],
[3.5, 4.5, 10.87, 0.05, 0.06, 10.93, -2.01, 0.06, 0.03, -2.05, -4.48, 0.12, 0.07, -4.57, nan, nan, nan, nan, nan, nan, nan, nan],
[4.5, 5.5, 10.68, 0.24, 0.18, 10.95, -2.14, 0.10, 0.04, -2.20, -4.99, 0.49, 0.26, -5.46, nan, nan, nan, nan, nan, nan, nan, nan],
[5.5, 6.5, 11.03, 0.12, 0.16, 11.20, -2.05, 0.10, 0.16, -2.20, -5.75, 0.32, 0.15, -6.09, nan, nan, nan, nan, nan, nan, nan, nan],
])
_SF_DATA = np.array([
[0.2, 0.5, 10.79, 0.11, 0.13, 10.84, -1.42, 0.06, 0.06, -1.54, -3.08, 0.14, 0.14, -3.34, -0.49, 0.62, 0.63, -0.94, -3.26, 0.29, 0.29, -3.12],
[0.5, 0.8, 10.81, 0.06, 0.05, 10.80, -1.41, 0.04, 0.04, -1.46, -3.08, 0.09, 0.12, -3.16, -0.51, 0.31, 0.32, -0.56, -3.05, 0.07, 0.07, -2.97],
[0.8, 1.1, 10.79, 0.06, 0.06, 10.81, -1.49, 0.05, 0.05, -1.51, -3.29, 0.09, 0.09, -3.33, 0.01, 0.32, 0.27, -0.10, -3.25, 0.08, 0.08, -3.21],
[1.1, 1.5, 11.08, 0.04, 0.06, 11.07, -1.37, 0.09, 0.06, -1.33, -3.47, 0.09, 0.15, -3.38, -0.82, 1.00, 0.46, 0.97, -4.12, 1.00, 0.33, -5.02],
[1.5, 2.0, 11.13, 0.04, 0.04, 11.09, -1.36, 0.03, 0.05, -1.33, -3.58, 0.08, 0.13, -3.50, -0.44, 1.11, 0.79, 1.00, -4.47, 0.90, 0.28, -4.87],
[2.0, 2.5, 11.08, 0.03, 0.02, 11.09, -1.55, 0.00, 0.00, -1.55, -3.79, 0.02, 0.02, -3.80, nan, nan, nan, nan, nan, nan, nan, nan],
[2.5, 3.0, 11.01, 0.02, 0.02, 11.02, -1.70, 0.01, 0.00, -1.70, -3.90, 0.03, 0.02, -3.92, nan, nan, nan, nan, nan, nan, nan, nan],
])
# For quiescent at z > 0.8: alpha1/log10_phi1 = NaN; single Schechter uses
# alpha2/log10_phi2 columns (as in the published table).
_QUIESCENT_DATA = np.array([
[0.2, 0.5, 10.87, 0.05, 0.05, 10.91, -1.85, 0.29, 0.32, -2.34, -4.88, 1.18, 0.33, -5.87, -0.48, 0.14, 0.11, -0.59, -2.86, 0.05, 0.05, -2.88],
[0.5, 0.8, 10.83, 0.03, 0.03, 10.84, -1.66, 0.45, 0.46, -2.48, -5.00, 2.28, 0.37, -6.39, -0.39, 0.09, 0.00, -0.44, -2.87, 0.01, 0.01, -2.86],
[0.8, 1.1, 10.74, 0.02, 0.02, 10.74, nan, nan, nan, nan, nan, nan, nan, nan, -0.17, 0.04, 0.04, -0.18, -3.09, 0.02, 0.02, -3.09],
[1.1, 1.5, 10.62, 0.02, 0.03, 10.62, nan, nan, nan, nan, nan, nan, nan, nan, 0.46, 0.07, 0.07, 0.45, -3.31, 0.01, 0.01, -3.30],
[1.5, 2.0, 10.82, 0.02, 0.02, 10.84, nan, nan, nan, nan, nan, nan, nan, nan, 0.10, 0.11, 0.11, 0.02, -3.88, 0.01, 0.01, -3.89],
[2.0, 2.5, 10.95, 0.04, 0.05, 10.98, nan, nan, nan, nan, nan, nan, nan, nan, 0.05, 0.18, 0.05, 0.00, -4.62, 0.04, 0.04, -4.62],
[2.5, 3.0, 10.93, 0.17, 0.23, 11.12, nan, nan, nan, nan, nan, nan, nan, nan, 0.45, 0.85, 0.40, 0.00, -5.47, 0.90, 0.23, -5.49],
])
# fmt: on
#: Column name → column index in the :data:`ZALESKY2025` parameter arrays.
ZALESKY2025_COLUMNS: dict[str, int] = {
"z_lo": 0,
"z_hi": 1,
"log10_mstar_med": 2,
"log10_mstar_ehi": 3,
"log10_mstar_elo": 4,
"log10_mstar_map": 5,
"alpha1_med": 6,
"alpha1_ehi": 7,
"alpha1_elo": 8,
"alpha1_map": 9,
"log10_phi1_med": 10,
"log10_phi1_ehi": 11,
"log10_phi1_elo": 12,
"log10_phi1_map": 13,
"alpha2_med": 14,
"alpha2_ehi": 15,
"alpha2_elo": 16,
"alpha2_map": 17,
"log10_phi2_med": 18,
"log10_phi2_ehi": 19,
"log10_phi2_elo": 20,
"log10_phi2_map": 21,
}
#: Digitized Tables A.1–A.3 of Zalesky et al. (2025).
#:
#: Keys are ``'total'``, ``'star_forming'``, and ``'quiescent'``.
#: Each value is a numpy array of shape (n_zbins, 22).
#: Use :data:`ZALESKY2025_COLUMNS` to map column names to indices.
#: For quiescent at z > 0.8 the single Schechter component occupies
#: the ``alpha2``/``log10_phi2`` columns (``alpha1``/``log10_phi1`` are NaN).
ZALESKY2025: dict[str, np.ndarray] = {
"total": _TOTAL_DATA,
"star_forming": _SF_DATA,
"quiescent": _QUIESCENT_DATA,
}
[docs]
def zalesky2025_smf(
log10_mstar: np.ndarray | jnp.ndarray,
z_bin: int,
population: str = "total",
use_map: bool = False,
) -> jnp.ndarray:
"""Evaluate the Zalesky et al. (2025) best-fit SMF for a redshift bin.
Looks up the Schechter parameters from Tables A.1–A.3 and evaluates
the appropriate single or double Schechter model via
:func:`schechter_mass` or :func:`double_schechter_mass`.
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:`ZALESKY2025` for the bin ordering (Table A.1–A.3).
population : {'total', 'star_forming', 'quiescent'}
Galaxy population.
use_map : bool
If True use MAP (maximum a posteriori) values; otherwise median.
Returns
-------
phi : jnp.ndarray
φ(log10 M) [Mpc^-3 dex^-1].
Examples
--------
Evaluate the total SMF in the 0.2 < z ≤ 0.5 bin:
>>> import numpy as np
>>> log_m = np.linspace(9.0, 12.5, 200)
>>> phi = zalesky2025_smf(log_m, z_bin=0, population='total')
"""
if population not in ZALESKY2025:
raise ValueError(f"population must be one of {list(ZALESKY2025)!r}, got {population!r}")
data = ZALESKY2025[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}"
)
col = ZALESKY2025_COLUMNS
row = data[z_bin]
_ = 3 if use_map else 0 # offset unused: lookups use named keys in ZALESKY2025_COLUMNS
log10_mstar_char = row[col["log10_mstar_map"] if use_map else col["log10_mstar_med"]]
alpha1 = row[col["alpha1_map"] if use_map else col["alpha1_med"]]
log10_phi1 = row[col["log10_phi1_map"] if use_map else col["log10_phi1_med"]]
alpha2 = row[col["alpha2_map"] if use_map else col["alpha2_med"]]
log10_phi2 = row[col["log10_phi2_map"] if use_map else col["log10_phi2_med"]]
m = jnp.asarray(log10_mstar, dtype=float)
# Single Schechter: quiescent z > 0.8 uses alpha2/phi2; others use alpha1/phi1
if np.isnan(alpha1) and np.isnan(log10_phi1):
# quiescent single Schechter (alpha2/phi2 columns)
return schechter_mass(m, log10_phi2, log10_mstar_char, alpha2)
if np.isnan(alpha2) or np.isnan(log10_phi2):
# total / star-forming single Schechter (alpha1/phi1 columns)
return schechter_mass(m, log10_phi1, log10_mstar_char, alpha1)
# Double Schechter
return double_schechter_mass(m, log10_mstar_char, log10_phi1, alpha1, log10_phi2, alpha2)