Accretion History & Rate#

This module implements a range of models for the growth of dark matter halos over cosmic time. It provides three related quantities, each through a wrapper function that selects between a number of fitting functions from the literature:

The formation redshift \(z_{\rm f}\), the redshift at which the main progenitor of a halo had assembled a fixed fraction (usually half) of its final mass. See formationRedshift().
The mass accretion history (MAH), \(M(z)\), the mass of the main progenitor as a function of redshift, anchored at a chosen reference redshift. See massAccretionHistory().
The mass accretion rate (MAR), the dimensionless rate of mass growth, evaluated either instantaneously or averaged over one dynamical time. See massAccretionRate().

The MAH models often depend on a formation redshift, and MAR models in turn may depend on a MAH.

Basics#

The three wrapper functions partially share a common interface. For example, the half-mass formation redshift of a halo is computed as:

from colossus.cosmology import cosmology
from colossus.halo import accretion_history as acc_hist

cosmology.setCosmology('planck18')
zf = acc_hist.formationRedshift(1E12, model = 'bera26')

The mass accretion history can be evaluated either by passing a formation redshift, by using a default model thereof, or by passing the name of another model:

z = np.linspace(0.0, 5.0, 50)
Mz_1 = acc_hist.massAccretionHistory(1E12, z, model = 'wechsler02', zf = 2.5)
Mz_2 = acc_hist.massAccretionHistory(1E12, z, model = 'wechsler02', zf = 'bera26')

Similarly, the mass accretion rate can be predicted directly by models or be derived from models for the accretion history:

gamma_1 = acc_hist.massAccretionRate(1E12, 0.0, model = 'bera26')
gamma_2 = acc_hist.massAccretionRate(1E12, 0.0, model = 'wechsler02', zf = 'correa15')

Each model was calibrated for one particular spherical overdensity mass definition (see Halo Mass Definitions for details). If the user passes a mass in a different definition through the mdef parameter, the wrapper converts it to the model’s native definition before evaluating the fit. If mdef is omitted, the mass is assumed to already be in the model’s native definition and no conversion is performed.

Note

The user must ensure that the cosmology is set consistently. Most models were calibrated for a particular cosmology, and neither the wrapper functions nor the individual model functions issue a warning if the set cosmology does not match the model cosmology.

Formation redshift models#

The following models can be passed as the model parameter to the formationRedshift() function:

ID	Default mdef	M-range (z=0)	z-range	Reference
lacey93	vir	Any	Any	Lacey & Cole 1993
vandenbosch02	vir	Any	Any	van den Bosch 2002
mcbride09	vir	1E11 < M < 1E14	0 < z < 6	McBride et al. 2009
giocoli12	vir	Any	Any	Giocoli et al. 2012
vandenbosch14	vir	1E11 < M < 1E14	Any	van den Bosch et al. 2014
correa15	vir	Any	Any	Correa et al. 2015
bera26	200m	Any	Any	Bera & Diemer 2026

The properties of these models are recorded in models_zf.

Mass accretion history models#

The following models can be passed as the model parameter to the massAccretionHistory() function. The z0 > 0? column indicates whether a model can be anchored at a non-zero observation redshift (otherwise it must be anchored at \(z = 0\)):

ID	z0 > 0?	M-range (z=0)	z-range	Reference
wechsler02	yes	M > 1E12	Any	Wechsler et al. 2002
vandenbosch02	no	Any	Any	van den Bosch 2002
mcbride09	yes	1E11 < M < 1E14	0 < z < 6	McBride et al. 2009
vandenbosch14	yes	1E11 < M < 1E14	Any	van den Bosch et al. 2014
correa15	yes	Any	Any	Correa et al. 2015

Most models rely on a formation redshift and are thus given a default formation redshift model. For example, the vandenbosch02 MAH model uses the vandenbosch02 formation redshift model by default, but other models can be chosen. The model properties are set in models_mah.

Mass accretion rate models#

The following models can be passed as the model parameter to the massAccretionRate() function:

ID	Default delta_t	inst	M-range (z=0)	z-range	Reference
wechsler02	crossing	yes	M > 1E12	Any	Wechsler et al. 2002
mcbride09	crossing	yes	1E11 < M < 1E14	0 < z < 6	McBride et al. 2009
vandenbosch14	crossing	yes	1E11 < M < 1E14	Any	van den Bosch et al. 2014
correa15	crossing	yes	Any	Any	Correa et al. 2015
rodriguezpuebla16	peri	yes	1E10 < M < 1E15	0 < z < 6	Rodriguez-Puebla et al. 2016
yung24	peri	no	–	z > 6	Yung et al. 2024
yung25	peri	no	–	2 < z < 18	Yung et al. 2025
bera26	crossing	no	Any	Any	Bera & Diemer 2026

The Default delta_t column gives the time interval used when delta_t is not set, and the inst column indicates whether the instantaneous rate (delta_t = 'inst') is available. For models whose instantaneous rate is not given in closed form, but where the mass accretion history is available (e.g. vandenbosch14), the MAR is computed as a numerical derivative of the corresponding mass accretion history.

The mass and redshift ranges quoted in the tables above are the approximate regimes over which each fit was calibrated; Any denotes the analytic excursion-set models, which can in principle be evaluated at any mass and redshift. Please consult the cited papers for the exact limits. The model properties are set in models_mar.

Module contents#

`FormationRedshiftModel`()	Characteristics of formation redshift models.
`models_zf`	List of formation redshift models.
`formationRedshift`(M[, z0, model, mdef, ...])	Wrapper function for formation redshift.
`zfModelLacey93`(M, z0[, f_form, C, zmax, ...])	The formation redshift model of Lacey & Cole 1993.
`zfModelVanDenBosch02`(M, z0[, f_form, C, ...])	The formation redshift model of van den Bosch 2002.
`zfModelMcBride09`(M)	The formation redshift model of McBride et al. 2009.
`zfModelGiocoli12`(M, z0[, f_form, zmax, ...])	The formation redshift model of Giocoli et al. 2012.
`zfModelVanDenBosch14`(M, z0[, f_form, ...])	The formation redshift model of van den Bosch et al. 2014.
`zfModelCorrea15`(M)	The formation redshift model of Correa et al. 2015.
`zfModelBera26`(M, z0[, f_form, C_0, C_1, ...])	The formation redshift model of Bera & Diemer 2026.
`AccretionHistoryModel`()	Characteristics of mass accretion history models.
`models_mah`	List of mass accretion history models.
`massAccretionHistory`(M0, z[, model, zf, ...])	Wrapper function for mass accretion history.
`mahModelWechsler02`(M0, z[, zf, z0])	The mass accretion history model of Wechsler et al. 2002.
`mahModelVanDenBosch02`(M0, z[, zf])	The mass accretion history model of van den Bosch 2002.
`mahModelMcBride09`(M0, z[, zf, beta, gamma, z0])	The mass accretion history model of McBride et al. 2009.
`mahModelVanDenBosch14`(M0, z[, statistic, ...])	The mass accretion history model of van den Bosch et al. 2014.
`mahModelCorrea15`(M0, z[, z0, ps_args, ...])	The mass accretion history model of Correa et al. 2015.
`AccretionRateModel`()	Characteristics of mass accretion rate models.
`models_mar`	List of mass accretion rate models.
`massAccretionRate`(M, z0[, model, delta_t, ...])	Wrapper function for mass accretion rate.
`marModelWechsler02Inst`(M, z0, zf)	The instantaneous accretion rate for the Wechsler et al. 2002 model.
`marModelCorrea15Inst`(M, z0[, ps_args, ...])	The instantaneous accretion rate for the Correa et al. 2015 model.
`marModelMcBride09Inst`(M, z0[, statistic])	The instantaneous accretion rate for the McBride et al. 2009 model.
`marModelRodriguezPuebla16`(M, z0[, inst])	The mass accretion rate model of Rodriguez-Puebla et al. 2016.
`marModelYung24`(M, z0)	The mass accretion rate model of Yung et al. 2024.
`marModelYung25`(M, z0)	The mass accretion rate model of Yung et al. 2025.
`marModelBera26`(M, z0[, tdyn_def, ps_args, ...])	The mass accretion rate model of Bera & Diemer 2026.

Module reference#

class halo.accretion_history.FormationRedshiftModel#

Characteristics of formation redshift models.

The models_zf dictionary contains one item of this class per available model.

Attributes:

func: callable: The model function.
ps_dependence: bool: Whether the model function accepts a ps_args kwarg.
sigma_dependence: bool: Whether the model function accepts a sigma_args kwarg.
deltac_dependence: bool: Whether the model function accepts a deltac_args kwarg.
requires_z0: bool: Whether the model function takes z0 as a positional argument.
default_mdef: str: The mass definition used by default when the wrapper function is called without an explicit mdef.

halo.accretion_history.models_zf = {'bera26': <halo.accretion_history.FormationRedshiftModel object>, 'correa15': <halo.accretion_history.FormationRedshiftModel object>, 'giocoli12': <halo.accretion_history.FormationRedshiftModel object>, 'lacey93': <halo.accretion_history.FormationRedshiftModel object>, 'mcbride09': <halo.accretion_history.FormationRedshiftModel object>, 'vandenbosch02': <halo.accretion_history.FormationRedshiftModel object>, 'vandenbosch14': <halo.accretion_history.FormationRedshiftModel object>}#: List of formation redshift models.

halo.accretion_history.formationRedshift(M, z0=0.0, model='bera26', mdef=None, c_model='diemer19', ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={}, **kwargs)#

Wrapper function for formation redshift.

Here we define the formation redshift as the half-mass redshift. Each formation redshift fit was calibrated against halos defined in a specific spherical overdensity mass definition (its default_mdef). If the user supplies M in a different definition, this function converts M to the model’s calibration mdef via changeMassDefinitionCModel() before evaluating the fit. If mdef is omitted, M is assumed to already be in the model’s calibration mdef and no conversion is performed.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
model: str: Identifier of the formation redshift model. See models_zf for possible values.
mdef: str: Mass definition in which M is given. If omitted, we assume that this is the same as the model’s calibration mdef (no conversion). If supplied and different from the model’s calibration mdef, M is converted internally using changeMassDefinitionCModel().
c_model: str: Concentration model used by the mass conversion (see mdef). Defaults to HALO_CONCENTRATION_MODEL.
ps_args: dict: Extra arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().
kwargs: dict: Extra arguments passed to the individual model function.

Returns:

zf: array_like: Formation redshift.

halo.accretion_history.zfModelLacey93(M, z0, f_form=0.5, C=np.float64(0.4769362762044699), zmax=30.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The formation redshift model of Lacey & Cole 1993.

The formation redshift is computed by inverting the EPS first-crossing condition

\[\delta_c(z_f) = \delta_c(z) + C \sqrt{2 [\sigma^2(f M) - \sigma^2(M)]}\]

for \(z_f\). With the default f_form = 0.5 and C = erfinv(0.5) this reproduces the Lacey & Cole 1993 median half-mass formation redshift. Other formation redshift models (van den Bosch 2002, Bera & Diemer 2026, etc.) are special cases obtained by varying f_form and C.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
f_form: float: Mass fraction defining formation (default 0.5).
C: float: The prefactor \(C\) in the barrier equation. Defaults to erfinv(0.5) ~ 0.4769.
zmax: float: Upper bracket for the root finder.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelVanDenBosch02(M, z0, f_form=0.254, C=0.477, zmax=30.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The formation redshift model of van den Bosch 2002.

This model (their Appendix A6) uses the original choice of f_form = 0.254 and C = 0.477. It is a special case of the Lacey & Cole 1993 construction and just forwards to zfModelLacey93().

The Eke-Cole-Frenk correction \(\delta_c \propto \Omega_{\rm m}(z)^{0.0055}\) is the original van den Bosch 2002 convention, so corrections=True is applied to deltac_args by default. The user can override by passing deltac_args={'corrections': False}.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
f_form: float: Mass fraction defining formation (default 0.254).
C: float: The prefactor \(C\) in the barrier equation (default 0.477).
zmax: float: Upper bracket for the root finder.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelMcBride09(M)#

The formation redshift model of McBride et al. 2009.

A fitting function for the half-mass formation redshift of halos observed at \(z = 0\), \(z_f = -0.24 \log_{10}(M / 10^{12} M_{\odot}) + 1.26\), where the mass is in physical solar masses (the input mass is divided by h internally).

Parameters:

M: array_like: Halo mass at \(z = 0\) in \(M_{\odot}/h\); can be a number or a numpy array.

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelGiocoli12(M, z0, f_form=0.5, zmax=30.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The formation redshift model of Giocoli et al. 2012.

The main-progenitor formation redshift of Giocoli, Tormen & Sheth 2012 (their Eqs. 9, 11, 13). The Giocoli+12 barrier \(\delta_c(z_f) = \delta_c(z) + \tilde\omega \sqrt{\sigma^2(fM) - \sigma^2(M)}\) with \(\tilde\omega = \sqrt{2 \ln(\alpha_f + 1)}\) is the Lacey & Cole 1993 barrier with \(C = \tilde\omega / \sqrt{2} = \sqrt{\ln(\alpha_f + 1)}\), so this routine just forwards to zfModelLacey93() with the appropriate prefactor.

The prefactor \(\alpha_f = 0.815 e^{-2 f^3} / f^{0.707}\) is a Giocoli+12 fit calibrated against EPS merger trees over the range \(0.1 \lesssim f_{\rm form} \lesssim 0.9\). The default f_form = 0.5 corresponds to the half-mass formation redshift; other choices give the formation redshift relative to the chosen mass fraction.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
f_form: float: Mass fraction defining formation (default 0.5). Any value in (0, 1) is accepted; the Giocoli+12 fit for \(\alpha_f\) was calibrated for roughly 0.1 to 0.9.
zmax: float: Upper bracket for the root finder.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelVanDenBosch14(M, z0, f_form=0.5, statistic='median', g_val=0.4, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The formation redshift model of van den Bosch et al. 2014.

The formation redshift is found by solving \(F(f) = \tilde\omega(f, z)\) for \(z\), where \(F\) is the universal mass accretion history function (their Eq. 20) and \(\tilde\omega\) the corrected time variable (their Eq. 19).

The exponent g_val sets the strength of the \(G^\gamma\) correction in Eq. (19). The default g_val = 0.4 is the value used by van den Bosch et al. 2014; g_val = 0.2 gives results similar to the Bera & Diemer 2026 model.

The Eke-Cole-Frenk correction \(\delta_c \propto \Omega_{\rm m}(z)^{0.0055}\) is the original van den Bosch et al. 2014 convention, so corrections=True is applied to deltac_args by default. The user can override by passing deltac_args={'corrections': False}.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
f_form: float: Mass fraction defining formation (default 0.5).
statistic: str: Which statistic the universal MAH function describes; can be median or mean.
g_val: float: The exponent of the \(G\) correction factor in Eq. (19).
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelCorrea15(M)#

The formation redshift model of Correa et al. 2015.

A fitting function for the formation redshift (Eq. 23 of their Paper I), \(z_f = -0.0064 (\log_{10} M)^2 + 0.0237 \log_{10} M + 1.8837\). The fit is for M in physical solar masses; the input mass is divided by h internally.

Parameters:

M: array_like: Halo mass at \(z = 0\) in \(M_{\odot}/h\); can be a number or a numpy array.

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

halo.accretion_history.zfModelBera26(M, z0, f_form=0.5, C_0=0.646, C_1=-0.062, zmax=30.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The formation redshift model of Bera & Diemer 2026.

Uses the Lacey & Cole 1993 inversion (see zfModelLacey93()) with the prefactor \(C(n_{\rm eff}) = C_0 + C_1\, n_{\rm eff}(M)\) and a default f_form = 0.5.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
f_form: float: Mass fraction defining formation (default 0.5).
C_0: float: Constant part of the C(n_eff) coefficient.
C_1: float: Slope of C with respect to the power-spectrum slope n_eff.
zmax: float: Upper bracket for the root finder.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

zf: array_like: Formation redshift; has the same dimensions as M.

class halo.accretion_history.AccretionHistoryModel#

Characteristics of mass accretion history models.

The models_mah dictionary contains one item of this class for each available model.

Attributes:

func: callable: The model function.
mdefs: list of str: Allowed spherical-overdensity mass definitions for the model. ['*'] means that the model accepts any definition.
default_mdef: str: Mass definition assumed by the model when no mdef is passed.
mdef_dependence: bool: Whether the model function depends on the mass definition directly (rather than only via the wrapper function’s mdef conversion).
ps_dependence: bool: Whether the model function accepts a ps_args kwarg.
sigma_dependence: bool: Whether the model function accepts a sigma_args kwarg.
deltac_dependence: bool: Whether the model function accepts a deltac_args kwarg.
uses_zf: bool: Whether the model takes a per-halo formation redshift through the zf kwarg.
default_zf_model: str: For models that take zf, the default formation redshift model to use when the user does not supply one explicitly.
accepts_z0: bool: Whether the model can be anchored at an arbitrary redshift z0 (instead of z = 0).

halo.accretion_history.models_mah = {'correa15': <halo.accretion_history.AccretionHistoryModel object>, 'mcbride09': <halo.accretion_history.AccretionHistoryModel object>, 'vandenbosch02': <halo.accretion_history.AccretionHistoryModel object>, 'vandenbosch14': <halo.accretion_history.AccretionHistoryModel object>, 'wechsler02': <halo.accretion_history.AccretionHistoryModel object>}#: List of mass accretion history models.

halo.accretion_history.massAccretionHistory(M0, z, model='wechsler02', zf=None, mdef=None, z0=0.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={}, **kwargs)#

Wrapper function for mass accretion history.

Parameters:

M0: array_like: Halo mass at the anchor redshift z0 in \(M_{\odot}/h\). Can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z). Must satisfy z >= z0 for the generalized anchor case.
model: str: Name of an entry in models_mah.
zf: None, str, float, or array_like: Formation redshift, which is accepted by models for which uses_zf = True. If None, the model’s default_zf_model is used to compute a formation redshift. If the name of a model in models_zf, that model is used instead. If a scalar, that formation redshift is used for all halo masses. If an array, that array must have the same dimensions as M0. This parameter is ignored by certain models that compute their formation redshift internally. If z0 > 0, the formation redshift model is evaluated at z0.
mdef: str: Mass definition in which M0 is given. If no mdef is given, the function falls back to the model’s default_mdef attribute. Some models can accept a range of mdefs.
any definition.
z0: float: Anchor redshift such that M0 is interpreted as M(z = z0). Only models with accepts_z0 = True support z0 > 0, while others raise an Exception.
ps_args: dict: Extra arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().
kwargs: dict: Extra arguments passed to the individual model function.

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\). The shape follows the input, namely a 2D array if both M0 and z are arrays, 1D if one of them is, and a scalar if both inputs are scalars.

halo.accretion_history.mahModelWechsler02(M0, z, zf=None, z0=0.0)#

The mass accretion history model of Wechsler et al. 2002.

An exponential mass accretion history, generalized to an arbitrary anchor redshift:

\[M(z) = M(z_0) \exp[-\alpha (z - z_0)], \qquad \alpha = \frac{\ln 2}{z_f - z_0}\]

where \(z_f\) is the formation redshift, evaluated by massAccretionHistory() at z = z_0. For z < z_0 the function returns NaN (the model only describes the main progenitor at z >= z_0).

Parameters:

M0: array_like: Halo mass at the anchor redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z).
zf: array_like: Formation redshift per halo (or a scalar broadcast over the halos), already computed by massAccretionHistory().
z0: float: Anchor redshift; M0 is interpreted as \(M(z = z_0)\).

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\); a 2D array with dimensions (number of halos, number of redshifts).

halo.accretion_history.mahModelVanDenBosch02(M0, z, zf=None)#

The mass accretion history model of van den Bosch 2002.

An exponential mass accretion history with the same functional form as mahModelWechsler02(), \(M(z) = M_0 \exp(-\alpha z)\), but with the corrected van den Bosch 2002 Appendix B2 conversion between the formation redshift and the rate constant, \(\alpha = (z_f \cdot 1.43)^{-1.05}\). By default zf comes from the van den Bosch 2002 formation redshift relation (f_form = 0.254, C = 0.477).

The model was calibrated for haloes observed at \(z = 0\), so M0 must be M(z=0). The wrapper function massAccretionHistory() always anchors at z=0 anyway, which matches this constraint.

Parameters:

M0: array_like: Halo mass at \(z = 0\) in \(M_{\odot}/h\); can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z).
zf: array_like: Formation redshift per halo (or a scalar broadcast over the halos), already computed by massAccretionHistory().

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\); a 2D array with dimensions (number of halos, number of redshifts).

halo.accretion_history.mahModelMcBride09(M0, z, zf=None, beta=0.1, gamma=None, z0=0.0)#

The mass accretion history model of McBride et al. 2009.

An analytic mass accretion history, generalized to an arbitrary anchor redshift:

\[M(z) = M(z_0) (1 + z - z_0)^\beta \exp[-\gamma (z - z_0)].\]

The shift \(z \to z - z_0\) makes the model self-similar about the anchor; the McBride+09 calibration is recovered at \(z_0 = 0\). The formation redshift \(z_f\) here is also measured relative to \(z_0\). If gamma is supplied, it is used directly; otherwise the rate constant is derived from the half-mass condition, \(\gamma = (\ln 2 + \beta \ln(1+z_f)) / z_f\).

Note

The McBride+09 \(z_f\) fit was calibrated for haloes observed at z = 0. When z0 > 0 and the supplied zf matches the McBride+09 zf formula (i.e. the user is relying on the model’s default), a warning is emitted suggesting an alternative formation-redshift model. Comparing against the McBride formula directly means the warning does not fire when the user has already opted into an alternative (e.g. zf='bera26').

Parameters:

M0: array_like: Halo mass at the anchor redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z).
zf: array_like: Formation redshift per halo (or a scalar broadcast over the halos), already computed by massAccretionHistory(). Ignored if gamma is supplied.
beta: float: The power-law exponent \(\beta\) (default 0.1).
gamma: float: The rate constant \(\gamma\). If None (default), it is derived from zf via the half-mass condition.
z0: float: Anchor redshift; M0 is interpreted as \(M(z = z_0)\).

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\); a 2D array with dimensions (number of halos, number of redshifts).

halo.accretion_history.mahModelVanDenBosch14(M0, z, statistic='median', g_val=0.4, z0=0.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The mass accretion history model of van den Bosch et al. 2014.

The universal mass accretion history model, where the MAH has a universal form when expressed as a function of a scaled time variable (their Eqs. 19 and 20). By default the MAH is anchored at z = 0. Pass z0 > 0 to anchor at the observation redshift, e.g. for use with high-z halos.

Parameters:

M0: array_like: Halo mass at the anchor redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z).
statistic: str: Which statistic the universal MAH function describes; can be median or mean.
g_val: float: The exponent of the \(G\) correction factor in Eq. (19).
z0: float: Anchor redshift; M0 is interpreted as \(M(z = z_0)\).
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\); a 2D array with dimensions (number of halos, number of redshifts).

halo.accretion_history.mahModelCorrea15(M0, z, z0=0.0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The mass accretion history model of Correa et al. 2015.

The generalized mass accretion history of their Paper III, Eq. 2:

\[\tilde M(z, M(z_0), z_0) = M(z_0)(1 + z - z_0)^{\tilde\alpha}\,e^{\tilde\beta (z - z_0)}\]

with \(\tilde\alpha\) and \(\tilde\beta\) given by Eqs. 3 and 4 of the paper, evaluated at the anchor redshift \(z_0\). The formation redshift fit (Eq. 8) is applied to M(z0) interpreted as the present-day-mass placeholder, and the variance \(\sigma\) is the standard EPS z = 0 quantity (per Eq. 6). For z < z_0 the function returns NaN.

Parameters:

M0: array_like: Halo mass at the anchor redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z: array_like: Redshift array at which to evaluate M(z).
z0: float: Anchor redshift; M0 is interpreted as \(M(z = z_0)\).
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

mz: array_like: The mass accretion history M(z) in \(M_{\odot}/h\); a 2D array with dimensions (number of halos, number of redshifts).

class halo.accretion_history.AccretionRateModel#

Characteristics of mass accretion rate models.

The models_mar dictionary contains one item of this class for each available model.

Attributes:

func: callable: The model function.
mdefs: list of str: Allowed spherical-overdensity mass definitions for the model. ['*'] means that the model accepts any definition.
default_mdef: str: Mass definition assumed by the model when no mdef is passed.
default_delta_t: str: Default time interval ('inst' or a dynamical time definition) used when the wrapper function is called without an explicit delta_t.
mdef_dependence: bool: Whether the model function depends on the mass definition directly (rather than only via the wrapper function’s mdef conversion).
ps_dependence: bool: Whether the model function accepts a ps_args kwarg.
sigma_dependence: bool: Whether the model function accepts a sigma_args kwarg.
deltac_dependence: bool: Whether the model function accepts a deltac_args kwarg.
supports_inst: bool: Whether the model provides a closed-form instantaneous (delta_t='inst') accretion rate. Note that models with has_mah_model=True support delta_t='inst' even if they have no closed form expression because the wrapper function computes it as a numerical derivative of the MAH.
supports_tdyn: bool: Whether the model provides a dynamical-time-averaged accretion rate (delta_t set to a dynamical-time definition such as 'crossing', 'peri', or 'orbit').
returns_mar: bool: Whether the model returns the dimensionless mass accretion rate \(\dot M / (M H(z_0))\) rather than \(\Gamma\). When True, the wrapper function converts it to \(\Gamma_{\rm inst}\) or \(\Gamma_{\rm dyn}\).
has_mah_model: bool: If True, the t_dyn-averaged rate is computed from the corresponding MAH model via the unified MAR-from-MAH algorithm (finite log-ratio between \(M(z_0)\) and \(M(z_i)\) where \(t_i = t_{\rm obs} - t_{\rm dyn}\)).
uses_zf: bool: Whether a MAH-based model function depends on a formation redshift zf. If True, the wrapper function calculates zf based on the user inputs (see massAccretionHistory()).

halo.accretion_history.models_mar = {'bera26': <halo.accretion_history.AccretionRateModel object>, 'correa15': <halo.accretion_history.AccretionRateModel object>, 'mcbride09': <halo.accretion_history.AccretionRateModel object>, 'rodriguezpuebla16': <halo.accretion_history.AccretionRateModel object>, 'vandenbosch14': <halo.accretion_history.AccretionRateModel object>, 'wechsler02': <halo.accretion_history.AccretionRateModel object>, 'yung24': <halo.accretion_history.AccretionRateModel object>, 'yung25': <halo.accretion_history.AccretionRateModel object>}#: List of mass accretion rate models.

halo.accretion_history.massAccretionRate(M, z0, model='bera26', delta_t=None, mdef=None, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={}, **kwargs)#

Wrapper function for mass accretion rate.

This function undertakes a number of conversions. If a model has returns_mar=True, the model function returns the dimensionless mass accretion rate \(\dot M / (M H(z_0))\), which this function converts to \(\Gamma_{\rm inst}\) or \(\Gamma_{\rm dyn}\) according to delta_t. If has_mah_model=True, we compute an averaged accretion rate numerically by evaluating the MAH at the respective times.

Some models provide an instantaneous rate (supports_inst=True), while for MAH-based models we evaluate it numerically. An instantaneous rate cannot be computed for models that return only an averaged accretion rate.

The output of the model of McBride et al. 2009 depends on the statistic keyword argument. If statistic=None, we translate the MAH model into a MAR. If statistic='mean'/'median', we use the direct fits for the accretion rate (see marModelMcBride09Inst()).

Parameters:

M: array_like

Halo mass at redshift z0 in \(M_{\odot}/h\). Can be a number or a numpy array.

z0: float

Observation redshift.

model: str

Name of an entry in models_mar.

delta_t: str

The time interval over which the rate is computed, can be:

None: the model’s default is used.
crossing (default): Averaged over one dynamical time at the observation redshift, defined as the crossing time from one side of the halo to another. When SO definitions are used, that crossing time is not a function of halo mass but only cosmology and redshift.
peri: Same as crossing, but the dynamical time is halved (i.e., the time to pericenter). Note that in other contexts (e.g., ROCKSTAR output), this is the definition called one dynamical time.
orbit: Same as crossing, but using the orbital time. See the dynamicalTime() function for more details on the time definitions.
inst: The instantaneous rate for \(\Delta_t \rightarrow 0\). For some models, this definition is available in closed-form, e.g., for a differentiable MAH fit. For any MAH-based model without a closed form it is computed numerically as a forward finite-difference derivative. We cannot reliably calculate the instantaneous rate for models that predict \(\Gamma\) over a wider time interval, since \(\Gamma\) evolves with redshift.

ps_args: dict

Extra arguments passed to matterPowerSpectrum().

sigma_args: dict

Extra arguments to sigma().

deltac_args: dict

Extra parameters passed to collapseOverdensity().

kwargs: dict

Extra arguments passed to the individual model function.

Returns:

gamma: array_like: Dimensionless accretion rate for each halo.

halo.accretion_history.marModelWechsler02Inst(M, z0, zf)#

The instantaneous accretion rate for the Wechsler et al. 2002 model.

The closed-form \(\Gamma_{\rm inst}\) for the Wechsler+02 exponential MAH \(M(z) = M(z_0) e^{-\alpha(z-z_0)}\) with \(\alpha = \ln 2 / (z_f - z_0)\):

\[\Gamma_{\rm inst} = \alpha (1 + z_0).\]

Note

This function provides only the instantaneous rate (delta_t='inst'). The dynamical-time averaged rate is handled by massAccretionRate() itself via the generic MAR-from-MAH algorithm.

Parameters:

M: array_like: Halo mass at redshift z0. This parameter is ignored because the exponential form of the MAH model means that its derivative is independent of the actual mass. The parameter appears in the function signature only to make it compatible with other models.
z0: float: Observation redshift.
zf: array_like: Formation redshift per halo, already computed by massAccretionRate().

Returns:

gamma: array_like: The dimensionless instantaneous accretion rate \(\Gamma_{\rm inst}\); has the same dimensions as M.

halo.accretion_history.marModelCorrea15Inst(M, z0, ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The instantaneous accretion rate for the Correa et al. 2015 model.

The closed-form \(\Gamma_{\rm inst}\) for the Correa+15 Paper III generalized MAH \(M(z) = M(z_0)(1+z-z_0)^{\tilde\alpha} e^{\tilde\beta (z-z_0)}\):

\[\Gamma_{\rm inst} = -(1 + z_0)(\tilde\alpha + \tilde\beta),\]

evaluated at the anchor where \(d\ln M/du = (\tilde\alpha + \tilde\beta)(1+z_0)\).

Note

This function provides only the instantaneous rate (delta_t='inst'). The dynamical-time averaged rate is handled by massAccretionRate() itself via the generic MAR-from-MAH algorithm.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

gamma: array_like: The dimensionless instantaneous accretion rate \(\Gamma_{\rm inst}\); has the same dimensions as M.

halo.accretion_history.marModelMcBride09Inst(M, z0, statistic=None, **mah_kwargs)#

The instantaneous accretion rate for the McBride et al. 2009 model.

The instantaneous \(\Gamma_{\rm inst}\) is computed from one of two sources, selected via statistic:

statistic=None (default) — the analytic McBride+09 MAH \(M(z) = M_0 (1+z-z_0)^\beta e^{-\gamma (z-z_0)}\) together with the half-mass formation redshift relation gives the closed form \(\Gamma_{\rm inst} = (1+z_0)(\gamma - \beta)\). The McBride+09 \(z_f\) fit is calibrated for present-day haloes only, so this branch is restricted to z0 = 0 unless the user supplies an alternative zf.
statistic='mean' or 'median' — the calibrated \(dM/dt\) fit from Eqs. 8 (mean) or 9 (median) of McBride+09 gives \(\Gamma_{\rm inst} = (dM/dt) / (M H(z_0))\). The fit is just a function of mass and redshift, so this branch works at any z0 >= 0.

Note

This function provides only the instantaneous rate (delta_t='inst'). The dynamical-time averaged rates are handled by massAccretionRate() itself: statistic=None via the generic MAR-from-MAH algorithm, and statistic='mean'/'median' via backward integration of the \(dM/dt\) fit over one dynamical time.

zf is only used when statistic=None; it has no effect for 'mean'/'median'.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
statistic: str: The source of the instantaneous rate; can be None (analytic MAH), 'mean', or 'median' (dM/dt fits).
mah_kwargs: dict: Extra arguments consumed by the statistic=None branch (zf, beta, gamma); see mahModelMcBride09() for their meaning.

Returns:

gamma: array_like: The dimensionless instantaneous accretion rate \(\Gamma_{\rm inst}\); has the same dimensions as M.

halo.accretion_history.marModelRodriguezPuebla16(M, z0, inst=False)#

The mass accretion rate model of Rodriguez-Puebla et al. 2016.

This function returns the dimensionless mass accretion rate \(\dot M / (M H(z_0))\) (massAccretionRate() converts it to \(\Gamma_{\rm inst}\) or \(\Gamma_{\rm dyn}\) according to delta_t).

inst=False (a dynamical-time-averaged delta_t) uses the single-power-law dynamical-time-averaged fit (Table 2 row 2 of the paper).
inst=True (delta_t='inst') uses the double-power-law instantaneous fit of Eq. 14 with parameters from Eqs. 15-17.

Both fits give \(\dot M\) in \(h^{-1} M_\odot / yr\); M is in \(M_{\odot}/h\), and the ratio \(\dot M / M\) is normalized by \(H(z_0)\) to give the dimensionless rate.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
inst: bool: If True, use the instantaneous fit; otherwise the dynamical-time-averaged fit.

Returns:

mar: array_like: The dimensionless mass accretion rate \(\dot M / (M H(z_0))\); has the same dimensions as M.

halo.accretion_history.marModelYung24(M, z0)#

The mass accretion rate model of Yung et al. 2024.

This function returns the dimensionless mass accretion rate \(\dot M / (M H(z_0))\) (their Eq. 5, normalized by \(M H(z_0)\)); massAccretionRate() converts it to \(\Gamma_{\rm dyn}\).

The fit was calibrated against the rate averaged over one dynamical time at \(z \geq 6\). An instantaneous evaluation is not provided by the authors, so this model only supports the t_dyn-averaged rate (massAccretionRate() does not allow delta_t='inst' for this model). A warning is emitted when z0 < 6 since the fit is extrapolating outside its calibrated range.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.

Returns:

mar: array_like: The dimensionless mass accretion rate \(\dot M / (M H(z_0))\); has the same dimensions as M.

halo.accretion_history.marModelYung25(M, z0)#

The mass accretion rate model of Yung et al. 2025.

This function returns the dimensionless mass accretion rate \(\dot M / (M H(z_0))\) (their Eqs. 4 and 5, normalized by \(M H(z_0)\)); massAccretionRate() converts it to \(\Gamma_{\rm dyn}\).

Same functional form as Yung+24 but revised coefficients from a fit to data over \(2 \leq z \leq 18\), with additional VSMDPL and gureft data at \(z > 15\). The fit was calibrated against the rate averaged over one halo virial dynamical time. An instantaneous evaluation is not provided by the authors, so this model only supports the t_dyn-averaged rate (massAccretionRate() does not allow delta_t='inst' for this model). A warning is emitted when z0 < 2 since the fit is not intended for use at lower redshift.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.

Returns:

mar: array_like: The dimensionless mass accretion rate \(\dot M / (M H(z_0))\); has the same dimensions as M.

halo.accretion_history.marModelBera26(M, z0, tdyn_def='crossing', ps_args={'model': 'eisenstein98', 'path': None}, sigma_args={}, deltac_args={})#

The mass accretion rate model of Bera & Diemer 2026.

The fit returns \(\Gamma_{\rm dyn}\) averaged over one dynamical time; massAccretionRate() does not allow delta_t='inst' for this model. The fit was calibrated for tdyn_def='crossing' — other values are accepted but a warning is emitted since the coefficients are used as-is without recalibration.

Parameters:

M: array_like: Halo mass at redshift z0 in \(M_{\odot}/h\); can be a number or a numpy array.
z0: float: Observation redshift.
tdyn_def: str: The definition of the dynamical time; see dynamicalTime(). The fit was calibrated for 'crossing'.
ps_args: dict: Arguments passed to matterPowerSpectrum().
sigma_args: dict: Extra arguments passed to sigma().
deltac_args: dict: Extra parameters passed to collapseOverdensity().

Returns:

gamma_dyn: array_like: The dimensionless accretion rate \(\Gamma_{\rm dyn}\) averaged over one dynamical time; has the same dimensions as M.