Colossus tutorial: mass accretion history, formation redshift, and accretion rate¶

Welcome to the Colossus halo mass accretion tutorial. The mass history of halos can be characterized in three closely related ways, all of which are implemented in the colossus.halo.mass_accretion module:

  • Formation redshift $z_{\rm f}$: The redshift at which the main progenitor of a halo had assembled a fixed fraction of its present-day mass. In this module, we use the half-mass redshift unless otherwise specified.
  • Mass accretion history $M(z)$: The mass of the main progenitor as a function of redshift. Given that some halos form earlier than others, any functional form for the accretion history relies on a formation redshift as input.
  • Mass accretion rate $\Gamma$: We express accretion rates in dimensionless units, which can either be $\Gamma_{\rm dyn} = \Delta \ln M / \Delta \ln a$, which is averaged over one dynamical time, or an instantaneous rate, $\Gamma_{\rm ins} = \dot M / (M H(z))$. These different quantities can be chosen by setting the time interval (delta_t='crossing' or delta_t='inst').

Each of the three quantities has its own registry of models and an umbrelly function that can handle all implemented models. It is recommended to use those functions rather than the individual model functions.

In [1]:
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

As always with Colossus, we begin by setting a cosmology.

In [2]:
from colossus.cosmology import cosmology
cosmology.setCosmology('planck18');

Import the mass-accretion module:

In [3]:
from colossus.halo import accretion_history as macc

Formation redshift¶

The available formation-redshift models are registered in the models_zf dictionary:

In [4]:
for model_name in macc.models_zf:
    print(model_name)

lacey93
vandenbosch02
mcbride09
giocoli12
vandenbosch14
correa15
bera26

The recommended entry point is formationRedshift(), which feeds through the power-spectrum, $\sigma(M)$, and collapse-overdensity argument dictionaries automatically. For example, the Correa et al. 2015 model gives for a halo of $10^{12}\,M_{\odot}/h$:

In [5]:
macc.formationRedshift(1E12, z0 = 0.0, model = 'correa15', mdef='vir')
Out[5]:
np.float64(1.2242758876617554)

The function also accepts an array of masses:

In [6]:
M = 10**np.arange(10.0, 15.0, 1.0)
macc.formationRedshift(M, z0 = 0.0, model = 'bera26')
Out[6]:
array([1.45592165, 1.22708689, 0.99682154, 0.77285832, 0.5610913 ])

Let's plot the formation redshift predictions of all models as a function of halo mass. Different models use different mass definitions (e.g., vir or 200m), but the function automatically converts the input halo masses.

If appropriate, the function can take additional arguments that are forwarded to the individual model functions. For example, the van den Bosch et al. 2014 model has a parameter called g_val, the exponent of the Parkinson et al. 2008 correction term. The default value is g_val = 0.4, but a modified value g_val = 0.2 brings the model closer to the Bera & Diemer 2026 calibration.

In [7]:
M_arr = 10**np.linspace(8.0, 15.0, 60)

plt.figure(figsize = (6.5, 5))
plt.xscale('log')
plt.xlabel(r'$M_{\rm vir}\ [M_{\odot}/h]$')
plt.ylabel(r'$z_{\mathrm{f}}$')

for model_name in macc.models_zf:
    zf = macc.formationRedshift(M_arr, z0 = 0.0, model = model_name, mdef = 'vir')
    plt.plot(M_arr, zf, lw = 2.0, label = model_name.replace('_', '\\_'))

# Explicit g_val = 0.2 variant of vdb14 — closely tracks bera26.
zf_vdb_mod = macc.formationRedshift(M_arr, z0 = 0.0, model = 'vandenbosch14', g_val = 0.2)
plt.plot(M_arr, zf_vdb_mod, '--', lw = 2.0, color = 'C3',
         label = r'vandenbosch14, $g = 0.2$')

plt.ylim(0.0, 3.0)
plt.legend();
No description has been provided for this image

Mass accretion history¶

The MAH sub-module follows the same pattern as formation redshift. A number of models for the MAH are available:

In [8]:
for model_name in macc.models_mah:
    print(model_name)

wechsler02
vandenbosch02
mcbride09
vandenbosch14
correa15

By default, MAH models are anchored at $z = 0$. Almost all models take a formation redshift as a parameter, but it can be given in a number of ways:

  • None: If no information is given, each model has a default model for the formation redshift that it falls back to.
  • A formation redshift model name: This model is used instead of the default.
  • A number: The MAHs for all masses given (if more than one) are evaluated with the same formation redshift.
  • An array: This array is assumed to contain a formation redshift for each halo mass for which the MAH is to be evaluated.

Below are some examples:

In [9]:
z = np.linspace(0.0, 6.0, 13)
macc.massAccretionHistory(1E12, z, model = 'wechsler02')  # default zf = 'bera26'
Out[9]:
array([1.00000000e+12, 7.06325800e+11, 4.98896136e+11, 3.52383213e+11,
       2.48897355e+11, 1.75802623e+11, 1.24173929e+11, 8.77072495e+10,
       6.19498932e+10, 4.37568079e+10, 3.09065624e+10, 2.18301024e+10,
       1.54191646e+10])

Swap in a different formation redshift model by name:

In [10]:
macc.massAccretionHistory(1E12, z, model = 'wechsler02', zf = 'mcbride09')
Out[10]:
array([1.00000000e+12, 7.52582586e+11, 5.66380549e+11, 4.26248139e+11,
       3.20786927e+11, 2.41418655e+11, 1.81687476e+11, 1.36734830e+11,
       1.02904252e+11, 7.74439483e+10, 5.82829669e+10, 4.38627460e+10,
       3.30103388e+10])

Same call with a scalar zf override:

In [11]:
macc.massAccretionHistory(1E12, z, model = 'wechsler02', zf = 2.0)
Out[11]:
array([1.00000000e+12, 8.40896415e+11, 7.07106781e+11, 5.94603558e+11,
       5.00000000e+11, 4.20448208e+11, 3.53553391e+11, 2.97301779e+11,
       2.50000000e+11, 2.10224104e+11, 1.76776695e+11, 1.48650889e+11,
       1.25000000e+11])

Models that do not use a formation redshift internally (e.g., van den Bosch et al. 2014) simply ignore the zf argument. Let's compare a few MAH models for the same halo mass at $z = 0$. Note that MAH models do not care about the SO mass definition. Although different models were calibrated for different definitions, they really quantify the relative mass compared to some final mass. Thus, we assume that they work for any input mass definition.

In [12]:
z = np.linspace(0.0, 8.0, 50)
M0 = 1E12

plt.figure()
plt.yscale('log')
plt.xlabel(r'$z$')
plt.ylabel(r'$M(z)\ [M_{\odot}/h]$')

for model in macc.models_mah:
    mz = macc.massAccretionHistory(M0, z, model = model)
    plt.plot(z, mz, label = model.replace('_', '\\_'))

plt.xlim(0.0, 8.0)
plt.legend();
No description has been provided for this image

Some MAH models do not need to be anchored at $z = 0$. Those models accept an optional z0 keyword, namely wechsler02, vandenbosch14, and correa15. For example, let's plot some models anchored at $z_0 = 2$:

In [13]:
z = np.array([3.0, 4.0, 6.0, 8.0])
for m in ['wechsler02', 'vandenbosch14', 'correa15']:
    mz = macc.massAccretionHistory(1E12, z, model = m, z0 = 2.0)
    print(f'{m:14s}  M(z) = {mz}')

wechsler02      M(z) = [4.35690432e+11 1.89826152e+11 3.60339681e+10 6.84018953e+09]
vandenbosch14   M(z) = [4.85674116e+11 2.20830589e+11 5.30062558e+10 1.51943677e+10]
correa15        M(z) = [4.69914833e+11 2.22120804e+11 4.99333253e+10 1.12651703e+10]

Models that don't support z0 raise an exception if a non-zero value is passed.

In [14]:
mz = macc.massAccretionHistory(1E12, z, model = 'mcbride09', z0 = 2.0)

/var/folders/9c/lqfp3yps2951jvnbfxhl8svh0000gr/T/ipykernel_27746/916919872.py:1: UserWarning: mcbride09 MAH at z0 = 2 with the default McBride+09 zf fit: the zf fit was calibrated for z = 0 haloes only. Consider passing an alternative formation-redshift model via zf (e.g. zf='bera26') for z0 > 0.
  mz = macc.massAccretionHistory(1E12, z, model = 'mcbride09', z0 = 2.0)

Mass accretion rate¶

There are a number of different definitions of the MAR in the literature, and the main MAR function can convert them into one another (within reason). Fundamentally, we express all rates in dimensionless units of $\Gamma \equiv \Delta \ln M / \Delta \ln a$ (as opposed to rates in $M_\odot/{\rm yr}$, for example). However, the time interval over which the rate is measured does matter. The delta_t keyword allows for a number of choices:

  • crossing (default): Averaged over one dynamical time at the observation redshift, defined as the crossing time from one side of the halo to another (the definition of Diemer 2017). 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 mass_so.dynamicalTime() function for more details on the time definitions.
  • inst: The instantaneous rate for $\Delta t \rightarrow 0$. In this limit, the expression for $\Gamma$ becomes $\Gamma_{\rm inst} = \dot M / (M H(z))$. For some models, this definition is available in closed-form, e.g., for a differentiable MAH fit (Wechsler et al. 2002, McBride et al. 2009, Rodriguez-Puebla et al. 2016, Correa et al. 2015). For any MAH-based model without a closed form (e.g., van den Bosch et al. 2014) 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.

Another important factor in the dynamical time is the mass definition (vir, 200m and such). This definition matters only marginally for the $\Delta \ln M$ term because the absolute halo mass is scaled out, but the dynamical time does increase for larger radius definitions. Thus, the user can also pass an mdef parameter that is used to calculate the correct definition.

While the main MAR function attempts to convert any model to the desired definition of $\Gamma$, the individual model functions return only the definition they were calibrated for. McBride et al. 2009 suggested a number of models. Their MAH model translates to a MAR as discussed above, although this is calibrated only for $z = 0$ (for the default formation redshift model). They also present fitting functions for the the instantaneous mean and median MAR, which can be evaluated using statistic='mean'|'median'.

The full list of MAR models is:

In [15]:
for name, m in macc.models_mar.items():
    print(f'{name:20s}  closed-form inst = {m.supports_inst}')

wechsler02            closed-form inst = True
mcbride09             closed-form inst = True
vandenbosch14         closed-form inst = False
correa15              closed-form inst = True
rodriguezpuebla16     closed-form inst = True
yung24                closed-form inst = False
yung25                closed-form inst = False
bera26                closed-form inst = False

Let's evaluate the default model at $M = 10^{12}\,M_{\odot}/h$, $z = 2$:

In [16]:
macc.massAccretionRate(1E12, 2.0, model = 'bera26')
Out[16]:
2.7692072638395704

We can vectorize the halo mass input:

In [17]:
M = 10**np.linspace(8.0, 15.0, 8)
macc.massAccretionRate(M, 2.0, model = 'bera26')
Out[17]:
array([1.12176112, 1.35193105, 1.6709104 , 2.11802426, 2.76920726,
       3.72899418, 5.17632112, 7.3992663 ])

Asking for the instantaneous rate of a model that does not support it raises an error:

In [18]:
try:
    macc.massAccretionRate(1E12, 2.0, model = 'bera26', delta_t = 'inst')
except Exception as e:
    print('Got expected error:', e)

Got expected error: Model bera26 does not provide an instantaneous (delta_t='inst') accretion rate; it returns Gamma averaged over one dynamical time only.

MAR as a function of redshift at fixed mass¶

Let's make a plot of $\Gamma_{\rm dyn}$ for a Milky-Way-scale halo across the various models:

In [19]:
z = np.linspace(0.0, 14.0, 30)
M0 = 1E12

plt.figure()
plt.xlabel(r'$z$')
plt.ylabel(r'$\Gamma_{\rm dyn}$')

mar_kwargs = {'mcbride09': {'statistic': 'median'}}
for model in ['wechsler02', 'mcbride09', 'rodriguezpuebla16', 'correa15', 'vandenbosch14', 'bera26']: 
    g = np.array([macc.massAccretionRate(M0, zi, model = model, **mar_kwargs.get(model, {})) for zi in z])
    plt.plot(z, g, label = model.replace('_', '\\_'))

plt.xlim(0.0, 14.0)
plt.ylim(0.0, 10.0)
plt.legend();
No description has been provided for this image

MAR as a function of mass at fixed redshift¶

Let's see how the MAR evolves with redshift, using the default model (bera26):

In [20]:
M = 10**np.linspace(8.0, 15.0, 60)

plt.figure()
plt.xscale('log')
plt.xlabel(r'$M\ [M_{\odot}/h]$')
plt.ylabel(r'$\Gamma_{\rm dyn}$')
for z in [0.0, 1.0, 2.0, 4.0, 6.0, 8.0, 10, 12]:
    g = macc.massAccretionRate(M, z, model = 'bera26', mdef='200m', delta_t='crossing')
    plt.plot(M, g, label = r'$z = %.1f$' % z)
plt.xlim(1E8, 1E15)
plt.ylim(0.0, 10.0)
plt.legend();
No description has been provided for this image

Instantaneous vs. averaged rates¶

Below we compare the models that support both instantaneous and dynamical-time-averaged MARs:

In [21]:
z = np.linspace(0.0, 8.0, 30)
M0 = 1E12

fig, axes = plt.subplots(1, 2, figsize = (12, 4.5), sharey = True)
for ax, dt in zip(axes, ['crossing', 'inst']):
    mar_kwargs = {'mcbride09': {'statistic': 'median'}}
    for model in ['wechsler02', 'mcbride09', 'rodriguezpuebla16', 'correa15', 'vandenbosch14']:
        g = np.array([macc.massAccretionRate(M0, zi, model = model, delta_t = dt,
                                              **mar_kwargs.get(model, {})) for zi in z])
        ax.plot(z, g, label = model.replace('_', '\\_'))
    ax.set_xlabel(r'$z$')
    ax.set_title(r"delta_t = '%s'" % dt)
    ax.set_xlim(0, 8); ax.set_ylim(0, 8)
axes[0].set_ylabel(r'$\Gamma$')
axes[0].legend();
No description has been provided for this image

Consistency between the MAH and the MAR¶

Because the dynamical-time-averaged $\Gamma_{\rm dyn}$ is simply a finite difference of $\ln M$, the MAR and MAH are self-consistent. For example, finite-differencing the z₀-anchored vandenbosch14 MAH should reproduce the direct vandenbosch14 MAR call.

In [22]:
from colossus.halo import mass_so

cosmo = cosmology.getCurrent()
z_obs = 2.0
M_obs = 1E12

t_obs = cosmo.age(z_obs)
t_dyn = mass_so.dynamicalTime(z_obs, '200m', definition = 'crossing')
z_i = cosmo.age(t_obs - t_dyn, inverse = True)

mz = macc.massAccretionHistory(M_obs, np.array([z_i]), model = 'vandenbosch14', z0 = z_obs)
psi_i = float(mz[0]) / M_obs
gamma_from_mah = np.log(psi_i) / np.log((1 + z_obs) / (1 + z_i))
gamma_direct = macc.massAccretionRate(M_obs, z_obs, model = 'vandenbosch14')
print('Gamma from anchored MAH: %.6f' % gamma_from_mah)
print('Gamma from MAR wrapper:  %.6f' % gamma_direct)

Gamma from anchored MAH: 2.405216
Gamma from MAR wrapper:  2.405216

For the broader theoretical and empirical context, see the references in the mass_accretion module documentation and references therein.