1. Introduction
Within the hierarchical structure formation framework, dark matter halos serve as the building blocks in which the baryonic diffuse matter collapses to form galactic discs (White & Rees Reference White and Rees1978; Fall & Efstathiou Reference Fall and Efstathiou1980; Mo, Mao, & White Reference Mo, Mao and White1998; Zavala, Okamoto, & Frenk Reference Zavala, Okamoto and Frenk2008). There is an extensive literature on how dark matter halos gain their angular momentum through gravitational tidal interactions with the large scale mass distribution of the universe (Peebles Reference Peebles1969; Doroshkevich Reference Doroshkevich1970; Efstathiou & Jones Reference Efstathiou and Jones1980; White Reference White1984). This theory, known as Tidal Torque Theory, has been shown to be consistent with the results of cosmological N-body simulations and captures well the growth of the angular momentum of proto-halos (Barnes & Efstathiou Reference Barnes and Efstathiou1987; Sugerman, Summers, & Kamionkowski Reference Sugerman, Summers and Kamionkowski2000; Porciani, Dekel, & Hoffman Reference Porciani, Dekel and Hoffman2002a,Reference Porciani, Dekel and Hoffmanb). However, its validity breaks down in the non-linear regime, when halos have collapsed and virialised (White Reference White1984; Barnes & Efstathiou Reference Barnes and Efstathiou1987; Porciani et al. Reference Porciani, Dekel and Hoffman2002b; Bailin & Steinmetz Reference Bailin and Steinmetz2005), and this is understood to reflect the importance of merger and accretion in the assembly histories of halos (Gardner Reference Gardner2001; Maller, Dekel, & Somerville Reference Maller, Dekel and Somerville2002; Vitvitska et al. Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002; Benson, Behrens, & Lu Reference Benson, Behrens and Lu2020), and consequently the galaxies that they host.
In this context, cosmological simulations provide crucial insight into how halos and galaxies acquire their angular momentum over time. The properties of galaxies seem to be strongly correlated to those of the dark matter halo in which they are embedded (Mo et al. Reference Mo, Mao and White1998; Dutton et al. Reference Dutton, van den Bosch, Dekel and Courteau2007). For instance, a halo’s angular momentum can be connected to properties of galactic discs, such as their size and surface density (Fall & Efstathiou Reference Fall and Efstathiou1980), which in turn affects their star formation rates (van den Bosch Reference van den Bosch2000). This relationship between halo angular momentum and the intrinsic properties of galaxies is one of the core assumptions underpinning semi-analytic models (SAMs) of galaxy formation, which evolve baryons in the background of dark matter halos by means of a coupled set of differential equations which parameterise key astrophysical processes, such as cooling, star formation, and stellar feedback (Kauffmann et al. Reference Kauffmann, Colberg, Diaferio and White1999a,Reference Kauffmann, Colberg, Diaferio and Whiteb; van den Bosch Reference van den Bosch2000; Croton et al. Reference Croton2006). This approach has successfully reproduced many of the global properties of the galaxy population across cosmic time (e.g. Knebe et al. Reference Knebe2015).
Despite extensive work on the spin distribution at fixed epoch (e.g. Vitvitska et al. Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002; Hetznecker & Burkert Reference Hetznecker and Burkert2006; Bett et al. Reference Bett, Eke, Frenk, Jenkins, Helly and Navarro2007; Macciò et al. Reference Macciò, Dutton and van den Bosch2008; Macciò, Dutton, & van den Bosch Reference Macciò, Dutton, van den Bosch, Moore, Potter and Stadel2007; Knebe & Power Reference Knebe and Power2008; Muñoz-Cuartas et al. Reference Muñoz-Cuartas, Macciò, Gottlöber and Dutton2011; Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014; Zjupa & Springel Reference Zjupa and Springel2017), only a limited subset of studies has quantified its redshift evolution. In particular, Hetznecker & Burkert (Reference Hetznecker and Burkert2006) present one of the few explicit evolutionary analyses – but, as we find in this paper, their results are affected by limited mass resolution.Footnote
a
Similarly, Ahn et al. (Reference Ahn, Kim, Shin, Kim and Choi2014) track
$\lambda(z)$
and
$\lambda^\prime(z)$
directly and find broadly consistent qualitative trends, but focus primarily on mass dependence and provide no analytical fitting relations for either the mean or the dispersion of the log-normal distribution as a function of redshift, making their results difficult to implement directly in models that require sampling from a log-normal at arbitrary z. Zjupa & Springel (Reference Zjupa and Springel2017) confirm the log-normal form of the spin distribution at multiple redshifts in the Illustris simulation (their Figures 9, 11, 18), but do not extract or model the redshift evolution of the distribution parameters quantitatively. Consequently, to our knowledge, no existing work provides closed-form, resolution-validated expressions for both
$\mu(\ln\lambda, z)$
and
$\sigma(z)$
simultaneously – the combination required to fully specify the spin distribution at any redshift for use in forward modelling. This is the gap our paper addresses.
The most common means to characterise a dark matter halo’s angular momentum is the dimensionless spin parameter
$\lambda$
. Peebles (Reference Peebles1969) defined the spin parameter as,
where G is Newton’s gravitational constant;
$J=\left|\sum_{i=0}^{N_p}\vec{r_i}\times \vec{p_i}\right|$
is the magnitude of the angular momentum of
$N_p$
particles in the halo, each of mass
$M/N_p$
where M is the mass of the halo; and
$|E| = |E_{\text{kin}} + E_{\text{pot}}|$
is the absolute value of the halo’s total energy where
$E_{\text{kin}}$
is its kinetic energy and
$E_{\text{pot}}$
is its gravitational potential energy. We can interpret Equation (1) as the halo’s rotational velocity (
$v_{\text{rot}}\sim\,J/MR$
) normalised by its velocity dispersion (
$\sigma\sim\,\sqrt{GM/R}\sim\sqrt{|E|/M}$
).
This classic Peebles (Reference Peebles1969) definition has been supplemented subsequently by the Bullock et al. (Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001) definition,
where
$V=\sqrt{GM/R}$
is the circular velocity of the halo. We can interpret Equation (2) as the halo’s angular momentum (J) normalised by the angular momentum we would expect if halo material was on circular orbits (
$\sim \sqrt{GMR}$
).
$\lambda'$
has the advantage that it is computationally cheaper to obtain than
$\lambda$
because it avoids the expensive calculation of the gravitational potential energy,
$E_{\text{pot}}$
(Vitvitska et al. Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002; Knebe & Power Reference Knebe and Power2008). Moreover, its value is less sensitive to the definition of a halo’s boundary (Bullock et al. Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001), which also makes it a versatile measure for the spin parameter of galaxies (Burkert & D’Onghia Reference Burkert, D’Onghia, Block, Puerari, Freeman, Groess and Block2004; Tonini et al. Reference Tonini, Lapi, Shankar and Salucci2006; Obreja, Buck, & Macciò Reference Obreja, Buck and Macciò2022).
There is a substantial body of work on the growth of the spin parameter of halos driven by the large-scale gravitational tidal field in the linear regime (e.g. Peebles Reference Peebles1969; Doroshkevich Reference Doroshkevich1970; White Reference White1984; Catelan & Theuns Reference Catelan and Theuns1996). As halos enter the non-linear regime, their self-gravity dominates the large-scale tidal torques (e.g. Sugerman et al. Reference Sugerman, Summers and Kamionkowski2000) and merger and accretion events are predicted to become important (e.g. Maller et al. Reference Maller, Dekel and Somerville2002; Vitvitska et al. Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002; Peirani, Mohayaee, & de Freitas Pacheco Reference Peirani, Mohayaee and de Freitas Pacheco2004; Hetznecker & Burkert Reference Hetznecker and Burkert2006; Sharma, Steinmetz, & Bland-Hawthorn Reference Sharma, Steinmetz and Bland-Hawthorn2012) – although the precise role of mergers in setting the spin distribution remains debated (e.g. D’Onghia & Navarro Reference D’Onghia and Navarro2007; Wang & White Reference Wang and White2009).
Maller et al. (Reference Maller, Dekel and Somerville2002) used a SAM to demonstrate that both tidal torques and orbital angular momentum transfer during mergers can produce similar halo spin distributions. Vitvitska et al. (Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002) used N-body simulations to model spin growth as a random walk driven by uncorrelated orbital angular momenta of merging subhalos. This framework has been refined subsequently by Benson et al. (Reference Benson, Behrens and Lu2020), who showed that accurately reproducing the N-body spin distribution requires including a weak correlation between the spins of progenitor halos and the angular momenta of their infalling subhalos, and by Hou, Luo, & Lacey (Reference Hou, Luo and Lacey2025), who used merger trees drawn from the Millennium and Millennium-II simulations to construct a stochastic model in which spin vector evolution is driven by the accretion of mass and angular momentum. Peirani et al. (Reference Peirani, Mohayaee and de Freitas Pacheco2004) found that halo angular momentum evolves in two distinct phases, with mergers increasing spin and smooth accretion reducing it, but despite these variations, the overall spin parameter distribution remains largely unchanged with redshift, apart from a slight increase in the mean value. Sharma et al. (Reference Sharma, Steinmetz and Bland-Hawthorn2012) showed that mergers redistribute angular momentum, causing the inner halo to lose spin while the outer regions gain it, leading to a net spin reduction of 40–80% per Gyr.
Interestingly, Hetznecker & Burkert (Reference Hetznecker and Burkert2006) explicitly tracked the growth of both the Peebles (
$\lambda$
) and Bullock (
$\lambda'$
) spin parameters with redshift and reported contrasting behaviour in their respective evolutions. They found that
$\lambda'$
remains nearly constant but
$\lambda$
increases with redshift. Evolution is driven by minor mergers and accretion rather than major mergers. Although major mergers boost
$\lambda$
– by a factor of
$\sim$
$1.3$
– their infrequency makes them a secondary influence. Similar trends were reported in Macciò et al. (Reference Macciò, Dutton and van den Bosch2008) and Ahn et al. (Reference Ahn, Kim, Shin, Kim and Choi2014).
It is worth noting that D’Onghia & Navarro (Reference D’Onghia and Navarro2007) argued that the apparent spin-up associated with major mergers is largely an artefact of unrelaxed systems, and that the virialisation process itself leads to a net decrease in spin through the ejection of high-angular-momentum material beyond the virial radius. More fundamentally, Wang & White (Reference Wang and White2009) demonstrated that universal halo properties, including the spin distribution, are reproduced even in a hot dark matter universe where structure forms through monolithic collapse rather than hierarchical merging, suggesting that mergers may not be the necessary ingredient for establishing the log-normal form of the spin distribution.
Motivated by the contrasting findings in the literature, we revisit in this paper the evolution of the Peebles
$\lambda$
and Bullock
$\lambda^\prime$
spin parameters across redshift, with two specific aims: to provide a statistically robust, resolution-validated characterisation of how the mean and dispersion of the log-normal spin distribution evolve from
$z=5$
to
$z=0$
; and to deliver closed-form fitting functions that allow physically motivated spin values to be drawn at any redshift within this range without recourse to simulations. Although previous studies have noted that
$\lambda$
and
$\lambda^\prime$
evolve differently (e.g. Hetznecker & Burkert Reference Hetznecker and Burkert2006; Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014), no existing work provides analytical expressions for both the mean and the width of the distribution simultaneously – the combination required to fully specify the spin distribution at arbitrary redshift. We provide these here and discuss the physical origin of the differences between the two definitions in Section 4.
The main focus of this work is to quantify, for a dark matter halo sample found in cosmological simulations, the redshift evolution of both spin definitions (Equations 1 and 2). The approach will be to study the statistical properties of the spin in the halo population. To this end, this work is structured as follows: in Section 2, we introduce our suite of simulations, the halo finder used to identify dark matter halos and provide a brief discussion on the spin and the halo sample’s basic properties. Secondly, in Section 3 we showcase our main results, based upon studying in detail the redshift evolution of the spin probability distribution, eventually providing fitting functions that describe this evolution. Lastly, in Section 4 we present a discussion of our results and state the main takeaways of our experiments.
Simulation parameters in the B-suite (columns 1–3), number of halos
$N_h(z)$
found by AHF with more than 500 particles at redshifts
$z = 0.0$
,
$1.1$
and
$5.0$
(columns 4–6, respectively), and number of field halos
$N_h^{f}(z)$
after additionally filtering out subhalos (columns 7–9).

2. Data and methodology
A detailed understanding and characterisation of the spin parameter of dark matter halos requires numerical simulations that accurately model the relevant physical processes – in this case, the gravitational collapse and dynamical evolution of halos in a cosmological background – with sufficient resolution to produce reliable results and recover statistical samples of halos throughout a wide range of masses and scales; a halo finding code that correctly identifies and characterises the overdensities as bound objects; and appropriate post-processing pipelines for their analysis. In this section we introduce all these technical details and provide some basic discussion of the resulting dataset, setting the ground for the main results of our work.
2.1. Simulations and halo finder
Simulation Details: We have run four dissipationless
$\Lambda$
-Cold Dark Matter (
$\Lambda \mathrm{CDM}$
) cosmological simulations using GADGET4 (Springel et al. Reference Springel, Pakmor, Zier and Reinecke2021). Initial conditions at redshift
$z=63$
were created on-the-fly using the internal N-GENIC functionality in GADGET4 and configured to generate and use a realisation of the Eisenstein & Hu (Reference Eisenstein and Hu1999) power spectrum of primordial density fluctuations. We stored 48 output snapshots between
$z=5$
and
$z=0$
, uniformly distributed in
$\log_{10}(a)$
space. The mass resolution of our runs is insufficient to explore higher redshifts reliably (i.e.
$z\geq\,5$
), but we consider the redshift range large enough for the conclusions drawn from our data.
Each simulation follows the evolution of
$N=256^3$
equal mass particles and vary only in the box size B, whose values are given in column 1 of Table 1. The particle masses varies as
$m_{\mathrm{sim}} = \rho_{\mathrm{crit}} \Omega_0 \left(B/N\right)^3$
, where
$\rho_{\mathrm{crit}} = 3H_0^2/8\pi G$
, as shown in column 2 of Table 1. Even with this modest value of N we can probe a wide range of halo masses from
$z=5$
to the present day. We choose a gravitational softening of
$\epsilon = 0.025\,B/N^{1/3}$
, which is intermediate between the converged and optimal softenings of Ludlow, Schaye, & Bower (Reference Ludlow, Schaye and Bower2019); the values used in each of the simulations is shown in Table 1. We adopt cosmological parameters of
$\Omega_0 = 0.315$
,
$\Omega_{\Lambda}=0.685$
,
$h=67.4$
and
$\sigma_8 = 0.810$
and
$n_s=0.965$
, in agreement with the latest observations (Planck Collaboration et al. Reference Collaboration2020).
It is worth noting that, although there is an abundance of state-of-the-art simulation datasets publicly available, there are good reasons for running our own suite of simulations tailored to our specific needs. First, we require a uniform minimum particle number across the different box sizes to maintain a common resolution threshold over a broad halo-mass range and redshift interval. Second, we need both spin definitions to be computed in a consistent manner – crucially, the Peebles spin needs gravitational potential energy, which is not generally available in public catalogues, and post-processed potential energy calculations across dozens of snapshots and simulation suites are computationally prohibitive. Third, our science goal demands explicit resolution-convergence tests and volume-dependence checks, which we provide in the Appendix; this shows that resolution and box size influence the inferred
$\lambda^\prime$
evolution, underscoring the need for controlled systematics.
Halo Finding: Dark matter halos were identified using the Amiga Halo Finder algorithmFootnote
b
(AHF; see Knollmann & Knebe Reference Knollmann and Knebe2009). The halo mass is determined by truncating the halo edge at a distance R from the centre, where the mean density within the halo volume goes below a predefined threshold of
$\Delta(z)$
multiplied by a reference density
$\rho_{\mathrm{ref}}(z)$
. For our study, we will focus on the
$R_{200}$
definition, which uses
$\Delta(z)=200$
and
$\rho_{\mathrm{ref}}(z) = \rho_{\mathrm{crit}}(z)$
, determining the halo masses as (Knebe et al. 2008; Knebe et al. Reference Knebe2013)
Note that we show in Appendix A our results for when
$R_{\mathrm{vir}}$
is defined by
$\Delta(z)$
as provided by the spherical top-hat collapse model and
$\rho_{\mathrm{ref}}(z) = \Omega_m(z) \rho_{\mathrm{crit}}(z)$
is the reference density (cf. Knebe et al. Reference Knebe2013). We find that while the results show only a minor quantitative dependence on the definition of the halo’s boundary, the qualitative behaviour and key conclusions remain consistent with those presented in the main analysis.
2.2 Halo sample
This section provides a succinct summary of the properties of our halo sample. We investigate three key aspects: first, we ensure that our sample spans a broad mass range with sufficient resolution (Section 2.2.1); second, we assess the correlation between spin and mass and compare to previous studies (Section 2.2.2); and third, we characterise the spin probability distribution and compare with theoretical expectations (Section 2.2.3). This last point is particularly important because the statistical properties of the spin distribution play a crucial role in our analysis. By addressing these aspects, we establish a solid framework for interpreting the results presented in Section 3.
2.2.1. Cumulative halo mass function
One of the main tools to study the consistency of a cosmological simulation and its mass range is the cumulative halo mass function
$N(\!\gt M_{\mathrm{200}})$
, hereafter referred to as cHMF.
Figure 1 shows the volume-normalisedFootnote
c
cHMF of the dark matter halos found in each of our simulations: blue, orange, green and red correspond, respectively, to B20, B40, B80, and B160 (as designated in Table 1), at redshifts
$z=0.0$
(solid),
$1.1$
(dashed curves), and
$5.0$
(dotted). The theoretical results, calculated using the HMFcalc code (Murray, Power, & Robotham Reference Murray, Power and Robotham2013) at the corresponding redshifts, are shown as solid black curves. Furthermore, the mass of halos with 500 particles in each simulation is marked by vertical lines, coloured according to their corresponding simulation. This minimum particle number will be justified in more detail below, and it sets the mass limit of objects considered in our study.
Volume-normalised cumulative halo mass functions of our
$N=256^3$
suite of simulations with
$B=20$
(blue), 40 (orange), 80 (green), and
$120 \, \mathrm{Mpc} \, h^{-1}$
(red) at redshifts
${z = 0.0}$
(solid),
${z = 1.1}$
(dashed) and
${z = 5.0}$
(dotted). The solid black curves indicate the analytical cumulative mass distributions obtained using the HMFcalc code (Murray et al. Reference Murray, Power and Robotham2013) at each redshift, whereas the vertical lines indicate the mass of halos with 500 particles in each simulation, coloured accordingly. The leftmost end of the curves correspond to the minimum 20-particle limit for a halo.

The main point of Figure 1 is to show that – stacked together – our suite of simulations covers a range of halo masses from approximately
$10^{10}h^{-1}\,\mathrm{M}_{\odot}$
to close to
$10^{15}h^{-1}\,\mathrm{M}_{\odot}$
– smaller simulation volumes (e.g. B20) provide insight into the low-mass halo regime, whereas larger volume (e.g. B160) probe the highest mass halos. The combination of the B-suite results in a continuous halo mass function which covers a wide mass range and follows well the HMFcalc theoretical expectations. Note, however, that we apply two additional selection criteria.
Halos must contain at least 500 particles: This ensures sufficient resolution of a halo’s internal properties to compute its spin parameter (see, for instance, Appendix A in Allgood et al. Reference Allgood, Flores, Primack, Kravtsov, Wechsler, Faltenbacher and Bullock2006, for a study of the resolution effects on different halo measurements).Footnote d
Substructure halos are removed: Tidal stripping of substructure and its gravitational interactions with the host halo affect its angular momentum and spin. In addition, unambiguous definition of material associated with the halo is difficult for substructure.
For the remainder of this work, we combine halos from all simulations that satisfy these criteria into a single catalogue, which we summarise in Table 1 (cf.
$N_h^{f}$
, columns 7–9).
2.2.2. Spin-mass correlation
If we are to reliably identify redshift-dependent trends in the spin parameter, then we must account for a potential degeneracy with halo mass. This is because the maximum mass decreases with increasing redshift (see Figure 1), and so an observed evolution of spin with redshift may reflect the underlying trend with halo mass rather than a genuine evolutionary effect. To address this, in Figure 2 we investigate the correlation between the Peebles spin (
$\lambda$
), left panels), and Bullock spin (
$\lambda'$
, right panels) and halo mass at three representative redshifts –
$z = 0$
,
$z = 1.1$
, and
$z = 5$
(top to bottom) – corresponding to those in Figure 1. Each panel also includes the Spearman rank correlation coefficient for the halo sample at that redshift.
Correlation between spin and mass at different redshifts z, coloured by density of data points. Each row indicates, from top to bottom respectively, redshifts
${z} = 0.0$
,
$1.1$
, and
$5.0$
, whereas left and right columns indicate the
$\lambda$
and
$\lambda '$
data, respectively. The text on the bottom right of each plot indicates the redshift z and the spearman rank coefficient
${r_s}$
.

Figure 2 reveals that there is minimal correlation between spin and halo mass, as indicated by the low Spearman coefficients. The correlation is weakly negative, except for
$\lambda$
at
$z = 0$
, and tends to strengthen with increasing redshift. This trend is consistent for
$\lambda$
, while for
$\lambda'$
, the strongest correlation appears at
$z = 1.1$
. As noted by Bett et al. (Reference Bett, Eke, Frenk, Jenkins, Helly and Navarro2007), the small magnitude and substantial scatter in
$\lambda$
and
$\lambda'$
makes this correlation difficult to detect. Our results are consistent with Macciò et al. (Reference Macciò, Dutton and van den Bosch2008), Knebe & Power (Reference Knebe and Power2008), Muñoz-Cuartas et al. (Reference Muñoz-Cuartas, Macciò, Gottlöber and Dutton2011), Ishiyama et al. (Reference Ishiyama2013), which gives us confidence that any redshift evolution we observe in the mean value of the spin parameter is not driven by a correlation between spin and mass of the underlying halo population.
2.2.3. Spin distribution
There is a rich body of literature (e.g. Barnes & Efstathiou Reference Barnes and Efstathiou1987; Bullock et al. Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001; Vitvitska et al. Reference Vitvitska, Klypin, Kravtsov, Wechsler, Primack and Bullock2002; Maller et al. Reference Maller, Dekel and Somerville2002; Hetznecker & Burkert Reference Hetznecker and Burkert2006; Bett et al. Reference Bett, Eke, Frenk, Jenkins, Helly and Navarro2007; Knebe & Power Reference Knebe and Power2008; Antonuccio-Delogu et al. Reference Antonuccio-Delogu, Dobrotka, Becciani, Cielo, Giocoli, Macciò and Romeo-Veloná2010; Davis & Natarajan Reference Davis and Natarajan2010; Muñoz-Cuartas et al. Reference Muñoz-Cuartas, Macciò, Gottlöber and Dutton2011; Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014; Zjupa & Springel Reference Zjupa and Springel2017) that has demonstrated that the spin probability distribution is well fit by a log-normal function, which can be understood physically as a consequence of the stochastic nature of angular momentum growth.Footnote
e
We verify this in Figure 3, in which we show the probability distribution of
$\ln(\lambda)$
(top panel) and
$\ln(\lambda')$
(bottom panel) at redshifts
$z=0.0$
,
$1.1$
, and
$5.0$
(black, red, and magenta, respectively).
Probability distribution of
$\ln (\lambda)$
(top panel) and
$\ln \lambda'$
(bottom panel) for our
$N=256^3$
B-suite of simulations at redshifts
${z} = 0.0$
(black filled squares),
${z} = 1.1$
(red filled squares), and
${z} = 5.0$
(magenta filled squares), all using Poisson uncertainty. The solid, dashed, and dotted curves are Gaussian fits to the distributions at each redshift discussed respectively, coloured according to their associated dataset. The coloured vertical lines show the mean values of their same-coloured dataset at each redshift.

Here the uncertainties are calculated using Poisson statistics. Both panels also show Gaussian fits to each dataset as continuous, dashed and dotted curves, respectively (colour-coded by redshift). Vertical lines indicate the mean value of the data at each redshift, for comparison with the Gaussian fits. In all subsequent analysis we use the mean and standard deviation as derived directly from the data rather than from the Gaussian fits.
Figure 3 confirms that both
$\ln(\lambda)$
and
$\ln(\lambda')$
are well fit by Gaussians, especially at low redshifts. There are small offsets between mean values measured from the data (marked by the vertical lines) and the values at which the Gaussians peak, which we would expect to be smaller for larger halo samples. At higher redshifts, the larger offsets could be an effect of the binning or an intrinsic property of the data itself, but there is still very good consistency between the data and the fits even in this regime.
3. The redshift evolution of the spin parameter
In this section, we present our results for the redshift evolution of both the mean (Section 3.1) and standard deviation (Section 3.2) of
$\ln(\lambda)$
and
$\ln(\lambda')$
. We define the probability distribution of the mean by,
\begin{equation} \mathcal{L}^{(\prime)}(z) = \frac{1}{N_h^{f}(z)} \sum_{i=1}^{N_h^{f}(z)} \ln \left[\lambda_i ^{(\prime)} (z) \right],\end{equation}
and of the standard deviation by,
\begin{equation} \mathcal{S}^{(\prime)}(z) = \sqrt{\frac{1}{N_h^{f}(z)-1} \sum_{i=1}^{N_h^{f}(z)} \left| \mathcal{L}^{(\prime)}(z)-\ln \left[ \lambda_i^{(\prime)}(z) \right] \right|^2}.\end{equation}
Here the ‘primed’ notation indicates whether we are dealing with
$\lambda$
or
$\lambda'$
, and
$N_h^{f}(z)$
is the number of field halos at a given redshift. We also demonstrate that our convergence criterion of 500 is sufficient.
3.1. Mean value
$\mathcal{L}$
Redshift evolution: Figure 4 shows
$\mathcal{L}^{\prime}(z)$
(upper set of curves) and
$\mathcal{L}(z)$
(lower set of curves) plotted against redshift, z. The different colours correspond to different minimum number of particles per halo, with
$N_h^{f}(z) \in [100, 1\,000]$
(see caption for further details).
Redshift evolution of
$\mathcal{L}(z)$
and
$\mathcal{L}'(z)$
, for our
$N=256^3$
B-suite of simulations using different
$N_{p}^\mathrm{min}$
filtering in the range
$[100,1\,000]$
(blue, magenta, and cyan solid curves, and orange, black, and red dashed curves; see legend). The Poisson uncertainty of the
$N_p^\mathrm{min}=500$
data is shown as shaded grey regions, and Equations (6) and (7) fit to
$\mathcal{L}(z)$
and
$\mathcal{L}'(z)$
respectively are plotted as solid black curves. The Figure includes two separate datasets displaying the same style but different behaviour: the bottom set corresponds to
$\mathcal{L}(z)$
, whereas the top set corresponds to
$\mathcal{L'}(z)$
. Using only host halos.

The uncertainties on
$\mathcal{L}(z)$
, which we compute as
are shown as shaded grey regions, only for the datasets with
$N_p^{\text{min}} = 500$
. For these same datasets, we also include fits to analytical functions – described below – shown as a continuous thick black line.
We observe a striking contrast between the evolution of the mean values estimated from the logarithm of the spin parameter definitions of
$\ln(\lambda)$
and
$\ln(\lambda')$
. The mean
$\ln(\lambda)$
increases approximately linearly with decreasing redshift – consistent with previous findings (e.g. Hetznecker & Burkert Reference Hetznecker and Burkert2006) – while the mean
$\ln(\lambda')$
increases gradually with decreasing redshift to
$z \sim 2$
, after which it begins to decline – which is inconsistent with Hetznecker & Burkert (Reference Hetznecker and Burkert2006). However, as we demonstrate in Appendix B, this inconsistency in
$\lambda'$
is an artefact of mass resolution and simulation volume. If we perform simulations that match those of Hetznecker & Burkert (Reference Hetznecker and Burkert2006), then we recover their results. This gives us confidence that we understand this discrepancy, and that our choice of a minimum number threshold of 500 particles is sufficient for convergence.
We note that Ahn et al. (Reference Ahn, Kim, Shin, Kim and Choi2014) have studied the redshift evolution of the Peebles (Reference Peebles1969) and Bullock et al. (Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001) spins, as well as a modified form (an empirical, redshift-dependent, correction) and find that spin decreases with increasing redshift (cf. their Figure 12). These authors track the evolution of
$\lambda$
and
$\lambda'$
rather than
$\ln(\lambda)$
and
$\ln(\lambda')$
, and so a like-for-like comparison is not straightforward, but it is noteworthy that they report a steady monotonic increase with decreasing redshift in
$\lambda$
, whereas
$\lambda'$
increases monotonically with decreasing redshift before it plateaus at
$z\simeq1$
at lower halo masses. This is broadly consistent with our results in Figure 4.
Fitting functions: The fitting functions used to model the redshift evolution of the mean logarithmic spin parameter are shown in Figure 4 as solid black curves. These fits were applied to the
$N_p^\mathrm{min} = 500$
dataset (represented by the dashed black curves), ensuring the stability and reliability of the derived parameters.
For the mean Peebles (Reference Peebles1969) spin parameter,
$\mathcal{L}(z)$
, we adopt a simple linear model:
where l is the sole free parameter, and
$\mathcal{L}(0)$
is fixed by the simulation data. This choice reflects our primary interest in the evolutionary trend of the spin parameter rather than its absolute normalisation.
In contrast, the mean Bullock et al. (Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001) spin parameter,
$\mathcal{L}'(z)$
exhibits a non-monotonic behaviour, with a turnover at lower redshifts,
$z\simeq1$
. To capture this feature, we extend the linear model by introducing a logarithmic correction term:
where l
′, a, and b are free parameters. The logarithmic correction term in Equation (7) captures the turnover behaviour of
$\lambda'$
at low redshift, allowing the model to recover the simulation-derived value at z = 0 and accurately trace its non-linear evolution.
Both fitting functions are well matched to the simulation data, as shown in Figure 4. The corresponding fit parameters are listed in Table 2. Note that the validity of these fits has been tested only within the redshift range
$0.0 \leq z \leq 5.0$
and for halo masses in the range
$9.9 \leq \log_{10}(M / M_{\odot} h^{-1}) \leq 15.3$
. There are relatively large uncertainties in some of the parameters, but this could be reduced with a larger halo sample.
The absolute change in
$\mathcal{L}(z)$
is modest compared to the intrinsic scatter (
$\approx$
$0.5$
), but it is measured with high statistical precision. The cumulative shift between
$z=5$
and
$z=0$
is
$\approx$
$0.3$
in
$\ln\lambda$
, or about 0.7 of a standard deviation. We describe the evolution as measurable and statistically robust rather than large in amplitude.
Sensitivity of results to halo definition: We have assessed the extent to which these results are sensitive to mass resolution and our choice of halo definition. We have already noted above the effects of mass resolution in our comparison to Hetznecker & Burkert (Reference Hetznecker and Burkert2006) (see Appendix B). Our choice of halo definition could lead to potential pseudo-evolution, which arises because halos are defined relative to a background density which is a function of redshift (cf. Diemer, More, & Kravtsov Reference Diemer, More and Kravtsov2013). We investigate this by using
$M_{\mathrm{vir}}$
rather than
$M_{200}$
as halo mass definition (see Appendix A) and show that the qualitative trends are consistent even though the precise values of the fit parameters differ. These tests confirm that our fitting functions reliably describe the spin evolution, although the fit parameters themselves depend on halo definition.
3.2. Standard deviation
$\mathcal{S}$
Redshift evolution: Figure 5 shows how the standard deviations in
$\ln(\lambda')$
,
$\mathcal{S}^{\prime}(z)$
(upper curves), and
$\ln(\lambda)$
,
$\mathcal{S}(z)$
(lower curves), vary with redshift, z for field halos with more than 500 particles. The shaded regions correspond to the uncertainties, which we calculate as
\[\frac{\mathcal{S}^{(\prime)}(z)}{\sqrt{2 \left(N_h^{f}(z)-1\right)}} .\]
These uncertainties represent sampling errors on the measured mean and dispersion of
$\ln\lambda$
and
$\ln\lambda^\prime$
arising from finite halo counts.
Redshift evolution of
$\mathcal{S}(z)$
(orange curve) and
$\mathcal{S'}(z)$
(blue curve) respectively, for our
$N=256^3$
B-suite of simulations, with Poisson uncertainty marked by shaded regions correspondingly coloured, and Equation (6) fit to the data as black solid curves. Using only host halos with
$N_{p} \geq 500$
.

The behaviour of
$\mathcal{S}(z)$
and
$\mathcal{S}'(z)$
are similar, both showing a weak linear increase with decreasing redshift. The normalisation of
$\mathcal{S}'(z)$
is higher and it has a steeper slope, compared to
$\mathcal{S}(z)$
. This is consistent with (Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014), who reported a positive correlation between the average spin and its dispersion.
Fitting functions: The fitting functions used to model the redshift evolution of the standard deviations in the logarithmic spin parameters are shown in Figure 5 as solid black curves. We use the same functional form as in the case of
$\mathcal{L}(z)$
(Equation 6) to describe
$\mathcal{S}(z)$
and
$\mathcal{S'}(z)$
,
where
$s^{\,(\prime)}$
is a free parameter of the model. We note that we have experimented with fitting functions that incorporate a degree of curvature, but we find that the linear fit is sufficiently accurate to not warrant the greater complexity and number of free parameters. Fit parameters are given in Table 2. As for the fits to
$\mathcal{L}(z)$
and
$\mathcal{L'}(z)$
, so too for
$\mathcal{S}(z)$
and
$\mathcal{S'}(z)$
– the validity of these fits has been tested only within the redshift range
$0.0 \leq z \leq 5.0$
and for halo masses in the range
$9.9 \leq \log_{10}(M / M_{\odot} h^{-1}) \leq 15.3$
. Parameter uncertainties could be reduced with a larger halo sample.
4. Discussion and conclusions
This study has investigated the redshift evolution of the dark matter halo spin parameter in a
$\Lambda$
CDM cosmology, using a suite of N-body simulations with varying box sizes and a consistent halo selection criterion (
$\geq$
500 particles, excluding substructure). We studied the two most commonly used spin parameters in the literature – those of Peebles (Reference Peebles1969),
$\lambda$
, and Bullock et al. (Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001),
$\lambda'$
– and provided a detailed statistical characterisation of their behaviour over the redshift range
$0 \leq z \lesssim 5$
.
As a first step, we verified that there is minimal correlation between spin and halo mass, and that the spin distributions at all redshifts are well described by log-normal functions, in agreement with previous studies. This was important to ensure that the observed redshift evolution is not a byproduct of mass-dependent selection effects and to validate our use of Gaussian statistics in our fitting of the data.
Our main results can be summarised as follows.
-
1. The mean of
$\ln(\lambda)$
increases approximately linearly with redshift, while
$\ln(\lambda')$
shows non-monotonic behaviour – increasing from
$z=5$
to
$z\sim1-2$
before declining to
$z=0$
. These trends are robust across simulation volumes and halo definitions. We have compared our results with those of earlier studies (e.g. Hetznecker & Burkert Reference Hetznecker and Burkert2006; Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014) and demonstrated that, where there are inconsistencies (e.g. in evolution of
$\lambda'$
), accounting for mass resolution can reconcile differences. -
2. The dispersion of spin values evolve with redshift. Both
$\ln(\lambda)$
and
$\ln(\lambda')$
show increasing standard deviation over time, as we would expect given the growing complexity of halo assembly histories due to mergers and accretion. However,
$\ln(\lambda')$
consistently displays higher dispersion, suggesting it is more sensitive to dynamical processes, albeit more stable against the halo edge definition as discussed in Section 1.
The contrasting evolutionary behaviours of
$\lambda$
and
$\lambda^\prime$
can be understood physically by examining how their definitions respond to the evolving internal structure of halos. The Peebles
$\lambda$
involves the total energy
$|E|$
, which is sensitive to the gravitational potential and hence to the halo concentration c. The Bullock
$\lambda^\prime$
instead normalises by the circular velocity
$V=\sqrt{GM/R}$
at the halo boundary, making it less sensitive to the internal mass distribution. Following Mo et al. (Reference Mo, Mao and White1998), the two definitions are approximately related by
$\lambda / \lambda' \approx \sqrt{2}(M_{\mathrm{vir}}/M_{200})^{1/3} f(c)$
, where f(c) encodes the concentration dependence. Since halo concentrations are known to increase with decreasing redshift (Macciò et al. Reference Macciò, Dutton and van den Bosch2008; Muñoz-Cuartas et al. Reference Muñoz-Cuartas, Macciò, Gottlöber and Dutton2011), this predicts that
$\lambda$
and
$\lambda^\prime$
should evolve differently – with
$\lambda$
more sensitive to the secular growth of concentration over time. This is consistent with our finding that
$\mathcal{L}(z)$
evolves linearly while
$\mathcal{L}'(z)$
shows a non-monotonic turnover near
$z \simeq 1$
, where concentration growth slows. The results of Appendix A support this interpretation – switching from
$R_{200}$
to
$R_{\mathrm{vir}}$
alters the effective overdensity threshold and hence the inferred concentration and changes the quantitative fit parameters while leaving the qualitative trends intact; this indicates that the differential evolution is physical rather than a numerical artefact.
We provide three fitting functions (Equations 6–8) that accurately describe the redshift evolution of the mean and standard deviation of
$\lambda$
and
$\lambda^\prime$
. Our fitting functions offer a practical, physically motivated, and computationally inexpensive way to include a redshift-dependent spin prescription into semi-empirical (e.g. Behroozi, Hearin, & Moster Reference Behroozi, Hearin and Moster2022) and SAMs (e.g. Lagos et al. Reference Lagos2018), as well as general theoretical modelling. We note that our N-body simulations, and consequently our fitting functions, do not account for the presence of baryons and galaxy formation processes (e.g. radiative cooling) in halo assembly, and so the trends we have measured will be modified when baryons are accounted for, as in hydrodynamical galaxy formation simulations. Zjupa & Springel (Reference Zjupa and Springel2017) have demonstrated that there are strong trends with mass and redshift in the spins of gas and stars within halos, which also impacts the inner dynamical structure of the halos. For this reason, we caution that the fitting functions we provide apply only to the dark matter component of the halo spin, and that additional modelling of baryonic effects will be required before they can be fully applied to the baryonic components in SAMs.
However, this does not diminish the practical value of our fitting functions for several reasons. First, even if individual halo and galaxy spins are only weakly correlated, halo spin enters SAMs through the specific angular momentum budget available for disc formation, setting the predicted disc scale lengths and structural properties. Second, the weak correlation measured in full-physics simulations applies primarily at low redshift, where our results show the evolutionary offset is modest; at
$z \gt 2$
, where disc formation and angular momentum transfer mechanisms are less well understood, the halo-galaxy spin connection remains poorly constrained. Third, regardless of the strength of the halo-galaxy spin connection, providing accurate, resolution-validated evolutionary tracks for a fundamental halo property – where none previously existed – is a worthwhile contribution to the community. Our fitting functions impose no additional computational cost and can trivially replace the common assumption of a constant spin distribution.
From a theoretical perspective, the linear growth and lower dispersion of
$\ln(\lambda)$
suggest it may serve as a more stable tracer of halo angular momentum evolution. However, the more complex evolution of
$\ln(\lambda')$
and its larger dispersion means that trends in spin should show more variation and greater sensitivity to halo dynamics, and it should be more accessible to observations, given its definition (dependence on V in Equation 2 rather than E in Equation 1); see, for instance, (Obreja et al. Reference Obreja, Buck and Macciò2022)
The differences between the Peebles and Bullock spin parameters,
$\lambda$
and
$\lambda'$
, though explored in only a limited number of previous studies (Hetznecker & Burkert Reference Hetznecker and Burkert2006; Ahn et al. Reference Ahn, Kim, Shin, Kim and Choi2014), are shown in this work to be both significant and persistent across redshift. While both are proxies for halo rotational support, they encapsulate fundamentally different physical assumptions – particularly in their treatment of energy and angular momentum. These differences manifest in distinct evolutionary trends and statistical properties, reinforcing the need for careful consideration when applying spin parameters in theoretical and observational contexts. Future work will aim to disentangle the contributions of these assumptions and assess their impact on galaxy formation models.
Acknowledgements
We thank the anonymous referee for thoughtful feedback that has helped to refine and sharpen the presentation of our results. Part of this work was performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR programme receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government.
Data availability statement
All data used in our analysis are available upon request.
Funding statement
AK is supported by the Ministerio de Ciencia e Innovación (MICINN) under research grants PID2021-122603NB-C21 and PID2024-156100NB-C21. He further thanks Fugazi for ‘waiting room’. CP acknowledges the support of the Australian Research Council (ARC) Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. RMP acknowledges the support of the Australian Government through the Australian Research Council Centre of Excellence for Dark Matter Particle Physics (CDMPP, CE200100008). AU acknowledges the support of an Australian Government Research Training Program Scholarship.
Competing interests
None.
Appendix A. Virial Radius Results
To study the effect of halo definition on our results, we have repeated the analysis of Section 3 using the virial radius,
$R_{\mathrm{vir}}$
, instead of
$R_{200}$
as our halo boundary. Here we assume a virial overdensity that is redshift dependent (cf. Bryan & Norman Reference Bryan and Norman1998),
where
where the reference density is the background density of the universe,
$\rho_{\mathrm{ref}}(z) = \Omega(z) \rho_{\mathrm{crit}}(z)$
. The masses are therefore determined by:
We expect spin values to be affected by this change in halo definition, but for the effect to be small. This should provide a test of stability of our results because they should not depend strongly on how the halo boundary is defined.
Our main results were the redshift and
$N_p^{\min}$
dependence of
$\mathcal{L}^{(\prime)}(z)$
, and the redshift dependence of
$\mathcal{S}^{(\prime)}(z)$
. In Figure A1 we show the evolution of
$\mathcal{L}(z)$
and
$\mathcal{L'}(z)$
(upper panel) and
$\mathcal{S}(z)$
and
$\mathcal{S'}(z)$
(lower panel), following the conventions of Figures 4 and 5. We provide all fitting parameters resulting from our analysis in Table A1.
Analog to Table 1 in the case of
$R_{\mathrm{vir}}$
.

There is very good qualitative agreement between these results for
$R_{\text{vir}}$
and those for
$R_{200}$
. The
$\lambda$
values tend to be slightly larger and the curves are weakly convex, whereas the
$\lambda '$
values tend to be slightly lower with no change in shape. The fitting functions provided by Equations (6) and (7) work well for both spin definitions. There now exists a small deviation between the
$\lambda$
values and the linear fit, although it is within the data’s uncertainty.
In the bottom panel we see the effect of using
$R_{\mathrm{vir}}$
is mild on
$\mathcal{S}^{(\prime)}(z)$
. The quantitative values are very similar, with the largest difference occurring for
$\lambda$
, for which the evolution is now very weak.
If we compare the fitting parameters of Tables 2 and A1, we observe the differences are small. The largest differences occur for
$\mathcal{S}^{(\prime)}(z)$
, however the general magnitude of the values is within a reasonable range. Our results are therefore robust to the choice for the halo boundary.
Appendix B. Reproduction of Hetznecker & Burkert Results
One of the earliest comparative studies of the spin parameters defined by Peebles (Reference Peebles1969) and Bullock et al. (Reference Bullock, Dekel, Kolatt, Kravtsov, Klypin, Porciani and Primack2001) was conducted by Hetznecker & Burkert (Reference Hetznecker and Burkert2006). Given the relevance of their findings, it was important for us to verify whether our simulations could reproduce their results before proceeding with our own analysis. Hetznecker & Burkert (Reference Hetznecker and Burkert2006) employed an ensemble of 82 cosmological N-body simulations with differing realisations of the initial power spectrum;
$N=64^3$
particles;
$B=20\,\mathrm{Mpc}\,h^{-1}$
; and cosmological values
$\Omega_0=0.3$
,
$\Omega_{\Lambda}=0.7$
,
$h=65$
, and
$\sigma_8=1.13$
. Their analysis was restricted to halos containing more than 500 particles.
To replicate their setup, we ran two additional low-resolution simulation ensembles using
$B=20\, \mathrm{Mpc}\, h^{-1}$
and the cosmology described in Section 2. The first ensemble consisted of 40 runs with differing realisations of the initial power spectrum with
$N=64^3$
particles, and the second comprised 10 different realisations with
$N=128^3$
. Our goal was to reproduce the results of Hetznecker & Burkert (Reference Hetznecker and Burkert2006) using the lowest-resolution set and to observe convergence towards the results presented in the main body of this work as resolution increased. For this purpose, we treat our B-suite simulations as two distinct cases: the
$B=20\,\mathrm{Mpc}\, h^{-1}$
simulation in isolation, and the full B-suite.
Redshift evolution of the mean
$\lambda$
(top panel) and
$\lambda '$
(bottom panel), normalised by their values at redshift
$z=0.0$
, for our
$B=20 \, \mathrm{Mpc} \, h^{-1}$
$N=64^3$
40-seed (blue solid curve),
$128^3$
10-seed (red dashed curve), and
$256^3$
(cyan dotted curve) simulation sets, and our
$N=256^3$
B-suite (purple solid curve). Using only host halos with
$N_{p} \geq 500$
. The Hetznecker & Burkert (Reference Hetznecker and Burkert2006) data is also shown as black filled squares.

Figures B1 shows the redshift evolution of the mean spin parameters
$\lambda$
(top panel) and
$\lambda'$
(bottom panel), normalised by their respective
$z=0.0$
values. The figure includes results from the
$N=64^3$
ensemble (solid blue curve), the
$N=128^3$
ensemble (dashed red curve), a single
$N=256^3$
particle
$B=20\,\mathrm{Mpc}\, h^{-1}$
simulation (cyan dotted curve), and the full
$N=256^3$
particle B-suite (solid purple curve). For reference, the original Hetznecker & Burkert (Reference Hetznecker and Burkert2006) data are shown as black filled squares.
Note that the
$N=256^3$
B-suite results correspond to the
$N_p^{\mathrm{min}}=500$
threshold introduced in Section 3.1 and illustrated in Figure 4. The agreement between our
$N=64^3$
ensemble and the Hetznecker & Burkert (Reference Hetznecker and Burkert2006) data is expected and reassuring.
We find that increasing resolution affects the behaviour of the spin parameters. The effect is modest for
$\lambda$
, and more pronounced for
$\lambda'$
, which we find to be more sensitive to both particle number N and box size B. This is consistent with its dependence on halo boundary definitions and dynamical state. As shown by Davis & Natarajan (Reference Davis and Natarajan2010), beyond
$N=256^3$
particles, the spin probability distribution becomes largely insensitive to further increases in resolution, suggesting that our results are robust.
An aspect not widely explored in the literature (but see Power & Knebe Reference Power and Knebe2006) is the influence of simulation volume. As seen in Figure B1, box size B has a noticeable impact on the evolution of
$\lambda'$
in particular. This warrants further investigation in future work. Overall, we successfully reproduce the behaviour reported by Hetznecker & Burkert (Reference Hetznecker and Burkert2006) and demonstrate convergence towards our main results with increasing resolution and simulation volume.






































































