What sets the splashback radius of dark matter haloes: accretion history or other properties? Free

ABSTRACT

The density profiles of dark matter haloes contain rich information about their growth history and physical properties. One particularly interesting region is the splashback radius, R_sp, which marks the transition between particles orbiting in the halo and particles undergoing first infall. While the dependence of R_sp on the recent accretion rate is well-established and theoretically expected, it is not clear exactly what parts of the accretion history R_sp responds to, and what other halo properties might additionally influence its position. We comprehensively investigate these questions by correlating the dynamically measured splashback radii of a large set of simulated haloes with their individual growth histories as well as their structural, dynamical, and environmental properties. We find that R_sp is sensitive to the accretion over one crossing time but largely insensitive to the prior history (in contrast to concentration, which probes earlier epochs). All secondary correlations are much weaker, but we discern a relatively higher R_sp in less massive, older, more elliptical, and more tidally deformed haloes. Despite these minor influences, we conclude that the splashback radius is a clean indicator of a halo’s growth over the past dynamical time. We predict that the magnitude gap should be a promising observable indicator of a halo’s accretion rate and splashback radius.

1 INTRODUCTION

In the hierarchical picture of structure formation, dark matter haloes grow by accreting a mixture of smaller haloes and smooth matter (depending on the nature of the dark matter particles; e.g. Lacey & Cole 1993). We observe indirect evidence for this process in the form of rising galaxy luminosity functions (e.g. Muzzin et al. 2013), higher clustering amplitudes at low redshift (e.g. Reid et al. 2012), and the average infall velocities around haloes (redshift-space distortions, e.g. Percival & White 2009) – but the actual process of accretion on to individual haloes remains invisible.

However, dark matter haloes carry an inevitable dynamical signature of their accretion: a caustic that forms at the apocentre of particles on their first orbit after infall. In idealized, spherically symmetric shell models, the caustic manifests as an infinitely sharp drop in density whose position depends only on the mass accretion rate (MAR) (Fillmore & Goldreich 1984; Bertschinger 1985; Shi 2016). The ‘splashback radius’, R_sp, has also been found in simulations (Adhikari, Dalal & Chamberlain 2014; Diemer & Kravtsov 2014; More, Diemer & Kravtsov 2015), where it is smoothed out due to the complexities of halo structure, as well as in satellite density and weak lensing observations (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019; Zürcher & More 2019; Murata et al. 2020; Adhikari et al. 2021; Shin et al. 2021). These detections open the possibility of constraining halo properties (such as the accretion rate) from the position of the splashback feature, or from the density profiles in general (Xhakaj et al. 2020, 2022).

Before we can attempt to infer accretion rates from such observations, we need to understand which of many possible halo properties besides the accretion rate might affect the splashback radius: the full mass accretion history (MAH) is more complicated than the power-law growth envisioned in spherical collapse models (Wechsler et al. 2002; Tasitsiomi et al. 2004; McBride, Fakhouri & Ma 2009); dark energy provides an effective force that pulls particles away from the halo centre; non-sphericity and tidal forces from nearby neighbours might change particle orbits; and the orbits are not purely radial as envision in one-dimensional shell models. Given these complications, it is not clear yet how ‘clean’ a probe of the accretion rate R_sp really is, and what definition of the accretion rate it is most sensitive to. Past work has been restricted to the dependencies of R_sp on halo mass, redshift, and cosmology via Ω_m (Diemer et al. 2017; Diemer 2020), as well as a relatively narrow range of definitions of the accretion rate (Diemer & Kravtsov 2014; Lau et al. 2015).

In this paper, we investigate a much larger array of possible correlations, including the full MAH, different definitions of the recent accretion rate, structural parameters such as concentration and ellipticity, and dynamical parameters such as angular momentum. We focus on the normalized splashback radius to scale out the absolute halo mass and size, X_sp ≡ R_sp/R_200m, where R_200m encloses an average overdensity of 200 ρ_m(z). The splashback mass generally follows similar, but weaker, trends compared to R_sp (Diemer et al. 2017). For individual haloes, R_sp cannot be measured from the noisy density profiles but must instead be inferred from the 3D density structure (Mansfield, Kravtsov & Diemer 2017) or from particle dynamics (Diemer 2017). We opt for the latter because it has produced the most comprehensive splashback halo catalogues to date (Diemer 2020), and because we are interested in how halo physics influences the underlying particle orbits. For comparison, the radius where the profile slope is steepest is commonly used in the observations discussed above, but suffers from a complex and stochastic connection to the underlying dynamics (O’Neil et al. 2021; Wang, Wang & Mo 2022). We limit the scope of our investigation to realistic ΛCDM cosmologies, although R_sp also exhibits interesting dependencies on the nature of dark matter and dark energy (Adhikari et al. 2018; Contigiani, Vardanyan & Silvestri 2019; Banerjee et al. 2020). Furthermore, we do not concern ourselves with whether the halo properties in question are observable – most of them are not, at least not directly. However, we do show that R_sp primarily responds to the recent MAR, which means that the dependencies on other quantities can mostly be neglected.

The paper is organized as follows. In Section 2, we describe the simulations and methods. We analyse the correlations between the splashback radius and MAH in Sections 3, and those with other halo properties in Section 4. We close with a conclusion in Section 5. In Appendix  A, we show that the exact definition of R_sp does not affect our conclusions. In Appendix  B, we present a PCA of mass accretion histories.

2 SIMULATIONS AND METHODOLOGY

2.1 N-body simulations and halo finding

We use a subset of the Erebossuite of dissipationless N-body simulations (Table 1). This suite contains simulations of different box sizes of two ΛCDM cosmologies; we focus on the first, which is the same as that of the Bolshoi simulation (Klypin, Trujillo-Gomez & Primack 2011) and is consistent with WMAP7 (Komatsu et al. 2011), namely, flat ΛCDM with Ω_m = 0.27, Ω_b = 0.0469, σ₈ = 0.82, and n_s = 0.95. The second is a Planck-like cosmology with Ω_m = 0.32, Ω_b = 0.0491, h = 0.67, σ₈ = 0.834, and n_s = 0.9624 (Planck Collaboration et al. 2014). The power spectra that correspond to these values were generated using Camb (Lewis, Challinor & Lasenby 2000), and they were translated into initial conditions using the 2LPTic code (Crocce, Pueblas & Scoccimarro 2006). The simulations were run with Gadget2 (Springel 2005).

We identify haloes and subhaloes using the phase–space friends-of-friends (Davis et al. 1985) halo finder Rockstar (Behroozi, Wechsler & Wu 2013a). The halo catalogues computed at different redshifts (Table 1) are the combined into time-connected merger trees using the Consistent-Trees code (Behroozi et al. 2013b). The halo catalogues and trees are described in detail in Diemer (2020).

2.2 Determining the splashback radius with SPARTA

While Rockstar computes most of the halo properties we investigate in this paper, our focus is the splashback radius. While this halo boundary is apparent as a sharp drop in the average density profiles of haloes, measuring it reliably for individual haloes demands more advanced algorithms that rely either on the full three-dimensional density field (Mansfield et al. 2017) or on particle dynamics. We opt for the latter and use the Sparta code, a flexible, parallel framework for the dynamical analysis of particle-based astrophysical simulations (Diemer 2017, 2020). Sparta follows each virtual particle as it falls into a halo and records the time and location of its first apocentre. It then reconstructs the halo’s splashback radius by smoothing the distribution of the particle apocentres in time and taking its mean (denoted R_{sp, mn}) or higher percentiles (e.g. |$R_{\rm sp,90~{{\ \rm per\ cent}}}$| for the 90th percentile). Similarly, mass definitions such as M_{sp, mn} and |$M_{\rm sp,90~{{\ \rm per\ cent}}}$| are computed as percentiles of the mass enclosed within the individual particle splashback events.

Besides changes in the radius and mass of host (isolated) haloes, using a different radius definition also changes which haloes we call subhaloes. The Moria extension to Sparta computes these new host-subhalo definitions and combines the results from Rockstar, Consistent-Trees, and Sparta in a single, convenient format for halo catalogues and merger trees. All data used in this paper is publicly available; we refer the reader to Diemer (2017, 2020) for details on the Sparta and Moria codes.

The different measures of the apocentre distribution provide different definitions of the splashback radius. For simplicity, we focus on |$R_{\rm sp,75~{{\ \rm per\ cent}}}$| in this work, but Fig. 1 and Appendix  A demonstrate that the correlations with other definitions are qualitatively very similar.

2.3 Halo selection

For each cosmology, we combine the boxes listed in Table 1, i.e. box sizes between 62.5 and 2000 h⁻¹Mpc for the fiducial WMAP7cosmology and a somewhat smaller range for the Planck cosmology. We focus on the 45 snapshots between z = 4.0 and z = 0.13 that are present in all WMAP7simulations. We do not compute results at z = 0 because the Sparta algorithm becomes less accurate for the final few snapshots of a simulation (Diemer et al. 2017). We combine haloes from the different boxes into a single sample per cosmology and per redshift, irrespective of which simulation box they originated from. We have confirmed that the results from the Planck cosmology are entirely consistent with those based on the WMAP7simulation set. We, thus report only results from the latter throughout the paper.

We make a number of cuts on the combined catalogues. First, we select only parent haloes using |$R_{\rm sp,75~{{\ \rm per\ cent}}}$| as the halo boundary. Second, we demand that haloes contain at least 1000 particles within R_200m (haloes more massive than M_200m = 1.7 × 10¹⁰M_⊙ for the smallest box and 5.6 × 10¹⁴M_⊙ for the largest box, respectively). Third, the Sparta algorithm must have successfully computed splashback radii; this criterion is not fulfilled for a small number of edge cases (Diemer 2020). Fourth, we exclude haloes whose total mass M_200m is larger than twice the gravitationally bound mass within the same radius. Most of these haloes are experiencing significant tidal disruption from nearby neighbours (Diemer 2022).

Finally, we remove so-called fly-by haloes from the sample: haloes that were previously a subhalo but that currently orbit outside of the halo boundary (e.g. Balogh, Navarro & Morris 2000; Mamon et al. 2004; Gill, Knebe & Gibson 2005). While the splashback radius should, in principle, contain all orbiting material, our chosen definition of |$R_{\rm sp,75~{{\ \rm per\ cent}}}$| means that we expect 25  per cent of particle apocentres to occur outside of this radius, and a similar number of subhalo apocentres as well. While the fraction of fly-by haloes is smaller for splashback definitions than for R_200m, it can still be significant (Diemer 2021). The mass accretion histories of fly-by haloes are typically strongly influenced by their time as a subhalo and would confuse our results. W,e thus exclude all haloes that were classified as a subhalo in the past and whose past host still exists as a host halo.

2.4 Definitions of halo properties

We consider a wide range of halo properties, which we express as dimensionless quantities to remove any dependencies on units or on the overall size scale of the halo. First, we normalize the splashback radius by R_200m and consider |$X_{\rm Rsp,75~{{\ \rm per\ cent}}} = R_{\rm sp,75~{{\ \rm per\ cent}}}/R_{\rm 200m}$| (and similarly for other percentile definitions, e.g. |$X_{\rm Rsp,50~{{\ \rm per\ cent}}}$|⁠). Our fiducial definition of the MAR follows Diemer (2017),

where a is the scale factor and t_dyn is the crossing time, which captures the time it takes a particle to cross a halo, t_{dyn, 200m} = 2R_200m/V_200m, where |$V_{\rm 200m} = \sqrt{GM_{\rm 200m}/R_{\rm 200m}}$|⁠. We note that t_{dyn, 200m} is independent of the halo radius and only depends on redshift and cosmology. We also investigate a number of other definitions of the MAR that are provided by Rockstar, which uses the virial mass, M_vir, instead of M_200m as well as different time intervals. The latter is the more important difference because the first-order difference in the mass is divided out (Xhakaj et al. 2019). In particular, Rockstar provides MARs calculated over timescales of t_{dyn, vir}, 0.5t_{dyn, vir}, and 100 Myr, which we denote Γ_{vir, Tdyn}, Γ_{vir, 0.5Tdyn}, and Γ_{vir, Inst}, respectively. We note that in the nomenclature of the Rockstar catalogues, a dynamical time is defined as a half-crossing time.

When quantifying halo mass, we not only wish to use a dimensionless quantity but also to remove the overall increase in halo mass with cosmic time. For this purpose, we use the peak height corresponding to M_200m,

where δ_c = 1.686 is the critical overdensity from the spherical top-hat collapse model (Gunn & Gott 1972) and D₊(z) the linear growth factor normalized to unity at z = 0. σ is the RMS density fluctuation within a sphere of radius R,

where M = (4π/3)ρ_m(z = 0)R³. Peak heights are calculated using the Colossus package (Diemer 2018).

We additionally consider the following quantities output by Rockstar: concentration, c_200m ≡ R_200m/r_s, where r_s is the scale radius determined from an NFW fit (Navarro, Frenk & White 1997); the normalized half-mass radius R_0.5M/R_200m; the strongest tidal force from any nearby halo averaged over the last half dynamical time, F_{tidal, tdyn}, in dimensionless units of R_vir/R_hill, where R_hill = D(M_vir/3M_{vir, neighbour})^1/3 and D is the distance between two haloes; the ratio of the kinetic and potential energies for bound particles, T/|U|; the halo spin parameter according to Bullock et al. (2001), λ_Bullock; ellipticity e ≡ (a − c)/(a + b + c), where a, b, and c are the principal axes of the halo dark matter ellipsoid and a > b > c; the half-mass scale a_0.5M, which captures how rapidly a halo accreted the second half of its M_vir; and the maximum circular velocity, v_max. We refer the reader to Behroozi et al. (2013a) for a detailed description of these parameters.

3 CORRELATIONS OF THE SPLASHBACK RADIUS WITH MAR AND HISTORY

In this section, we analyse the correlations of the splashback radius with a number of other halo properties (Section 3.1) and specifically focus on the well-known relation with MAR (Section 3.2) as well as the entire accretion history (Section 3.3). In Appendix  B, we present an alternative analysis of the connection to the accretion history using principal component analysis (PCA).

We use the dimensionless quantity |$X_{\rm Rsp,75~{{\ \rm per\ cent}}}$| as a stand-in for the splashback radius; we show that other definitions follow very similar trends in Appendix  A. We restrict ourselves to two redshifts snapshots, z = 0.2 and z = 2.0, which represent epochs before and after dark energy becomes important, respectively. Similarly, we present two halo samples: a ‘full sample’ of all haloes (Section 2.3) and a ‘high-mass sample’ with haloes of peak height ν > 2 at z = 0.2 (which corresponds to |$M_{\rm 200m} \lower.7ex{\,\, \buildrel\gt \over \sim \,\,}9 \times 10^{13} {\rm M}_{\odot }$|⁠).

3.1 General correlations between halo properties

Fig. 1 shows the correlation between the halo properties described in Section 2.4 for haloes at z = 0.2. We quantify the strength of the correlations using Pearson’s coefficient,

where x_i and y_i represent the ith data point of x and y. Here, ρ = 1 would indicate perfect positive correlation between x and y (red) whereas ρ = −1 would indicate perfect anti-correlation (blue). The lower triangle of Fig. 1 corresponds to the full sample and and the upper triangle to the high-mass sample (ν_200m > 2).

The normalized splashback radius (X_Rsp) is highly correlated with the variables related to the accretion history, namely, Γ_{200m, Tdyn}, Γ_{vir, Tdyn}, Γ_{vir, 0.5Tdyn}, and half-mass scale (a_0.5M). The instantaneous MAR over the last 100 Myr (Γ_{vir, Inst}) shows very low correlation to the splashback radius, as we discuss in detail in Section 3.2. X_Rsp is also noticeably correlated with concentration-related properties (c_200m and R_0.5M/R_200m) as well as T/|U|. These properties are tightly associated with the halo assembly history, as evidenced by their correlation to the MAR. The splashback radius is also moderately correlated with halo mass (v_max and ν_200m), but in Section 4, we show that most of this correlation is driven by the underlying correlation between mass an MAR, where higher-mass haloes accrete faster today. Interestingly, the splashback radius somewhat correlates with the tidal force around haloes (F_{tidal, tdyn}), especially at low mass. The spin parameter (λ_Bullock) and the ellipticity (e) also show interesting correlations, where fast-accreting haloes have higher angular momentum and are more elliptical. We discuss the residual impact of these properties on the splashback radius in Section 4.

In summary, at first sight X_Rsp correlates with all halo properties to various degrees. However, given their covariances with the MAR, these correlations may not be physically meaningful. Before removing the primary correlation, we now test which definition of the MAR is most predictive of X_Rsp.

3.2 Which definition of the accretion rate is most predictive?

When the splashback radius and its connection to MAR were first explored, it was not a priori clear which definition of the MAR would be most physically meaningful. In theoretical models with a fixed accretion rate, there is no ambiguity (Adhikari et al. 2014; Shi 2016), but in simulations we need to measure the halo growth over some time. The crossing time is a natural choice, since that is the timescale over which particles move to their first apocenter after infall. Indeed, Diemer & Kravtsov (2014) found that the MAR measured over approximately one crossing time was predictive of the density profiles, but the actual distribution of particles’ orbiting times is complex (Diemer 2017), and the optimal timescale has yet to be investigated in detail. Moreover, we can vary the spherical overdensity mass definition used, although this choice has a smaller impact (Xhakaj et al. 2019).

To visualize the relationship between the splashback radius and different definitions of the accretion rate, Fig. 2 shows |$X_{\rm Rsp,75~{{\ \rm per\ cent}}}$| as a function of Γ_{200m, Tdyn} (the change in M_200m over one halo crossing time) and the instantaneous rate measured over 100 Myr. As in all previous studies of the splashback radius, we find that the X_sp is a decreasing function of MAR (blue line and shaded area) and larger at higher redshift (z = 2.0, red line; e.g. More et al. 2015). The correlation is also tighter at high z (ρ = −0.646 versus −0.483), but that is partly due to the fixed mass limit of 1000 particles, which means that the average peak height in the high-redshift sample is higher. As evident from the yellow line and shaded area, the correlation for ν > 2 at low redshift is even stronger, ρ = −0.799. We observe a slight trend with mass, where R_sp is slightly smaller for higher ν at fixed Γ (in agreement with Diemer et al. 2017). We further explore this secondary correlation in Section 4.

The right-hand panel of Fig. 2 contrasts these results with the ‘instantaneous’ MAR measured over the past 100 Myr, Γ_{vir, Inst}. Given that the crossing time is about 5 Gyr at z ≈ 0, it is not surprising that this definition of the MAR is essentially uncorrelated with the splashback radius. However, what about intermediate time intervals? To quantify the tightness of the Γ-X_sp relation, Fig. 3 shows histograms of the residuals of |$X_{\rm Rsp,75~{{\ \rm per\ cent}}}$| with respect to the median relationships to four different definitions of the MAR (see Section 2.4). The panels show the full sample (top), the high-mass sample (middle), and the high-redshift sample (bottom). We confirm that Γ_{200m, Tdyn} is most tightly correlated to X_sp (blue histograms). Increasing the spherical overdensity threshold to M_vir (and thus shrinking the radius R_vir compared to R_200m) slightly increases the scatter (orange). Physically, this difference indicates that the mass at large radii, R_vir < r < R_200m, does influence the particle orbits and their apocentres. At z ≈ 2 the definitions become essentially identical, as M_vir ≈ M_180m.

Decreasing the time interval has a larger effect on the scatter, which grows to almost the same level as for Γ_{vir, Inst} (red) when halving the interval (green). Physically, this strong difference in scatter indicates that accretion during the first half-crossing time (from infall to pericentre) matters as much for the energetics of the particle orbits as accretion during their second half-orbit (from pericentre to apocentre). The effect of the time interval is even more pronounced in the high-mass sample. In the full sample, all definitions experience roughly equal tails towards small and large residuals, indicating that those are not caused by the MAR but by halo physics or by noise in the determination of R_sp (see also Figs 1 and 2).

In summary, we confirm a picture where splashback is most correlated with an accretion rate measured over one crossing time, and where the halo mass should be considered to a relatively large radius (R_200m). However, we have not probed even longer time intervals in the past, or whether it makes sense to summarize the accretion history in a single number. We, thus turn to the correlation of X_sp with the entire history now.

3.3 Which periods in the accretion history matter most?

Inspired by the results of the previous section, we use the dynamical time, t_{dyn, 200m}, as our unit of time (Section 2.4). The number of dynamical times between two redshifts z₁ and z₂ (with z₁ > z₂) is calculated as

We also consider the current dynamical time at z = 0.2, t[t_{dyn, 200m, z = 0.2}], which is roughly ∼4 Gyr. This definition is used to calculate Γ_{200m, Tdyn}. We calculate individual mass accretion histories by normalizing M_200m of haloes at each redshift by their maximum mass at any epoch, M_peak. In addition to this ‘mass history’, we also consider the ‘accretion history’ (ΔM/M_peak)/Δt_{dyn, 200m}, the change per time of M_200m/M_peak between adjacent snapshots.

Fig. 4 shows the mass and accretion histories from z = 0.2 (or t = 0) to z = 4 (top and third panel, respectively). Naturally, haloes gain mass on average. As expected in hierarchical structure formation, the high-mass sample (orange) forms more recently, on average accreting half of its current mass over the last ∼1.5 dynamical times. This recent growth is mirrored in the third panel, where we see the accretion rate reach up to half the peak mass per dynamical time at z ≈ 0.

In the second panel of Fig. 4, we quantify the correlation between |$X_{\rm Rsp,75~{{\ \rm per\ cent}}}$| and each epoch in the mass history using Pearson’s correlation coefficient. At first sight, it may seem counter-intuitive that this correlation is positive at all epochs, but a higher (normalized) mass at any epoch in the past means a lower MAR today, and thus a larger X_sp. Interestingly, the correlation peaks ∼1–2 dynamical times in the past, or about one current dynamical time (upper ticks), indicating that the accretion rate thereafter has a particularly strong impact on X_sp. The correlation of the high-mass sample (orange) peaks slightly later than for the full sample. While it is not clear whether this difference is statistically significant, we speculate that the spatial extent of the halo that really matters for particle dynamics is the splashback radius, which is slightly smaller on average for the faster-accreting high-mass sample. On the other hand, we cannot exclude other effects, e.g. from the haloes’ environments.

These arguments become easier to understand when considering the correlation between X_sp and the accretion history in the bottom panel of Fig. 4. The correlation is more or less flat at early times and turns negative between 1 and 2 dynamical times (or one current dynamical time) in the past. The peak in the correlation with the mass history corresponds to the zero-crossing of the accretion history correlations (dotted vertical lines). The rapid switch is particularly apparent in the high-mass sample, once again supporting our picture where the accretion during particles’ first orbit from infall to apocentre, but not before infall, sets X_sp. We also observe a weakening of the correlation towards t = 0, which indicates that accretion during the very last phase of the particle orbits comes too late to significantly influence their apocentres. It is unclear whether there is a fundamental mass dependence to these trends, or whether the high-mass sample is simply less affected by neighbouring haloes and stochastic noise in the particle orbits. Our results for the high-mass sample show that measuring the MAR over one current dynamical time is, indeed, a near-optimal estimator of the impact of the MAH on R_sp. In low-mass haloes, R_sp tends to be larger than R_200m, meaning that infalling particles have to travel slightly longer distances than implied by the R_200m definition of the MAR.

For comparison, we also consider the correlations of concentration with the mass and accretion histories (dashed lines in Fig. 4). Concentration has long been known to be tightly connected to the formation time of haloes (Navarro et al. 1997; Bullock et al. 2001; Wechsler et al. 2002; Ludlow et al. 2013). Naturally, formation time (which can be captured by the half-mass time, a_0.5M) and current accretion rate are strongly correlated (ρ = 0.85 for the high-mass sample in Fig. 1), but it has recently been shown that the density profiles of haloes retain separate information from these different epochs: the formation time sets the concentration, whereas the MAR affects the entire profile and particularly the position of the splashback radius (Diemer 2022; Lucie-Smith, Adhikari & Wechsler 2022). Fig. 4 confirms that the correlation between c_200m and accretion history (bottom panel) does not flatten like that of X_sp but instead rises monotonically going back in time (see Fig 1 of Wang et al. (2020) for a similar result). As a consequence, early-forming haloes have high concentrations, as expected. At very recent times, the correlation inverts, presumably because major mergers tend to temporarily increase the concentration (Wang et al. 2020).

In summary, we have, for the first time, analysed how the splashback radius is set by different epochs in haloes’ accretion histories, and confirmed that the late-time accretion over one crossing time is the key determinant. One caveat is that we have used averaged accretion histories and their correlations, leaving open the question of whether there are outlier haloes. In Appendix  B, we further explore the history-splashback connection for individual haloes using a PCA technique.

4 CORRELATIONS WITH OTHER HALO PROPERTIES

We have established that R_sp is primarily determined by the MAR, but Fig. 1 showed a number of other, potentially interesting correlations. In this section, we test whether these correlations are ‘real’ or purely driven by covariances with the MAR. Specifically, Fig. 5 is laid out like Fig. 2, but instead of the correlations between X_sp and MAR it shows secondary (or partial) correlations with other variables. Here, we have subtracted the median X_sp at the Γ_{200m, Tdyn} of each halo, i.e. the dominant correlation with MAR shown in Fig. 2. Based on the residual relation, we also compute the partial Pearson’s correlation coefficient, ρ_p(x, X_sp|Γ), between the residual of X_sp and another variable x. The partial correlation coefficients, ρ_p(x, X_sp|Γ), are calculated in a similar manner as the Pearson’s correlation coefficient (equation (4)), but using the residuals of the variables with respect to their primary relation to Γ.

Our first impression is that the secondary correlations are weak at best: the median relations (dashed lines) are more or less flat, and Pearson’s correlation coefficients are low. This conclusion also holds for the other halo properties that are not shown in Fig. 5, which we summarize in Fig. 6. Here, we have ordered the properties by the strength of their correlation for the full sample (blue). The high-mass sample exhibits somewhat different values (orange). The six properties shown in Fig. 5 were selected either because they are important in general or because they have an interesting, albeit small, partial correlation with X_sp. We now discuss these six halo properties in more detail (Sections 4.1–4.4), and explore whether the magnitude gap could be a relevant observable (Section 4.5).

4.1 Tidal force

We have already speculated that the environment of haloes might affect their splashback radius. Since we are most interested in forces acting on particles’ orbits, the tidal force from surrounding haloes might be a good way to capture these effects. Specifically, F_{tidal, tdyn} represents the strongest tidal force from any nearby halo averaged over the last 0.5t_{dyn, vir} in dimensionless units. The top left-hand panel of Fig. 5 shows that there is a positive correlation (ρ_p = 0.266), which would be compatible with a picture where particles experience outward forces that drive them to larger apocentres. The correlation seems to be most significant at the strongest tidal forces, and it vanishes for high-mass haloes (orange line), which tend to dominate their environment and are thus less susceptible to tidal influences. The correlation persists at higher redshift (red line).

4.2 Half-mass scale and concentration

As discussed in Section 3.3, the half-mass scale a_0.5M captures the general age of a halo as opposed to its late-time accretion, although the two are strongly correlated (Fig. 1). The age of haloes, in turn, affects their clustering, an effect known as ‘assembly bias’ (Gao, Springel & White 2005). We could thus imagine a number of mechanisms by which halo age might affect the splashback radius. The top-middle panel of Fig. 5 supports the picture from Fig. 4, where X_sp is almost exclusively set by late-time accretion but is slightly sensitive to halo growth at earlier times. The resulting trend is a slight negative correlation, where earlier-forming haloes have higher X_sp at fixed Γ. However, the scatter grows towards early a_0.5M because there is a longer time between the formation time and the last dynamical time, during which the haloes’ histories have time to diverge. The high-mass haloes (orange) with negative residuals show a similar trend as the full sample, but that trend reverses for late a_0.5M. We speculate that this counter-intuitive trend might be due to frequent major mergers in the recent past of very late-forming, massive haloes.

Given that the half-mass scale reliably predicts concentration (Section 3.3), we might expect that the mild secondary correlation of X_sp and a_0.5M would translate into a similar correlation with c_200m. However, the bottom left-hand panel of Fig. 5 indicates a non-significant correlation, meaning that almost the entire correlation observed in Fig. 1 is a product of the relationship between concentration and MAR.

4.3 Halo mass (peak height)

Since gravity is intrinsically scale-free, the dynamics of haloes, and thus their X_sp, should not change with halo mass. However, previous work found that X_sp is slightly smaller at fixed Γ for higher-mass haloes (Diemer et al. 2017; Diemer 2020). We confirm this picture in the top right-hand panel of Fig. 5, where the full sample shows a weak anti-correlation (ρ_p = −0.125) that appears to be caused by low-mass haloes (negative residuals). This trend is even slightly stronger at z = 2. The correlation almost vanishes for the high-mass sample. While surprising at first, this result seems in rough agreement with Fig. 3 in Diemer et al. (2017), which shows that the mass trend becomes very weak at ν > 2 and z = 0.

4.4 Ellipticity and kinetic energy

One uncertainty in the dynamical determination is the large spread of particle apocentres, and non-sphericity certainly contributes a significant fraction of that scatter (Diemer 2017). We could thus easily imagine that making a halo more elliptical would mean to distribute a part of its particles’ apocentres to larger radii, and thus to increase high percentiles of the distribution (such as the 75th percentile used here). While the formal correlation coefficient between X_sp and e in the bottom-centre panel of Fig. 5 is low, the median relation does exhibit such a trend. Moreover, in the high-mass sample the partial correlation coefficient is 0.181, comparable to those analysed in the previous sections. This result highlights why it is important to consider secondary correlations: overall, ellipticity is quite strongly anti-correlated with X_sp (Fig. 1) because a high e goes hand in hand with high accretion rate, and thus a low X_sp.

Another quantity that correlates strongly with both MAR and splashback radius is T/|U|, the ratio of the kinetic and potential energies of particles in a halo. As for ellipticity, however, the bottom right-hand panel of Fig. 5 shows that the secondary correlation with X_sp is positive, which we would expect given that higher kinetic energy at fixed potential energy (fixed mass) should mean larger orbits and thus larger apocentres.

4.5 The magnitude gap

While the correlations we have discussed in the previous sections tell us about the physics of halo formation, the respective quantities are generally not directly observable, at least not for individual haloes. One possible indicator of halo growth that is observable is the ‘magnitude gap’, the difference in magnitude between the brightest cluster galaxy (BCG) and its nth brightest satellite (Jones et al. 2003; D’Onghia et al. 2005; Dariush et al. 2010; Deason et al. 2013; Hearin et al. 2013; Kundert, D’Onghia & Aguerri 2017; Golden-Marx & Miller 2018; Vitorelli et al. 2018; Farahi, Ho & Trac 2020). Satellites get stripped and tidally disrupted over time, and these processes particularly affect the largest satellites, which suffer most dynamical friction (Chandrasekhar 1943). As time passes since a halo actively grew and accreted massive satellites, the brightness gap between the BCG and brightest satellite thus grows on average, meaning that a large magnitude gap is an indicator of an old halo age or a low (current) accretion rate. We would thus expect the size of the magnitude gap to positively correlate with X_sp.

To investigate this question without artificially ‘pasting’ galaxies on the haloes in our N-body simulations, we use a simple proxy for the magnitude gap: the ratio of M_200m of the parent halo and the peak mass of the most massive subhalo. Assuming that satellite magnitudes are decreasing (being brighter) with increasing subhalo mass, this quantity qualitatively reflects the the magnitude gap. Since we need the largest subhalo to be well-resolved, we apply a more stringent cut on the host halo mass, namely, 20 000, 20 000, 15 000, 10 000, and 5000 dark matter particles for the L0063, L0125, L0250, L0500, and L1000 simulations, respectively. The corresponding halo masses are 3.4 × 10¹¹, 2.8 × 10¹², 1.65 × 10¹³, 8.7 × 10¹³, and 3.5 × 10¹⁴M_⊙. We have checked that the results are converged with these limits.

Fig. 7 shows that the mass ratio does exhibit the expected positive correlation with X_sp. The size of the correlation coefficient (ρ = 0.28) is comparable to many of the correlations in Fig. 1, but it should be observable even in individual haloes (with the caveat of the translation from subhaloes to galaxies). The scatter in the relation is, however, comparable to the size of the overall trend. We have confirmed that the lower X_sp in the high-mass sample at fixed subhalo mass ratio is driven by the higher MARs in the high-mass sample. In Table 2, we show the Pearson’s correlation coefficients of the mass ratio to several representative halo properties.

In summary, given that the density profile of individual haloes is not observable with current weak lensing or satellite profiles, the magnitude gap is, to our knowledge, the most promising single observable to split haloes of fixed mass into samples that should exhibit different splashback radii. However, this connection should be checked in hydrodynamical simulations of galaxy formation, where actual magnitudes are available.

5 CONCLUSION

We have presented a thorough assessment of which properties of dark matter haloes set the (relative) position of their splashback radius. We have confirmed and refined the previous finding that the accretion history is the most important factor, but we have also discovered a number of interesting secondary effects. Our most important conclusions are as follows.

• The accretion rate over the last crossing time is by far the strongest predictor of the splashback radius, with a physically well-understood anti-correlation. Changing the mass definition from M_200m to higher overdensities such as M_vir, or reducing the time interval over which the rate is measured, increase the scatter in the relation.

• Correlating the splashback radius with the full MAH reveals that R_sp suddenly becomes sensitive to the recent accretion at one current dynamical time in the past (or 1.5 integrated dynamical times). As a result, R_sp is most positively correlated with the cumulative, fractional mass at that time. These results support a picture where particles are influenced by mass growth only between their infall and first apocentre.

• We confirm that halo concentration is also tightly related to the accretion history, but that it captures the overall age of the halo rather than its recent accretion.

• We study the effect of other halo properties on the splashback radius by considering their impact at fixed MAR. All secondary correlations are comparatively weak, but we confirm that more massive haloes have slightly smaller R_sp. An early formation time, strong tidal forces from neighbours, high ellipticity, and high kinetic-to-potential energy ratio can all slightly increase R_sp at fixed accretion rate.

• In a PCA of the mass history, the first component explains 69  per cent of the halo-to-halo variations and strongly correlates with the splashback radius. We build a polynomial fitting function based on the first three components, but this function predicts the splashback radius only slightly more accurately than a simple relation with accretion rate.

One caveat to our analysis is that Pearson’s correlation coefficient captures only the linear relationships between variables but ignores any non-linear effects. Similarly, we have neglected baryonic effects. We have focused on the particular splashback scale of haloes as determined by particle dynamics, but in future work this analysis could be extended to the full density profiles of haloes (e.g. Xhakaj et al. 2022).

ACKNOWLEDGEMENTS

TS is supported by the US Department of Energy under award DE-SC0018053. Many of the computations were run on the Midway computing cluster provided by the University of Chicago Research Computing Center and on the DeepThought2 cluster at the University of Maryland.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

References

APPENDIX A: CORRELATIONS WITH DIFFERENT SPLASHBACK DEFINITIONS

Throughout the paper, we have quantified correlations with the splashback radius as defined by the 75th percentile of particle apocentres, but our results hold for other definitions as we show in this Appendix. First, Fig. A1 repeats the z ≈ 0 panel of Fig. 2, but it shows the residuals of different percentile definitions (Section 2.2) over their median relationship with Γ_{200m, Tdyn}. The scatter monotonically increases with radius for two reasons. First, percentiles near the median (50th percentile) carry a smaller statistical error than those near the edge of the distribution (e.g. the 90th percentile). Secondly, and more importantly, the tail of the apocentre distribution is more sensitive to effects such as non-sphericity and environment (Diemer 2017). The 75th percentile represents an intermediate value with intermediate scatter, and our results for the connection between R_sp and the accretion history qualitatively hold for the other definitions.

We also need to check whether the secondary correlations found in Section 4 differ systematically for different definitions of the splashback radius. We summarize the results in Fig. A2, which shows Pearson’s correlation coefficients for different percentile definitions and various halo quantities. Overall, the correlations are very similar for different definitions, and the trends with percentile are monotonic for each halo property. Most notably, the correlation with accretion rate is even significantly stronger for the 50th percentile than for the 75th percentile that we analysed throughout the paper.

APPENDIX B: PCA

In Section 3.3, we explored the correlation between R_sp and different epochs in the average MAH. In this section, we explore whether we can explain the correlations based on a simplified version of individual haloes’ mass histories that are decomposed via a PCA. Since the mass history (first panel of Fig. 4) contain nearly identical information as the MAH (third panel), we perform the PCA only on the mass history, M(t)/M_peak. For n data points (haloes) each having m variables (epochs in their mass history), the PCA finds the vector in m-dimensional space along which the variance of the data points is maximized (the so-called ‘first principal direction’). The process is repeated for each subsequent dimension such that the new vector is orthogonal to all previous vectors. This procedure results in a vector space where the data can be exactly expressed as a series of coefficients times the component vectors, with the goal of reducing the dimensionality to the first few components (which ideally contain much of the ‘explained variance’). In other words, any halo’s mass history can be accurately approximated as a linear combination of the first few PCA components.

Fig. B1 shows the first three PCA components of the mass history, including the explained variance ratios of about 69 per cent, 12 per cent, and 7 per cent (almost 90 per cent total). Unsurprisingly, adding more components adds little to the accuracy of the PCA. The first component (red) modulates the overall shape of the mass history, whereas a high second (green) component indicates strong recent mass accretion. In contrast, the third (blue) component does not clearly indicate an excess or lack of recent accretion.

In Fig. B2, we investigate the relationship between the PCA coefficients of haloes and their relative splashback radii (top panels). We also explore how much of the correlation is mediated by the recent accretion rate by showing the correlation with Γ in the bottom panels. We quantify the strength of the correlation with Pearson’s coefficient (see legend). A high first coefficient is tightly correlated with a low accretion rate, and thus with a high splashback radius (left-hand panels of Fig. B2). The correlation is particularly clear in the high-mass sample (pink line), presumably because large haloes dominate their environment, and are thus less susceptible to other factors that might modulate R_sp. We find weaker and opposite trends in the second PCA component, where a high coefficient corresponds to a high recent accretion rate (middle panels). The third component is not predictive of recent accretion (Fig. B1) and also correlates with R_sp rather weakly (right-hand panels of Fig. B2), although it is a stronger predictor of R_sp in the high-mass sample. Overall, our results affirm once again that the majority of the effect on the splashback radius is mediated via the recent accretion over one dynamical time, and that the full mass history adds relatively little information.

We further assess the additional information content by constructing a simple fitting function that predicts the splashback radius from the PCA coefficients. In particular, we model |$X_{\rm Rsp,75 {{\ \rm per\ cent}}}$| with three second-order polynomials,

where x, y, and z are the coefficients of the first, the second, and the third principal components, and a_{x, i}, a_{y, i}, and a_{z, i} are free parameters. We determine the best-fit values of the nine parameters by minimizing χ². For comparison, we construct an even simpler second-order polynomial ‘model’ that is a function of only the accretion rate,

We quantify the success of these functions at predicting R_sp in Fig. B3, which shows the distribution of the fit residuals. While the PCA-based fit results in slightly lower average residuals than the purely accretion rate-based fit, the difference is marginal, especially given that equation (B1) contains nine free parameters and equation (B2) only three. We also compare to the fitting function of Diemer (2020) in yellow, which leads to slightly higher residuals than the PCA fit and a comparable performance equation (B2). We note that the fitting also quantifies trends with ν, Ω_m(z), and percentile, and thus sacrifices some accuracy for a limited number of free parameters.

In summary, the PCA provides an efficient way to parametrize mass histories. The first component captures |$\sim 70 {{\ \rm per\ cent}}$| of the halo-to-halo fluctuations and is tightly correlated to the splashback radius (ρ = 0.763 for the high-mass sample). This type of analysis could be applied to observational data, for example to constrain the accretion history of individual haloes based on their splashback radii (assuming those can be measured, e.g. from stellar haloes; Deason et al. 2021).