key: cord-0286582-tzecewzj
authors: Diemer, Benedikt
title: A dynamics-based density profile for dark haloes. I. Algorithm and basic results
date: 2021-12-07
journal: nan
DOI: 10.1093/mnras/stac878
sha: 485ec441feb8584b7acd0cef82bcc3577eed5f93
doc_id: 286582
cord_uid: tzecewzj

The density profiles of dark matter haloes can potentially probe dynamics, fundamental physics, and cosmology, but some of the most promising signals reside near or beyond the virial radius. While these scales have recently become observable, the profiles at large radii are still poorly understood theoretically, chiefly because the distribution of orbiting matter (the one-halo term) is partially concealed by particles falling into halos for the first time. We present an algorithm to dynamically disentangle the orbiting and infalling contributions by counting the pericentric passages of billions of simulation particles. We analyse dynamically split profiles out to 10 R200m across a wide range of halo mass, redshift, and cosmology. We show that the orbiting term experiences a sharp truncation at the edge of the orbit distribution. Its sharpness and position are mostly determined by the mass accretion rate, confirming that the entire profile shape primarily depends on halo dynamics and secondarily on mass, redshift, and cosmology. The infalling term also depends on the accretion rate for fast-accreting haloes but is mostly set by the environment for slowly accreting haloes, leading to a diverse array of shapes that does not conform to simple theoretical models. While the resulting scatter in the infalling term reaches 1 dex, the scatter in the orbiting term is only between 0.1 and 0.4 dex and almost independent of radius. We demonstrate a tight correspondence between the redshift evolution in LCDM and the slope of the matter power spectrum. Our code and data are publicly available.

Although dark matter accounts for the majority of mass in the Universe, it is notoriously hard to observe. The most promising place to at least infer some of its properties are dark matter haloes (hereafter simply haloes) because their roughly spherical nature makes them more orderly than filaments, walls, or (arguably) voids, and because the highest densities are reached at halo centres. However, those centres are also strongly affected by baryonic physics that we cannot claim to fully model at this point. Thus, it is critical to understand haloes in detail out to large radii, where they join into the surrounding large-scale structure.

The shapes of haloes are complex and not spherically symmetric, especially at their interface with the cosmic web (e.g., Bond et al. 1996; Jing & Suto 2002; Allgood et al. 2006; Mansfield et al. 2017 ). However, the full three-dimensional density structure of haloes is challenging to model theoretically (e.g., Vogelsberger et al. 2009 Vogelsberger et al. , 2011 , complex to quantify from simulations (Lilley et al. 2018) , and difficult to observe (Shin et al. 2018) . Instead, spherically averaged density profiles (hereafter abbreviated as 'profiles') provide a simple way to describe halo structure. This paper continues a long history of efforts to theoretically understand density profiles, be it via spherical collapse models (e.g., Gunn & Gott 1972) , from N-body simulations (e.g., Frenk et al. 1985 Frenk et al. , 1988 ; Barnes & Efstathiou 1987; Dubinski Email: diemer@umd.edu & Carlberg 1991), or with empirical fitting functions such as those of Einasto (1965) , Hernquist (1990) , and Navarro et al. (1997, hereafter NFW) .

This enormous body of work has established a number of basic facts. The inner profile, roughly within R 200m (the radius enclosing an average overdensity of 200 times the cosmic mean), is dominated by particles orbiting in the halo (e.g., Cuesta et al. 2008; Diemand & Kuhlen 2008) . Although often labelled the 'one-halo term', we call this component the 'orbiting term'. It gradually steepens with radius, a generic consequence of hierarchical collapse from initial peaks with a range of shapes (e.g., Huss et al. 1999; Dalal et al. 2010; Ludlow et al. 2013) . At large radii (roughly beyond R 200m ), the profiles become dominated by matter that is falling into the halo for the first time and, at even larger radii, by the statistical contribution from other haloes. These two components are often collectively mislabelled as the 'two-halo term'. We will refer to them as the 'infalling term' because that contribution dominates out to the radii we consider (10 R 200m ).

Despite this progress, two key questions remain. First, how is dark matter distributed near the halo centre, where baryons dominate? Much work has been focused on this question (Moore et al. 1999; Jing & Suto 2000; Navarro et al. 2004 Navarro et al. , 2010 Kazantzidis et al. 2006 ), but its solution demands an accurate understanding of complex baryonic physics (e.g., Pontzen & Governato 2012; Di Cintio et al. 2014; Velliscig et al. 2014; Schaller et al. 2015; Wang et al. 2020) . The second question is the focus of this paper: what is the structure near the transition between the orbiting and infalling terms, and beyond the halo edge? While feedback can eject gas to those radii, baryonic effects are much weaker near the transition region (Schneider et al. 2019; Aricò et al. 2020; Sorini et al. 2021 ) and do not fundamentally alter the features we are interested in (Lau et al. 2015; O'Neil et al. 2021) .

The halo outskirts (r > ∼ 0.5 R 200m ) have recently received renewed attention because they turn out to be sensitive to otherwise inaccessible halo properties such as the mass accretion rate (Diemer & Kravtsov 2014 , hereafter DK14; see also Xhakaj et al. 2020 Xhakaj et al. , 2021 , the nature of dark matter (Banerjee et al. 2020) , and deviations from General Relativity (Adhikari et al. 2018; Contigiani et al. 2019a ). At the same time, rapid observational progress has led to highaccuracy constraints based on profiles of satellites in individual clusters (Rines et al. 2013; Tully 2015; Patej & Loeb 2016) , in stacked clusters (More et al. 2016; Baxter et al. 2017; Nishizawa et al. 2018; Zürcher & More 2019; Murata et al. 2020) , and in weak-lensing profiles (Mandelbaum et al. 2006; Umetsu et al. 2011; Umetsu & Diemer 2017; Chang et al. 2018; Contigiani et al. 2019b; Shin et al. 2019 Shin et al. , 2021 . Deviations from simple fitting functions for the orbiting profile have long been known to bias weak lensing masses (Becker & Kravtsov 2011; Oguri & Hamana 2011) , but the recently acquired high signal-to-noise ratio in the outer profiles has made detailed theoretical modelling even more pressing.

The paramount difficulty in modelling the outer profile arises from the superposition of the orbiting and infalling terms, as illustrated in Fig. 1 . These terms (gray lines) trade off near the steepest part of the total profile (red). By construction, the steeply decreasing orbiting term leaves the infalling term to dominate at large radii, meaning that the true shape of the orbiting term near its truncation cannot be measured. In the figure, it is approximated by the DK14 fitting function, but its asymptotic functional form is essentially unconstrained. The shape near the truncation is particularly important because the orbiting profile contains valuable information about particle dynamics. For example, the edge of the apocentre distribution, also known as the splashback radius, marks a physically meaningful halo boundary (DK14, Adhikari et al. 2014 , More et al. 2015 . Approximations such as the Einasto or NFW profiles were designed to fit at r < ∼ R 200m and thus do not match the behaviour in the transition region ( Fig. 1) . Additional theoretical modelling and fitting functions have helped to elucidate the nature of the infalling term Prada et al. 2006; Tavio et al. 2008; Oguri & Hamana 2011) , but they could not quite fit the sharp drop near the transition region (DK14).

The key insight is that the profiles are intrinsically dynamical in nature. It has long been known that the orbiting profiles depend on the growth history of haloes (Navarro et al. 1996; Bullock et al. 2001; Wechsler et al. 2002; Zhao et al. 2003; Tasitsiomi et al. 2004; Kazantzidis et al. 2006; Ludlow et al. 2013) . Moreover, DK14 showed that the profile shape at all radii depends on the current mass accretion rate. While statistical mechanics has previously been the basis of theoretical models for the orbiting term (King 1966; Hjorth & Williams 2010; Pontzen & Governato 2013) , the dependence on growth rate goes beyond the velocity dispersion of particles and highlights that a full understanding of the transition region must depend on a dynamical distinction between infalling and orbiting particles. A cut in phase space cannot achieve this goal because the two populations overlap in the radial range of interest (see Figs. 2 and 3; Fukushige & Makino 2001; Cuesta et al. 2008; Diemand & Kuhlen 2008; Sugiura et al. 2020) , although approximate cuts are useful when classifying observed satellite distributions (e.g., Oman et al. 2013; Tomooka et al. 2020; Aung et al. 2021) .

In this paper, we present a novel, reliable algorithm to system- In general, we can measure only the total profile (red), which makes it difficult to discern the true shape of the orbiting term. The gray shading marks r < 0.1 R 200m , roughly 10 times the size of galaxies (Kravtsov 2013) , although baryons play an important role even at larger radii. The bottom panel shows the logarithmic slope, d ln ρ/d ln r.

For this example, we consider the typical (median) profile of a cluster halo with M 200m = 5 × 10 14 h −1 M , which corresponds to a median concentration of c 200m ≈ 6.5 (Diemer & Joyce 2019) . The gray/red profiles show the fitting function of DK14, which accurately fits the total profiles, but the behaviour of the orbiting term near its cut-off is poorly constrained. This edge represents the extent of the largest particle orbits and is also known as the splashback radius, R sp . This boundary is often operationally defined as the location where the logarithmic slope is steepest, R steep , but this radius is the result of a complicated interplay between the infalling and orbiting components. The overall shape of the profiles in the transition regions is not well matched by classical fitting functions such as the Einasto or NFW profiles (yellow).

atically perform a dynamical distinction for billions of particles. Specifically, we count the pericentric passages for each halo particle in large N-body simulations and construct separate orbiting and infalling profiles (see Bakels et al. 2021 for a similar algorithm applied to subhaloes and Sugiura et al. 2020 for an apocentre-based scheme). Our algorithm is implemented in the publicly available Sparta framework, which is described in Diemer (2017, hereafter D17) and Diemer (2020a, hereafter D20). The paper is organized as follows. In Section 2, we describe our simulations and the dynamical splitting algorithm. Additional convergence tests are shown in Appendix A. In Section 3, we present the most important insights into the separate orbiting and infalling profiles, referring the reader to online figures for a more complete accounting of the numerical data (see benediktdiemer.com/data). We discuss theoretical questions in Section 4 and summarize our find- 

10 −1 10 0 10 1 r (h −1 Mpc) Figure 2 . The dynamics and density structure of a halo with a relatively low accretion rate (Γ = 1.3, from TestSim200 at z ≈ 0). The halo's particles are split by our orbit counting algorithm into orbiting particles (which have had at least one pericentre, left column) and infalling particles (centre column). The top panels show the projected density field using an SPH-like smoothing kernel with a fixed size and a logarithmic color scale. The second row shows the corresponding phase-space density structure (on a linear color scale). The bottom row shows the corresponding density profiles. For comparison, we show R vir (solid lines) and R sp,90% (dashed lines) as measured by Sparta. For this halo, the infalling and orbiting terms are separated relatively clearly in phase space. ings in Section 5. This paper is the first of a series. The second paper (hereafter Paper II) will present improved fitting functions for the orbiting and infalling profiles. Paper III will investigate the parameter space of these functions and connect it to halo properties. Throughout the paper series, we follow the notation of D20 (see also Section 2.2). The code and profile data used for this paper are publicly available.

In this section, we give a detailed account of the numerical methods that underlie our results. We begin by introducing our simulation data in Section 2.1 and by reviewing the definitions of important halo-related quantities in Section 2.2. We briefly summarize the Sparta code in Section 2.3 but refer the reader to D17 and D20 for details. Our new algorithms consist of two main parts, namely counting the orbits of particles (Section 2.4) and constructing density profiles from those particles (Sections 2.5 and 2.6). We discuss how we select and combine halo samples in Sections 2.7 and 2.8, and we subject these methods to convergence tests in Appendix A.

We use the Erebos suite of dissipationless N-body simulations (Table 1). This suite contains simulations of different box sizes and resolutions as well as two ΛCDM cosmologies and self-similar universes. The first ΛCDM cosmology is the same as that of the Bolshoi simulation (Klypin et al. 2011 ) and is consistent with WMAP7 (Komatsu et al. 2011) , namely, flat ΛCDM with Ω m = 0.27, Ω b = 0.0469, σ 8 = 0.82, and n s = 0.95. The second is a Planck-like cosmology with Ω m = 0.32, Ω b = 0.0491, h = 0.67, σ 8 = 0.834, and n s = 0.9624 (Planck Collaboration et al. 2014 ). These two cosmologies span the currently favoured ranges of the parameters that are most important for structure formation, such as Ω m .

To investigate the universality of profiles, it is important to consider universes with very different power spectra. We rely on four Fig. 2 but for a rapidly accreting halo with Γ = 3. The phase-space structure of the orbiting particles is somewhat irregular. Moreover, the infalling term contains two massive and many smaller subhaloes. The velocity dispersion of the most massive subhalo (at a radius of about 1.3h −1 Mpc) is large enough to give some of its particles a positive total v r , highlighting why our algorithm cannot solely rely on the radial velocity as an indicator of pericentres.

self-similar Einstein-de Sitter universes with power-law initial spectra, P ∝ k n , with n = −1, −1.5, −2, and −2.5. In these simulations, length and time are scale-free, meaning that the density profiles are independent of redshift when rescaled to meaningful physical units. Their differences isolate the impact of the initial power spectrum (e.g., Efstathiou et al. 1988; Cole & Lacey 1996; Knollmann et al. 2008; .

The power spectra of the ΛCDM simulations were generated using Camb (Lewis et al. 2000) . The power spectra were translated into initial conditions using the 2LPTic code (Crocce et al. 2006) . The simulations were run with Gadget2 (Springel 2005) . We identify haloes and subhaloes using the phase-space friends-of-friends halo finder Rockstar (Behroozi et al. 2013a ) and construct merger trees using the Consistent-Trees code (Behroozi et al. 2013b ). These catalogues are described in detail in D20. The halo centre is taken to be the average position of the most bound particles, a procedure that minimizes the Poisson noise in the position (Behroozi et al. 2013a ). The velocity is taken to be the average of all particles within 0.1 R 200m .

We consider three types of halo boundary definitions: bound-only spherical overdensity (SO) radii, all-particle SO radii, and splashback radii. Given a halo radius R X , we denote the mass enclosed within it as M X and the corresponding number of particles as N X .

Bound-only radii are calculated by Rockstar, which removes gravitationally unbound particles from friends-of-friends groups and subgroups in six-dimensional phase space. SO radii are then computed by finding the outermost radius where the enclosed density falls below the given SO threshold. We mostly consider a threshold of 200 times the mean matter density of the Universe, ρ m ; we denote the corresponding boundary and enclosed mass as R 200m,bnd and M 200m,bnd . Alternatively, one can use a threshold based on the critical density of the Universe, ρ c , leading to definitions such as R 200c and M 200c . Finally, we compute the varying virial overdensity definition, R vir , using the approximation of Bryan & Norman (1998) .

However, throughout the paper we consider density profiles that include all particles, whether gravitationally bound or not. Thus, we mostly use the corresponding SO radii and masses, particularly Table 1 . The Erebos suite of N-body simulations. L denotes the box size in comoving units, N 3 the number of particles, m p the particle mass, and the force softening length in comoving units. The redshift range of each simulation is determined by the first and last redshifts z initial and z final , but snapshots were output only between z f−snap and z final . The earliest snapshots of some simulations do not yet contain any haloes, so that the first catalogue with haloes is output at z f−cat ; the Sparta and Moria data also begin at that redshift. The cosmological parameters are given in Section 2.1. Here, 'PL' indicates self-similar cosmologies with a power-law initial spectrum with slope n. The references correspond to Diemer et al. (2013b, DKM13) , Diemer & Kravtsov (2014, DK14) , Diemer & Kravtsov (2015, DK15) , and Diemer (2017, D17). Our system for choosing force resolutions is discussed in DK14. Some previous papers on the Erebos simulations (including DK15, D17, and D20) erroneously claimed to be fixed in physical units, but this error had no impact on their calculations. We consider only the density profiles of isolated (or host) haloes because those of subhaloes typically contain a large contribution from the host. We generally define subhaloes as those that lie within R 200m,bnd of another, larger halo (Behroozi et al. 2013b , D20), but we discuss the impact of other choices in Section 2.7.

When considering the question of universality, we need to fairly compare halo masses across redshifts and cosmologies. For this purpose, we convert mass to peak height, ν X , which captures the statistical significance of haloes in the linear overdensity field. It is formally defined as ν X = δ c /σ(M X ), where δ c (z) = 1.686 is the threshold overdensity in the top-hat collapse model (Gunn & Gott 1972) and σ(M X ) is the variance of the linear power spectrum measured in spheres of the Lagrangian radius of a halo, R L , the comoving radius that encloses the mass M X at the mean density of the Universe, M L = M X = (4π/3)ρ m (z = 0)R 3 L . We do not apply corrections for the finite volume of our simulations or for a slight evolution of δ c in ΛCDM cosmologies (Mo et al. 2010 ) because these effects lead to tiny changes in ν that are relevant when computing mass functions but irrelevant for our purpose of binning haloes (Diemer 2020b ). While we could compute peak height from any M X , we generally use ν 200m,all , which is denoted ν 200m or simply ν for brevity. We use the Colossus code (Diemer 2018) to compute peak heights, approximating the power spectrum with the transfer function of Eisenstein & Hu (1998) . We refer the reader to Diemer (2020b) for a much more detailed description of these calculations.

Finally, we define the logarithmic mass accretion rate per dynam-ical time following D17, using Γ as a shorthand,

Using mass definitions other than M 200m makes little difference because the overall normalization is divided out, but the definition of the dynamical time does matter (Xhakaj et al. 2019 ); we use the crossing time, (Diemer 2017) . The relatively large radius of R 200m reduces any baryonic effects on the accretion rate.

Sparta is a flexible, parallel C framework for the dynamical analysis of particle-based astrophysical simulations. Based on input from a halo finder (Rockstar and Consistent-Trees in our case), the code tracks the trajectories of haloes and their constituent particles; the trajectory relative to the halo centre is called an orbit. The user can add plugin-style modules that analyse each orbit. The aggregated results can then be further analysed on a halo-by-halo basis. For example, Sparta determines the splashback radius by tracking all particles that enter each halo and measuring their first apocentre. From the locations and times of these events, it computes the splashback radius of a given halo by smoothing the distribution of particle splashbacks in time and taking its mean (R sp,mn ) or higher percentiles (e.g., R sp,90% for the 90th percentile; D17).

Our new orbit counting algorithm is implemented in the same framework but counts pericentres instead of apocentres. The goal is to determine whether each particle in a halo is orbiting (part of the one-halo term) or falling in for the first time. We choose a simple criterion: a particle is orbiting if it has undergone at least one peri- 

The magenta dots in the top panels indicate the time and radius where the particle first crosses R 200m . The black dots show the pericentres detected by our algorithm, where we have interpolated the pericentre time to coincide with v r = 0. The first column shows a regular orbit with a roughly fixed apocentric distance. Our algorithm is not fooled by the short epoch with v r > 0 before infall because φ ≈ 1 indicates that the particle has not changed its angular direction. Conversely, the algorithm does detect an early pericentre where φ 1 in the second orbit before the particle has even crossed R 200m . The third column shows a highly tangential orbit, which never comes close to the halo centre, exhibits little variation in the radial velocity, and smoothly transitions from alignment to anti-alignment. At about 13 Gyr, our algorithm sets a lower limit flag (open circles) because φ < −0.8. In this particular case, a pericentre is also detected soon thereafter because the radial velocity does become slightly positive, but it is reassuring that the algorithm does not need to rely on this error-prone feature. The fourth column shows an orbit that is unresolved in time. While our algorithm correctly determines the first pericentre, many of the following pericentres cannot be identified in the noisy orbital data. centre, and infalling otherwise. 1 This definition reduces our task to counting pericentres in the orbits of each particle. In principle, a pericentre should manifest itself as a switch from negative to positive radial velocity with respect to the halo centre and as a minimum in radius (see example orbits in Fig. 4 ). In practice, however, a reliable algorithm needs to wrestle with three sources of confusion: noise in the trajectories, particles that fall in as part of subhaloes, and highly tangential orbits.

Before tackling these issues, we need to clarify our definition of a pericentre. Both theoretically and practically, it makes no sense to allow pericentres at arbitrarily far distances and long times before a particle enters a halo. For example, would we call a particle's motion an 'orbit' if it resides at, say, 4R 200m ? We could demand the particle to be gravitationally bound at the time of pericentre, but boundness is an ill-defined concept in cosmological simulations (see Behroozi et al. 2013c and the discussion in D20). We do not introduce a max-imum radius for pericentres because such limits lead to noticeable, unphysical features in the infalling density profile. Instead, we require that our algorithm must be able to discern between spurious and real pericentres even at large distances. In practice, however, we cannot measure pericentric radii greater than the radius where we begin to track the infall of particles, which we set to 2 R 200m in our Sparta runs (D17, D20; see Appendix A3).

We generally avoid using radius information and instead focus on the radial velocity, v r . All velocities are understood to be in physical units, meaning that they include the Hubble expansion flow. Besides a particle's velocity, we also consider its angular direction, which we quantify as the dot product of the current and initial unit position vectors, φ(t) ≡r r r(t) ·r r r(t ini ). This quantity deviates from unity as the particle's angular position moves away from its direction when it was first discovered (at 2 R 200m ) and reaches −1 at the opposite side of the halo (bottom panels of Fig. 4 ). We find that tracking particles only after they cross R 200m (rather than 2 R 200m ) can lead to a fairly different reference direction because it ignores the tangential motion that may have occurred previously (Appendix A3). As soon as a particle is being tracked, we begin looking for switches in the sign of v r from negative to positive. By itself, however, this criterion does not suffice. First, it picks up many spurious detections, e.g., due to halo finder noise in the centre positions of haloes. Second, in tangential orbits, the radial velocity may change very little during pericentric passages. Third, we expect particles to approach positive radial velocity on average as we approach the turnaround radius, which indicates that they are part of the Hubble flow rather than that they are orbiting the halo. To screen out unphysical pericentres, we demand that the first (but not any subsequent) pericentric passage should coincide with a switch in the angular position of the particle. For radial orbits, this switch occurs suddenly and in conjunction with the sign switch in v r . For tangential orbits, the velocity switch may be gradual, but the angular direction slowly approaches anti-alignment with the original infall direction (Fig. 4) . For the purposes of validating potential first pericentres, we demand that φ < 0.5 at the snapshot before or after v r changes sign. This criterion is robust to small shifts in the radial and angular trajectories. For example, if v r becomes positive but φ is too close to unity, a pericentre will be triggered at a subsequent snapshot as long as φ keeps decreasing and v r stays positive (as is be the case for pericentres that are reasonably well resolved in time).

By visual inspection, we find that this basic algorithm reliably identifies the vast majority of first pericentres. There are, however, a few special cases. First, pericentres in nearly circular orbits may not be accompanied by a sign switch in v r due to noise in the trajectory (third column of Fig. 4) . We identify this situation if no pericentre has been detected yet but φ falls below −0.8, that is, if the particle is almost exactly anti-aligned with its infall direction. 2 We cannot, however, detect the exact nature or time of such a pericentre and thus set a 'lower limit' flag that indicates that the current pericentre count of zero may be incomplete. When constructing density profiles, we will add particles with this flag into the orbiting term regardless of their pericentre count. Second, we cannot find the first pericentre of particles that are already inside a halo when it is first detected by the halo finder. Again, we set the lower limit flag, meaning that the orbiting profile could erroneously include infalling particles for the first few snapshots during which a halo is alive. Conversely, particles that have orbited and are outside R 200m when the halo is created are counted into the infalling population. Both problems are irrelevant in practice because the halo finder detects haloes with as few as 20 particles while we analyse only haloes with at least 500 particles (Section 2.7), which will have existed for many dynamical times. A third challenging situation occurs when orbits are unresolved in time, i.e., if the orbital time is shorter than the snapshot spacing of the simulation. In this case, the conditions for a pericentre detection may never be fulfilled but the velocity tends to switch sign in a seemingly random fashion (right column in Fig. 4 ). We screen for such orbits by tracking the total number of times the velocity changes from negative to positive. We only count instances that are plausible pericentres, that is, if they occur within R 200m and at least half a dynamical time (one radius-crossing time) after infall into the halo. If we find three such changes in velocity but have not yet detected a pericentre, we conclude that the orbit is unresolved and mark it as a lower limit. This case is so rare that makes no appreciable difference in the resulting profiles.

In summary, we record a first pericentre if we find a switch to positive velocity accompanied by φ < 0.5 and a lower limit if φ < −0.8, if a particle is found inside a newly discovered halo, or if we find three switches of the radial velocity. We have tested this algorithm extensively (Appendix A). We do not attempt to compute the ex-act time or radius of each pericentre because we are interested only in the first snapshot at which a particle is counted into the orbiting term. Our algorithm works with only three snapshots in time per orbit, which optimizes the memory footprint of the method. We have ensured that the outcome does not depend on the maximum radius to which particle trajectories are followed by Sparta. If a particle strays beyond this limit, its trajectory is interrupted but its orbit count is stored and reactivated in case the particle re-enters the halo.

The resulting orbit counts differ from halo to halo and depend on the simulation's box size and cosmology. We can get a general sense by considering the roughly 17 million particles tracked in Test-Sim100 at z ≈ 0 (that is, all particles that ever came within 2 R 200m of a halo). Out of those, 15% have orbited N orb = 0 times, out of which 12% are marked as lower limits (1.9% of the total population). The majority of the lower limit particles are truly orbiting, meaning that the population of uncertain classifications is very small. The most populated bin is N orb = 1 with 38% of the particles. As expected, the particle population in each orbit decreases in number thereafter, falling to less than 2% at N orb = 9 (the highest number tracked). The lower limit fraction is below 1% in all bins where N orb > 1. However, we emphasize that the results for N orb > 1 are not important for this paper and suffer from algorithmic uncertainties such as unresolved orbits (Fig. 4) . It would be a challenging undertaking to obtain a numerically reliable, converged distribution of orbit counts at all N orb (see Sugiura et al. 2020 , for an attempt).

Given the split into orbiting and infalling particles, we now construct the corresponding cumulative mass profiles (also within Sparta). We choose 80 logarithmically spaced radial bins spanning the range [0.01, 10] R 200m,all (see Figs. 2 and 3 for two examples). Ideally, we would extract profiles to even greater radii, but the runtime and memory consumption of Sparta increase steeply as the particle distribution stored on each process becomes the dominant factor. For example, the largest halo in the L0063-WMAP7 box has R 200m = 1.8 h −1 Mpc at z = 0, leading to an outer profile radius of 18 h −1 Mpc or almost 0.3 times the box size.

One algorithmic subtlety is that we only find out about pericentres one snapshot after they have occurred (because we require two consecutive positive velocities, Section 2.4). Particles that have just undergone pericentre would be erroneously counted into the infalling term, an error that affects the profiles over a wide range of radii because particles can travel a significant fraction of R 200m between snapshots (D17). We correct this error by adding such particles to the orbiting term of the profile recorded in the previous snapshot. This correction, however, fails at the very last snapshot of the simulation because we cannot record pericentres that lie in the future. Thus, the last (usually z = 0) orbiting term is very slightly underestimated; we use the profile at the second-to-last snapshot (z = 0.03 in our simulations) instead.

There are two numerical effects that render profiles unreliable at small radii, namely two-body relaxation and reduced centripetal forces due to force softening (e.g., Moore et al. 1998; Klypin et al. 2001; Power et al. 2003) . The recent calibration of Ludlow et al. (2019) suggests that two-body relaxation should affect the profiles in the WMAP7 simulations by less than 5% outside a minimum radius of r conv = 0.086 l(z), where l(z) = L/N/(1 + z) is the inter-particle separation in physical length units (Table 1 ). The pre-factor changes to 0.091 l(z) and 0.133 l(z) for the Planck and self-similar simulations, respectively. Similarly, Mansfield & Avestruz (2021) suggest that force resolution should reduce the profiles by less than 5% at radii larger than four force softening lengths (4 , Table 1 ). In all plots and calculations, we exclude bins with a lower edge smaller than the more stringent of the two limiting radii. We present the detailed calculations that lead to our criteria, as well as quantitative tests, in Appendix A.

The density profiles of individual haloes are subject to significant random noise, as visualized in Figs. 2 and 3 . The resulting scatter (quantified in Section 3.4) means that it is often difficult to derive meaningful conclusions from individual profiles. Thus, we will mostly consider the averaged (mean and median) density profiles of halo samples selected by certain properties. The radial bins in r/R 200m typically contain different numbers of valid profiles because the radial resolution cut occur at different radii for different halos in a given sample, since they have different physical sizes and originate from different simulations. We compute the uncertainty on the mean and median by bootstrap subsampling the profile values in each bin 500 times. The resulting uncertainty is typically very small due to the large number of profiles and due to their good numerical convergence (Appendix A).

A difficulty arises when computing the averages in regions where the number of particles per bin approaches zero. The mean and median statistics exhibit different issues in this case: the mean gradually becomes noisier, whereas the median density may suddenly fall to zero. To avoid such spurious results, we compute the expected number of particles in a radial bin based on the mean ρ/ρ m profile of the overall halo sample, from the bin volume in each individual halo, and from the particle mass of the given simulation. We drop the radial bin from the mean profile of the sample if it contains fewer than 1000 particles total. Technically, this criterion introduces a bias because we used the mean profile itself to estimate the number of particles, but the statistical fluctuation scales as 1/ √ N and is thus negligible for a 1000 particle limit. For the median profiles, we exclude individual halo profiles from a radial bin if their expected particle count falls below 10. This method is unbiased because it relies on the sample's mean profile rather than on the individual halo's profile. We have verified that our limits remove spurious medians and noisy bins without biasing the results.

In principle, we might want to rescale profiles by a radius other than R 200m,all , which would necessitate an interpolation to a new radial grid. In this paper, however, we restrict ourselves to profiles scaled by R 200m,all because the transition region is most universal in this radius, at least compared to other spherical overdensity definitions (DK14).

Besides omitting unreliable radial bins, we also exclude entire halos that are insufficiently resolved or otherwise biased. First, we consider only haloes with N 200m,all ≥ 500. This limit is slightly less strict than that of DK14, but we find good convergence (Appendix A) in agreement with Mansfield & Avestruz (2021) , who showed that profile-related quantities such as V max are converged at 500 particles in the Erebos simulations.

Second, we include only host (isolated) haloes because the density profiles of subhaloes are ill-defined due to the (often dominant) contribution from the host. However, the distinction is ambiguous given different definitions of the halo boundary (Diemer 2021) . We have experimented with the host-subhalo relations from the halo catalogues of D20. Using a small radius such as R 500c means including haloes that reside well within R 200m of larger haloes, and thus higher mean and median densities at large radii. Conversely, using a large radius such as R sp,90% lowers the average outer profiles significantly. Physically speaking, we are not interested in the contributions from other haloes as they are a manifestation of the halo-halo correlation function as much as of the density structure around isolated haloes. These contributions mostly affect the infall term of low-mass haloes because high-mass haloes are less likely to reside near another halo of comparable or larger size. We can thus gauge the effect of neighbouring haloes by comparing the mean or median profiles of samples with different peak height, and we find that excluding subhaloes out to larger radii brings the profiles into better agreement. On the other hand, using a radius as large as R sp,90% defines a significant fraction of all haloes as subhaloes (Diemer 2021) .

As a compromise, we use R 200m,bnd as our fiducial definition of subhaloes but introduce another criterion to screen for haloes with a large neighbour contribution: the fraction of particles that are not gravitationally bound. In particular, we exclude haloes for which M 200m,all /M 200m,bnd > 1.5, although we emphasize that there is no 'correct' value for this cut (Appendix A). Our choice represents a compromise that ensures relatively minor neighbouring-halo effects and excludes a small fraction of haloes: 2% of all haloes in the WMAP7 sample at z = 0, up to 5% at higher redshifts, and similar fractions in the self-similar simulations. This cut mostly affects haloes with high accretion rates.

In summary, we include only haloes with N 200m ≥ 500 particles and M 200m,all /M 200m,bnd ≤ 1.5. Out of their profiles, we include only bins with radii r > max(4 , r conv ) from the halo centre. In all figures, we omit radial bins that contain fewer than 80 individual halo profiles, a compromise between showing profiles to small radii and hiding noisy bins that distract the eye.

We now combine the profiles of haloes from different simulations into one sample per cosmology and redshift. These samples can then be further split by mass, mass accretion rate, and so on. For the ΛCDM simulations, we have computed profiles at z = [8, 7, 6, 5, 4, 3, 2, 1.5, 1, 0.8, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.03]. Some of the highest redshifts have been omitted for the largest box sizes, where no haloes exist yet at early times (Table 1) . In Appendix A2, we check that these samples are converged across simulations with different mass and force resolutions. The combined WMAP7 halo sample contains about 510 000 haloes at z = 0, decreasing to 38 000 at z = 6. The Planck sample contains 320 000 haloes at z = 0.

For the self-similar simulations, we combine different redshifts into one sample per simulation because cosmic time is not physically meaningful. We compute profiles at redshifts that correspond to even spacings in dynamical time (Diemer 2020b) , namely, for n = −1 at z = [15, 11.5, 8.8, 6.8, 5.1, 3.8, 3.3, 2.8, 2.3, 2.085], for n = −1.5 at z = [12.4, 9.6, 7.3, 5.6, 4.2, 3.1, 2.2, 1.5, 1.01], for n = −2 at z = [11.7, 9.0, 6.9, 5.2, 3.9, 2.9, 2.1, 1.4, 0.9, 0.531], and for n = −2.5 at z = [4.3, 3.2, 2.3, 1.6, 1.0, 0.6, 0.27, 0.03]. The profiles exhibit the expected self-similarity to an accuracy of about 10% (Appendix A2). The combined halo samples vary in size from 1.9 million for n = −1 to 180 000 for n = −2.5.

As an example, Fig. 5 shows the averaged halo profiles of the WMAP7 sample at z = 0. The remaining profile figures follow the same scheme, where the three larger panels show the profiles of all, orbiting, and infalling matter (ρ tot , ρ orb , and ρ inf , from top to bottom). Below those panels, we show the logarithmic slope, d ln ρ/d ln r, which is computed using a 4th-order Savitzky & Golay (1964) fil-ter with a window size of 15 bins. This filter smooths over noisy bins without changing the physical features of the profiles (DK14). The Savitzky-Golay filter does not work for the 7 left and rightmost bins (half the window size), where we estimate the derivative with a Gaussian smoothing as described in Appendix B of Churazov et al. (2010, using a window width of ∆ R = 0.2).

In this section, we discuss the shapes of the orbiting, infalling, and total profiles. We focus on averaged profiles because those of individual haloes are subject to significant random noise. We consider haloes across vast ranges of mass, redshift, and cosmology, and combine them into samples binned by mass, redshift, and accretion rate. Given the large parameter space, we summarize the most important trends with a few representative examples and refer the curious reader to a more comprehensive selection of online figures (Section 1). In particular, we discuss the profiles of halo samples selected by only mass and redshift in Section 3.1 and Fig. 6 . We additionally bin haloes by accretion rate in Section 3.2 and Fig. 7 . We study the impact of cosmology in Section 3.3 and Fig. 8 . And finally, we consider halo-to-halo scatter in Section 3.4 and Fig. 9 .

We begin by binning haloes in ΛCDM cosmologies in redshift and mass, which we express as peak height to obtain a redshiftindependent measure of the relative size of haloes (Section 2.2). As demonstrated in Fig. 5 , we visualize our results as the total, orbiting, and infalling profiles from top to bottom, with smaller panels showing their logarithmic slope. We follow this pattern in Fig. 6 , which summarizes our results for mass-selected profiles in the WMAP7 cosmology. All conclusions of this section equally apply to the Planck cosmology, whose profiles are visually indistinguishable. We focus on the WMAP7 simulations because they cover a larger range of box sizes and halo masses.

In the left column of Fig. 6 , we compare the median profiles of halo samples at z = 0 from small galaxy haloes (ν = 0.5 or M 200m ≈ 2 × 10 10 M ) to massive clusters (ν > 3 or M 200m > 1.1 × 10 15 M ). At higher redshift, a fixed peak height corresponds to lower mass. At first sight, the profiles look similar across the large mass range, but this impression is partly caused by the large dynamical range. The slope panels demonstrate that the shape of the profiles does depend on mass. While the total profiles of low-mass haloes smoothly steepen out to about R 200m , those of massive haloes are slightly shallower at small radii but then steepen significantly near R 200m (as shown by DK14 and subsequent works). At larger radii, the profiles approach the shallower slope of the infalling profile.

With our dynamical split, we can now identify the pattern that leads to the profile shape in the transition region. The orbiting term experiences a sharp truncation at a radius that varies with mass (or rather, average mass accretion rate as we show in Section 3.2). The logarithmic slope approaches arbitrarily steep values at this radius. Given the rapidly decreasing density of orbiting matter, the infalling term quickly comes to dominate, leading to a thin transition zone where the slope sharply drops and recovers. The infalling term itself depends on peak height, with higher mass haloes having steeper infalling profiles (again largely an effect of the average accretion rate). Moreover, the infalling profiles seem to approach different asymptotic density values at small radii, though the profiles tend to become . Example of dynamically split profiles for two ΛCDM halo samples. Each panel shows the median profiles (lines) and 68% scatter (shaded areas) for low-mass haloes (0.5 < ν < 1, blue) and high-mass haloes (ν > 3, orange) in the WMAP7 cosmology at z = 0. The larger three panels show the total, orbiting, and infalling profiles; the smaller bottom panels show the respective logarithmic slopes, d ln ρ/d ln r. The remaining figures in this paper are arranged similarly. The halo-to-halo scatter in the slopes is not shown because it is difficult to compute the slope of individual profiles in a robust manner. The example shown highlights a number of the most important trends. The scatter is larger around low-mass haloes because they are more influenced by nearby neighbours (Appendix A4). The median profiles of lowmass haloes steepen gradually, while those of high-mass haloes exhibit a clear edge at the distribution of orbiting matter. The position and sharpness of this feature depend on mass via the different mass accretion rates (DK14). While the total profiles reach slopes of about −4, the slope of the orbiting profile becomes arbitrarily steep at the truncation radius. The infalling profiles also vary systematically with halo mass and exhibit significant halo-to-halo scatter, especially for low-mass haloes. Mean, z = 0 0.5 < ν < 1 1 < ν < 1.5 1.5 < ν < 2 2 < ν < 3 3 < ν < 6 10 −1 10 0 10 1

Median, z = 4 10 −1 10 0 10 1

10 −1 10 0 10 1 r/R 200m Figure 6 . Average density profiles of haloes in the WMAP7 cosmology, selected by peak height and redshift. As in Fig. 5 and in all following figures, the large panels show the total, orbiting, and infalling profiles (from top to bottom). In the left column, we compare median profiles at z = 0 across masses (from ν = 0.5 or M 200m > 2 × 10 10 in blue to ν > 3 or M 200m > 1.1 × 10 15 M in yellow). The truncation of the orbiting term (the splashback radius) occurs at smaller radii in higher-mass samples because they also have higher mass accretion rates on average. The infalling profiles also evolve with mass, mostly for a similar reason (Section 3.2). The second column shows the mean profile of the same samples. In low-mass halos, the mean orbiting term reaches larger radii because of a few haloes where matter has been removed to large radii by interactions with other haloes. The third column shows the median profiles at z = 4, which appear fairly similar to their z = 0 counterparts. We confirm this impression in the fourth column, where we track the evolution of the median profiles for 2 < ν < 3 haloes across redshift (the trends are representative for other mass bins as well). The redshift evolution of the orbiting term is modest, but the normalization of the infalling term increases, leading to an overall shallower slope in the total profiles. unreliable (and are thus not plotted) before they clearly converge to a fixed density.

The radius where the slope is steepest is often taken to be a proxy for the splashback radius, but the profiles in Fig. 6 highlight why this is an imperfect definition: as both the orbiting and infalling terms depend on mass, the radius where the slope is steepest is only approximately equal to the radius where the orbiting term truncates. The difference is borne out when comparing the radius of steepest slope to dynamically calculated splashback radii (D20, see also O'Neil et al. 2022 ). The truncation radius could be taken as a new proxy for R sp , but we note that it is defined by all orbiting particles, whereas the splashback radius is defined by the apocentres of only the most recently accreted particles (Bertschinger 1985; More et al. 2015) . We will explore the detailed connection between the truncation and splashback radii in Paper III.

Our conclusions have been based on median profiles, but observations generally constrain mean profiles second column of Fig. 6 ). We expect differences between mean and median for asymmetric distributions, where outlier profiles can strongly influence the mean. Indeed, we notice significant differences in the mean orbiting and infalling profiles of low-mass haloes. The orbiting term continues past the truncation radius that is apparent in the median profiles because a few disrupted haloes have significant 'orbiting' material at large radii, which disproportionately contributes to the mean value at radii where the orbiting overdensity is small (say, below ρ m ). Similarly, the mean infalling profile of the lowest mass sample exhibits a sharp plateau feature near R 200m . We suspect that this feature is related to the definition of the halo boundary because excluding haloes with large fractions of unbound material reduces it substantially (Appendix A4). Moreover, the feature occurs only in low-mass haloes that are more strongly influenced by their environment than their high-mass counterparts. Luckily, most observations of density profiles currently target group and cluster scale haloes with ν > ∼ 2, where the mean and median profiles agree much better.

In the third column of Fig. 6 we investigate the redshift dependence of the profiles by comparing to haloes of the same peak heights but at z = 4. The lowest ν bins are not accessible at this redshift because they correspond to haloes below the resolution limit of our simulations. In general, the profiles appear similar, including the trends with peak height. For a more direct comparison between redshifts, we plot the 2 < ν < 3 profiles for a large range of redshifts in the right column (the evolution of the other ν bins proceeds similarly). Overall, the profiles become slightly shallower with redshift at all but the largest radii. However, the orbiting term remains remarkably unchanged with redshift (see also Le Brun et al. 2018) , meaning that the evolution of the total profiles is driven by an increase in the normalization of the infalling profiles. Since all profiles are expressed in units of ρ/ρ m , this increase in normalization goes beyond the cosmological increase of ρ m (z) with redshift and turns out to be related to the power spectrum (Section 3.3).

Most of the trends with halo mass are actually consequences of the dynamical state of haloes, namely, of their mass accretion rate Γ (equation 1). Thus, profiles are most universal across mass, redshift, and cosmology when they are binned by Γ. Fig. 7 shows the median profiles for the same ν bins as in Fig. 6 but additionally split into four bins in Γ. The halos in the lowest bin barely accrete (or even lose) matter because the pseudo-evolution due to the evolving overdensity threshold of 200ρ m (z) corresponds to roughly Γ ≈ 0.5 (Diemer et al. 2013a ). The profiles for haloes with Γ > 4 are very similar to those in the 2 < Γ < 4 bin, in agreement with Diemer et al. (2017) , who found that R sp /R 200m approaches a constant value at high Γ. We only show results for WMAP7 and z = 0, but the trends are the same for the Planck cosmology and for higher redshifts. Fig. 7 largely confirms the central conclusions of DK14, which is not surprising given that they were based on the same simulations. Most notably, the truncation radii of the orbiting term (or splashback radii) are smaller for higher Γ. This trend arises because a particle with a given kinetic energy at infall will reach apocentre at a smaller radius if mass is added during its orbit. This effect is independent of halo mass (More et al. 2015) , as evidenced by the alignment of the truncation radii for different ν. One exception is the lowest ν bin (blue), where neighbouring haloes strongly influence the profiles. Conversely, the truncation appears sharpest in the most massive halos because they dominate their environment. All of these trends makes sense in the context of spherically symmetric infall models, where the boundary of the orbiting term is delineated by a single, well-defined shell of particles. We review such models in Section 4.2.

The inner profile (r < R vir ) also changes significantly with Γ. In slowly accreting halos, the profile slope decreases smoothly with radius, whereas the profiles of fast-accreting halos approach a power law with a slope of about −2 across a wide range of radii. Such a profile is a natural prediction of fast collapse with purely radial orbits (Section 4.3). In slowly accreting halos (Γ < 1, left two columns), we observe wiggles in the slope at r < R 200m , which reflect a second caustic caused by the second orbit of recently accreted particles (Adhikari et al. 2014 ). In high-Γ halos, this signature is washed out by particles on their first orbit. Based on our new algorithm, we can for the first time discern the effects of mass and accretion rate on the infalling profiles. At low Γ, they are close to power laws (fixed slope) with a shallower slopes in smaller halos. When there is little physical accretion, environmental effects matter: the profiles of small halos contain significant contributions from close neighbours (see also Fig. A2 and Appendix A4). At high Γ, the profiles become virtually independent of ν because they are driven by freshly accreted matter. Those profiles exhibit a characteristic shape where they steepen as matter approaches the truncation radius, although this trend mostly disappears when considering overdensity rather than density (Paper II). The profiles flatten again as particles reaches the centre of the halo and appear to asymptote to a limiting central density at small radii, which is higher for higher ν haloes. To our knowledge, these features of the infalling profiles have not been predicted by theoretical models, which we further discuss in Section 4.4.

Having considered the impact of mass and accretion rate, we now turn to the effects of cosmology. Our unit system aims to be as inherently cosmology-independent as possible: peak height takes into account the varying linear growth factors in different cosmologies, we plot ρ/ρ m (z) to take advantage of a natural scaling with the mean matter density, and Γ is a logarithmic accretion rate that scales out the zeroth-order cosmology dependence of halo mass. By 'effect of cosmology', we thus mean differences due to cosmological parameters at fixed peak height and redshift.

Comparing our WMAP7 results from Sections 3.1 and 3.2 to the Planck cosmology, we find the orbiting terms to be statistically indistinguishable. The normalization of the infalling term appears to be a few percent lower in the Planck cosmology, regardless of halo selection. We can exclude a linear scaling with Ω m because the change in that parameter is 0.32/0.27 − 1 ≈ 18%, much larger than the differences we find. However, the change in σ 8 and other parameters might partially offset the effects of the matter density. We also note that the matter-matter correlation function, ξ mm (r), is higher (rather than lower) in the Planck cosmology at the relevant radii. Without a larger range of simulated cosmologies, we cannot present any firm conclusions regarding the dependencies on Ω m and σ 8 .

We can, however, test the origin of differences due to cosmology more generally by considering Einstein-de Sitter universes. As structure formation is perfectly self-similar in these simulations, the profiles at fixed peak height agree across cosmic time (Appendix A2), leaving only one free parameter: the slope of the power-law initial power spectrum, n. While self-similar universes are not realistic, the dependence of halo properties on n is a proxy for ΛCDM cosmologies, where n varies with spatial scale (and thus with mass and redshift). This correspondence has been demonstrated for concentration Ludlow et al. 2016) , the abundance of subhaloes (Diemer 2021 0.5 < Γ < 1 0.5 < ν < 1 1 < ν < 1.5 1.5 < ν < 2 2 < ν < 3 3 < ν < 6 10 −1 10 0 10 1 r/R 200m 1 < Γ < 2 10 −1 10 0 10 1 r/R 200m 2 < Γ < 4 10 −1 10 0 10 1 r/R 200m Figure 7 . Median profiles of samples split by both peak height (colours) and mass accretion rate (columns) at z = 0 in the WMAP7 cosmology. All trends are similar for other redshifts and cosmologies. Comparing to the profiles selected by only ν in the left column of Fig. 6 , a number of salient differences are readily apparent. The accretion rate is a stronger determinant of profile shape than peak height, given that the profiles in each column vary only mildly with ν but significantly across columns. The splashback feature becomes more pronounced with mass accretion rate and moves towards smaller radii. The inner profiles of slowly accreting halos exhibit smoothly decreasing slopes, whereas fast-accreting profiles resemble power laws. The wiggles in slope near the truncation radius in slowly accreting haloes (left two columns) are due to the second orbit of particles (Adhikari et al. 2014) . The shape of the infalling term is also strongly influenced by Γ, with a more variable slope in fast-accreting haloes.

The left column of Fig. 8 shows the median profiles in four selfsimilar universes with slopes between n = −1 and −2.5. The general shape of the profiles is very similar to ΛCDM. The main effect of a shallower n is to cause slightly steeper orbiting profiles out to the truncation radius and vice versa, a trend compatible with the literature (Crone et al. 1994; Reed et al. 2005; Knollmann et al. 2008; . The slope of the infalling term shows little correlation with n, but a shallower n means a lower normalization. These trends are reminiscent of the redshift trends we identified in Section 3.1. To establish this connection more quantitatively, we compute an effective slope in ΛCDM. The simplest definition is d ln(P)/d ln(k), but Diemer & Joyce (2019) showed that a tighter connection to concentration (and presumably to halo density profiles in general) is obtained by using

where R L is the Lagrangian radius of a halo of a given peak height and redshift (Section 2.2) and κ is a free parameter of order unity. Setting κ = 1 and ν = 1.75 to conform with the peak height bin shown in the left column of Fig. 8 , we obtain n eff = −1.78 at z = 0 and n eff = −2.4 at z = 4. In the second column of Fig. 8 , we thus compare the WMAP7 results at z = 0 and z = 4 to the self-similar simulations with n = −1.5 and n = −2.5 (the former provides a closer match than n = −2). The correspondence is striking, both In each column, we compare the median profiles for haloes with 1.5 < ν < 2 for different cosmologies. In the right two columns, we additionally select by mass accretion rate and only include haloes with 2 < Γ < 4. In the first and third columns, we compare profiles from the four self-similar simulations with different power spectrum slopes n. In the second and fourth columns, we repeat the lines for n = −1.5 (orange) and n = −2.5 (magenta) but contrast them with the WMAP7 ΛCDM cosmology at z = 0 (red dashed) and z = 4 (blue dashed). These comparisons are exemplary of a number of general trends. The slope of the self-similar profiles is systematically shallower for steeper slopes of the power spectrum. This trend holds even when selecting by Γ and affects both the orbiting and infalling terms. The ΛCDM profiles at z = 0 are strikingly similar to n = −1.5 and those at z = 4 to n = −2.5, which demonstrates that much of the redshift evolution in ΛCDM is due to the different part of the power spectrum probed by different halo masses.

in the orbiting and infalling profiles and both in their normalization and slope. In the third and fourth columns of Fig. 8 , we additionally select haloes with 2 < Γ < 4 to check that the trends with n are not simply due to different mass accretion rates. Our previous conclusions largely hold for this sample, as well as for other Γ bins. However, we observe that the truncation radius in ΛCDM at z = 0 is smaller than that in the corresponding self-similar cosmology. A similar conclusion was reached by Diemer et al. (2017) , who saw large splashback radii in self-similar universes with shallow n.

Overall, the agreement between ΛCDM and self-similar cosmologies demonstrates that the redshift evolution of profiles in ΛCDM is better understood as a reflection of the changing power spectrum probed by haloes of different masses. Similarly, differences in the ΛCDM parameters translate into different power spectra, whose slope affects the density profiles of haloes. The exact correspondence need not be simple though because the profiles can be affected by other haloes across a wide range of masses, e.g., due to satellite accretion. Definitions such as that of n eff in equation (2) collapse a complex dependence on the entire power spectrum into a single parameter. . Logarithmic scatter in dex around median profiles for massselected samples at z = 0 in the WMAP7 cosmology (left) and for the n = −2 self-similar cosmology (right). As in previous figures, the three rows show the scatter around the total, orbiting, and infalling profiles (from top to bottom). Plots for other redshifts and cosmologies appear very similar. The scatter in the total profiles is a strong function of radius: at r < ∼ R 200m , it is more or less constant and takes on values roughly between 0.15 and 0.35 dex, regardless of mass, redshift, or cosmology. Around the truncation of the orbital term, the scatter increases dramatically due to different truncation radii in different haloes and due to different shapes of the edge of the orbiting term. The scatter in the infalling term is strongly mass-dependent because lowmass haloes are much more affected by their surroundings (see also Fig. 5 ). Their scatter can reach 1 dex at r > ∼ R 200m , whereas the scatter around the most massive haloes remains below 0.5 dex (although it increases towards the centre of haloes).

We began the section by stating that individual halo profiles are widely scattered around their median, but how widely? Fig. 5 gives a visual example of the 68% scatter around ν-selected median profiles, showing a larger scatter in the infalling than orbiting profiles and in low-ν samples, which tend do be dominated by the environment rather than the haloes themselves. We confirm these impressions quantitatively in Fig. 9 , which shows the radial dependence of the 68% scatter in dex (that is, the logarithm of the ratio between the 84th to the 16th percentile). The scatter of the orbiting profiles is strikingly uniform across radii and peak heights, between 0.15 and 0.35 dex out to r ≈ 0.5 R 200m . Around the truncation radius, the scatter rapidly grows to arbitrary values as different haloes have slightly different truncation radii that lead to enormous fractional differences in density. Beyond the truncation radius, the scatter is dominated by that of the infalling profiles, which does depend on halo mass. While the highest ν bin experiences a scatter of about 0.2 dex at 10 R 200m , the lowest bin experiences about 1 dex. This picture is essentially independent of redshift and cosmology in ΛCDM. Perhaps more surprisingly, the scatter profiles are also very similar in the self-similar simulations (right column of Fig. 9 ). When selecting by mass accretion rate in addition to mass, the scatter decreases as expected, but the rate of decrease depends on ν and Γ. For ν > ∼ 0.5 and Γ < ∼ 2, the scatter falls to between 0.1 and 0.2 dex at small radii, whereas it remains similar to the purely mass-selected sample for the lowest ν bin. This trend suggests that low-mass haloes always suffer from larger scatter, even if selected by accretion rate. Moreover, the scatter gradually increases with Γ and reaches 0.2 to 0.4 dex for Γ > 4. Some of this increase is caused by massive satellites, since high accretion rates often include large individual accretion events. The scatter in the outer profiles is only weakly affected by a Γ selection.

Overall, we conclude that the scatter in the inner profiles (inside the truncation radius) is modest, between 0.1 and 0.4 dex depending on the selection, mass, and accretion rate. The corresponding ranges are shown in Fig. 5 but would barely be visible in Figs. 6, 7, and 8. On the other hand, the scatter in the infalling profiles is substantial, up to one dex for low-mass haloes (see also Avila-Reese et al. 1999) .

We have dynamically disentangled the orbiting and infalling components of halo density profiles and systematically investigated their shape. We have confirmed that the mass accretion rate of halos is the primary determinant of their profiles, although peak height, redshift, and cosmology also play a role (the latter two via the slope of the power spectrum). In this section, we discuss possible reasons for the physical origin of the profile shapes.

The two-scale nature of the orbiting term On a particle level, the orbiting term represents the superposition of the phase space occupied by its constituent orbits (e.g., Lithwick & Dalal 2011; Sugiura et al. 2020) , much like for stars in a galaxy (Schwarzschild 1979) . Assuming that the distribution of total particle energies is limited by some maximum, there must be a corresponding radial cut-off roughly at the scale where the halo's potential reaches that maximum (as visualized in Figs. 2 and 3) . Thus, the orbiting profiles exhibit two fundamental spatial scales. The first is some inner scale, whose existence has been observed ever since the early days of N-body simulations and which can be expressed as the scale radius r s . The second scale is the truncation radius, which we shall denote r t following DK14; it depends on the accretion rate (or, technically, its rescaled form r t /R 200m does). While the growth rate is correlated with concentration via the accretion histories of haloes (e.g. Wechsler et al. 2002) , significant scatter in the c-Γ and Γ-r t /R 200m relations means that haloes cannot generally be described by a single scale. For example, the truncation radii in Fig. 7 appear aligned due to the selection by Γ, but the profiles exhibit slightly different scale radii for different peak heights.

The two-scale nature of the orbiting term has a profound consequence: density profiles cannot be 'universal' in the sense that they are fully described by single-scale functions such as the NFW or Einasto profiles. While simulated profiles more or less follow those forms at r < ∼ R 200c , different dynamical states lead to different outer shapes. Conversely, lining up the profiles by their truncation radius cannot remove the dependence of the inner profile on concentration. This discussion begs the question of whether it makes sense to align the mean and median profiles at R 200m , rather than at r s , r t , or perhaps an entirely different radius (e.g., García et al. 2021) . We would expect an alignment at r t to strongly reduce the scatter in the truncation region but to increase the scatter at smaller radii. We leave a detailed investigation to future work. In Papers II and III, we quantify r s , r t , and their relation based on a new fitting function.

In practice, one caveat to the two-scale nature is that not all haloes experience a sharp cut-off because some particles truly escape after orbiting, which can spatially extend the orbiting term to arbitrary distances. These extended profiles are much more prominent in lowmass haloes (Fig. 6) and are mitigated by excluding haloes with a large unbound component (Appendix A4), indicating that the majority of escaping orbits are caused by close encounters with haloes of comparable size. Given that the vast majority of haloes do experience a sharp truncation, we leave a more detailed investigation of such unusual orbits for future work.

An intuitive theoretical explanation for the sharp truncation is provided by secondary infall models, where shells of dark matter are accreted onto an initial power-law overdensity (Fillmore & Goldreich 1984; Bertschinger 1985; Ryden & Gunn 1987; White & Zaritsky 1992) . These models predict an infinitely sharp caustic at the location of the last outermost orbiting shell (the splashback radius). While they naturally explain the general trend of this location with accretion rate (Adhikari et al. 2014; Shi 2016) , the agreement with simulations is only approximate Sugiura et al. 2020) . Similarly, the smoother truncation in realistic profiles is readily explained by non-sphericity (e.g., Lithwick & Dalal 2011; Vogelsberger et al. 2011; Adhikari et al. 2014 ), but there are currently no quantitative predictions as to what the realistic shape of the truncation should look like. The secondary infall model does, however, naturally predict the phase-space density profile of halos (Ludlow et al. 2010) , and it explains the second caustic in the slope in low-Γ halos as particles on their second orbit (Fig. 7 , Adhikari et al. 2014) . The infinitely sharp shell in the model is transformed into a smoother, more realistic profile by effects such as non-sphericity (e.g., fig. 3 in Ryden 1993) .

Alternatively, we might try to understand the shape of the truncation from fundamental physics. For example, similar cut-offs are seen in profile models that are based on assumptions about the distribution of particle energies (e.g., King 1966; Hjorth & Williams 2010 ), but we show in Paper II that the resulting models are not a good description of our orbiting profiles because energy-based cutoffs lead to infinitely steep profiles (Amorisco 2021). One reason for this disagreement might be that the halo particles do not really follow a single, simple distribution of energies (e.g., Pontzen & Governato 2013) . Instead, their energies are determined by the collapse of a non-isotropic peak and by subsequently accreted subhaloes with their own energy distribution. Nonetheless, energy-based arguments can augment our understanding of mechanisms that cause a truncation of the phase-space distribution and thus of density profiles, for example in tidally stripped haloes (Drakos et al. 2017; Amorisco 2021) .

In summary, the sharp truncation is not a surprise given the results of shell models and energy-based calculations, but a detailed model of its shape will likely need to rely on a more detailed modelling of the distribution of particle orbits.

The secondary infall model also predicts the full density profile inside and outside the splashback radius, but it is well known that these predictions are not realistic, at least not unless the model is coupled with collapse from cored peaks (e.g., Dalal et al. 2010; Ogiya & Hahn 2018; Delos et al. 2019) , realistic accretion histories (Lu et al. 2006) , or mergers (Ogiya et al. 2016; Angulo et al. 2017) . One reason is that, by itself, the model can predict only power-law density profiles because it evolves a self-similar initial density perturbation. 3 While simulated profiles generally exhibit evolving slopes, approximate power-law shapes have been found under certain conditions (e.g., Tasitsiomi et al. 2004; Anderhalden & Diemand 2013; Delos et al. 2018) . We also find strikingly constant slopes in the regime of very rapidly collapsing haloes (Section 3.2). In this section, we thus check whether the secondary infall model predicts the correct trends with accretion rate and power spectrum slope.

We assume an initial overdensity with δ = ρ(< r)/ρ m ∝ r −3ε , where ε ≡ −d ln δ/d ln M(< r). Fillmore & Goldreich (1984) showed that the resulting density profile is ρ ∝ r −γ with a slope

The perturbation cannot be steeper than δ ∝ r −3 or ε = 1 because that corresponds to a point perturbation (the case treated in Bertschinger 1985) , which gives γ = 9/4. The lower limit of γ > 2 can be understood as a consequence of the purely radial orbits assumed in the model: at radii much smaller than their apocentre, r apo , particles have approximately constant velocity and thus spend equal amounts of time in shells with volume proportional to r 2 , leading to ρ ∝ r −2 (Dalal et al. 2010) . Thus, the secondary infall model produces power-law profiles with a relatively narrow range of slopes, 2 < γ < 2.25. Comparing to the right columns of Figs. 7 and 8, we indeed find that the slopes of fast-accreting halos are close to γ ≈ 2. The shallower inner profiles in simulations demonstrate the presence of tangential (or at least non-radial) orbits (e.g., Ryden 1993; Subramanian et al. 2000; Lu et al. 2006 ). Regardless of the exact final slope, the secondary infall model predicts that the final profiles strongly depend on the initial conditions. The shape of the initial peaks is largely determined by the slope of the power spectrum (Bardeen et al. 1986) , and this shape imprints itself on the mass accretion histories and thus on the final profiles (Hoffman & Shaham 1985; Cole & Lacey 1996; Del Popolo et al. 2000; Manrique et al. 2003; Ricotti 2003; Ascasibar et al. 2004; Lu et al. 2006; Salvador-Solé et al. 2007; Dalal et al. 2010; Ludlow et al. 2013; Polisensky & Ricotti 2015) . Mathematically, we can establish the connection between n and γ by assuming that the structure of a density peak is well described by the two-point correlation function, δ(r) ∼ ξ(r) (e.g., Hoffman & Shaham 1985; Ogiya & Hahn 2018) . In self-similar cosmologies, ξ ∝ r −3−n and thus

which would give γ = 2/1.5/1 for n = −1/ − 2/ − 2.5. This prediction is clearly invalid in the regime of purely radial orbits and depends on the model assumptions (e.g., Subramanian et al. 2000) .

However, the general trend that higher n lead to steeper profiles (e.g., Hoffman & Shaham 1985; Łokas 2000) is borne out in our comparison of self-similar cosmologies (Fig. 8) , although the range of slopes is much smaller than suggested by equation (4). The power spectrum slope has also been shown to affect the curvature of the orbiting term, where a steeper n means a more rapidly steepening profile (Salvador-Solé et al. 2012; Udrescu et al. 2019; Brown et al. 2022 ). We will quantify this effect parametrically in Paper III. We can also establish a connection between the predicted profile slope and the accretion rate. Typical structures collapse when the density inside their Lagrangian radius crosses the threshold δ c ≈ 1.686. By the definition of Γ, the Lagrangian mass is M L ∝ a Γ ∝ R 3 L and thus R L ∝ a Γ/3 . The self-similarity of the model implies that Γ is constant, meaning that it does not matter that we actually evaluate it over a dynamical time rather than instantaneously. Combining the previous expression with the linear evolution of overdensities in an Ω m = 1 universe, δ ∝ a, and our definition of ε,

or, equivalently, d ln M/d ln r ∼ 3Γ/(3 + Γ) (Adhikari et al. 2014; Shi 2016) . For halos with Γ = 0/1.5/3/6 we predict γ = 3/2/1.5/1, or generally a shallower profile for higher accretion rates. This trend is borne out in Fig. 7 : all profiles reach γ ≈ 1.5 at the smallest radii we can test, but Γ ≈ 0 halos reach steepest slopes of γ = 3 at the truncation radius while Γ > 2 halos barely steepen past γ = 2. This small variation means that the profiles are, on average, fairly close to a power law, which reminds us of the argument for the asymptotic slope of −2: if a halo is accreting rapidly, many particles are on recently established orbits with apocentres near r t , meaning that most of the density profile falls into the r r apo regime. Moreover, the orbits of freshly accreted particles have had less time to become isotropized, e.g., via the radial orbit instability (e.g., Merritt & Aguilar 1985; MacMillan et al. 2006) .

In summary, the secondary infall model can give insights into aspects of the formation of the orbiting profile and, with some additional assumptions, predicts the correct trends with n and Γ as well as power-law profiles with γ = 2 for fast-accreting halos. However, the detailed shape of the profile is a function of the entire mass accretion and merger history (Ludlow et al. 2013 ) and thus carries the complexities of non-linear structure formation. A theoretical model that fully describes the resulting profiles has yet to be formulated.

The shape of the infalling term inside the transition region has hitherto been a mystery because its density is dominated by the orbiting term by many orders of magnitude. Our dynamical split has allowed us to systematically investigate the infalling term for the first time. Its exact shape at small radii could be fairly sensitive to the pericentre counting algorithm because the particle number per logarithmic radial bin declines towards the centre, leading to the omission of the infalling lines at radii smaller than 0.01 to 0.1 R 200m in Figs. 5-8. However, given our corrections for time binning effects (Section 2.5) and the results of convergence tests (Appendix A), we have no reason to suspect a systematic bias in the infalling profiles.

Overall, we find that the average infalling profiles are highly varied, with clear dependencies on mass, redshift, accretion rate, and cosmology (via the power spectrum). However, almost all infalling profiles steepen with radius and reach their steepest slope around the truncation radius before flattening at larger radii. The latter as-pect depends on the plotted quantity though: the infalling profiles in overdensity, ρ inf /ρ m − 1, are fairly well fit by a single power law (DK14 and Paper II). We have chosen to plot ρ/ρ m because the overdensity can become negative. Regardless, it is clear that the infalling profiles do not share a common slope.

What shapes might we expect for the infalling profiles? Perhaps the simplest expectation would be that infalling shells of dark matter contract gravitationally but do not cross, akin to applying the Zel'dovich (1970) collapse model to spherical peaks (e.g., Mohammed & Seljak 2014) . This assumption leads to ρ ∝ r −3/2 (Bertschinger 1985) , which is reflected in the overdensity slopes scattering roughly around −1.5 (DK14). We expect this picture to break down around the splashback (or truncation) radius, where shells cross for the first time. Indeed, we find that the profiles become shallower around this radius because infalling shells 'feel' an interior mass distribution that decreases as they approach the halo centre.

Another simple, cosmology-dependent model for the outer profiles is to multiply a mass-dependent bias with the matter-matter correlation function, ρ inf /ρ m = b(ν)ξ mm (r)+1 (e.g., Hayashi & White 2008; Tinker et al. 2010 Tinker et al. , 2012 . While the predicted profiles are relatively close to power laws, their slopes do not quite match simulations (DK14). Moreover, Fig. 6 makes it clear that there are significant slope variations with peak height, which cannot be captured by a scale-independent bias parameter. The predicted profile does evolve with redshift due to the evolving comoving scale probed by halos at fixed ν and due to the different time scalings of ρ m (z) and ξ mm (z), but the evolution does not match that of the simulated profiles. Nonetheless, the influence of cosmology on the outer profiles is strikingly borne out in our results, since the redshift evolution of the infalling profiles at fixed ν can be entirely explained by the effective slope of the power spectrum (Section 3.3).

In summary, our results demonstrate that the infalling profiles are far more complex than simple power laws. Not only do they flatten at small radii, their slopes at large radii vary as a function of numerous halo properties and cosmology. We will further explore these dependencies in Papers II and III.

We have systematically investigated the density profiles of dark matter haloes in N-body simulations between 0.01 and 10 R 200m , probing wide ranges of halo mass, redshift, and cosmology. For the first time, we have split the profiles into orbiting and infalling particles based on a novel algorithm that counts each particle's pericentres. Our main conclusions are as follows.

(i) The orbiting profiles experience a sharp truncation at the edge of the orbital apocentre distribution. This splashback radius is generally apparent in the total profiles but not always located exactly at the radius where the total slope is steepest.

(ii) The truncation and scale radii exhibit different dependencies on halo properties, meaning that the orbiting profiles cannot be described as a function of a single scale parameter as suggested by many common fitting functions.

(iii) The mass accretion rate is the main factor determining the shape of density profiles. Slowly accreting haloes exhibit gradually steepening orbiting profiles whereas fast-accreting profiles tend towards power-law-like behaviour with mildly varying slopes. Peak height is only a secondary determinant of the orbiting profiles but has a significant impact on the infalling profiles.

(iv) The slope of the infalling term strongly varies from halo to halo, as well as with peak height and mass accretion rate, demonstrating that there is no uniform outer halo profile.

(v) The profile evolution with redshift in ΛCDM can be understood as a variation in the effective slope of the power spectrum probed by haloes of different physical sizes, highlighting that density profiles are powerful diagnostics of cosmology. We can think of ΛCDM as part of a continuum of self-similar universes with different power spectra.

(vi) The logarithmic 68% scatter in the orbiting profiles is virtually constant inside the truncation radius and varies between 0.1 and 0.4 dex depending on the sample selection, peak height, and accretion rate. Lower mass haloes exhibit systematically larger scatter, especially in the infalling term where the scatter is mostly driven by the environment.

(vii) We demonstrate the convergence of our dynamical profile splitting algorithm with mass, force, and time resolution. However, we find that the density profiles are sensitive to definitions such as the halo boundary and the exclusion of merging systems.

By splitting density profiles into their orbiting and infalling components, we can now systematically understand the profile shape in the transition region. However, numerous fundamental aspects remain to be explored, for example, an accurate theoretical model of the orbiting term based on particle dynamics. We have made our code and data publicly available in the hope of stimulating further investigations.

I am deeply indebted to my PhD adviser, Andrey Kravtsov, for first suggesting halo density profiles as a research topic and for overriding my objection that everything worth knowing about them had already been explored. I thank Shaun Brown for feedback on a draft and the anonymous referee for their constructive suggestions. I am also grateful for enlightening discussions with Sten Delos (on theoretical collapse models) and with Phil Mansfield (on numerical convergence issues). Furthermore, I thank Susmita Adhikari, Han Aung, Rafael Garcia, Oliver Hahn, Daisuke Nagai, and Eduardo Rozo for productive conversations. This work was partially completed during the Coronavirus lockdown and would not have been possible without the essential workers who did not enjoy the privilege of working from the safety of their homes. 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. This research extensively used the python packages Numpy (Harris et al. 2020) , Scipy (Virtanen et al. 2020) , Matplotlib (Hunter 2007) , and Colossus (Diemer 2018) .

The Sparta framework is publicly available in a BitBucket repository, bitbucket.org/bdiemer/sparta. An extensive online documentation can be found at bdiemer.bitbucket.io/sparta. The Sparta output files (one file per simulation) are available in an hdf5 format at erebos.astro.umd.edu/erebos/sparta. A Python module to read these files is included in the Sparta code. Additional figures are provided online at benediktdiemer.com/data. The full particle data for the Erebos N-body simulations are too large to be permanently hosted online, but they are available upon request.

In this appendix, we derive the resolution limits applied to our halo selection (A1) and present tests that demonstrate the numerical convergence of our simulations and profile splitting algorithm (A2-A4). The figures shown are merely examples of a much larger range of data that we have inspected.

Any N-body simulation can produce reliable halo profiles only over a certain range of halo masses and spatial scales, and it is difficult to establish formal convergence criteria because different numerical effects can interact with each other (Moore et al. 1998; Klypin et al. 2001; Power et al. 2003) . Possible effects include spurious interactions between individual particles (known as two-body relaxation, Binney & Tremaine 2008), underestimated forces (and thus profiles) within a few force softening lengths of the centre of the halo, discreteness effects (Splinter et al. 1998; Joyce et al. 2009 ), and integration errors due to the finite time-steps. The latter two are thought to be subdominant, especially since the Gadget code uses adaptive timestepping that essentially ensures a certain number of time-steps per orbit at radii greater than the force resolution (Ludlow et al. 2019; Mansfield & Avestruz 2021) .

The remaining question is whether reduced centripetal forces or two-body relaxation are the more important, and how to limit their impact. As relaxation depends mostly on the particle number enclosed within a certain radius, the relative strength of the effects depends on /l, the ratio of the force softening and inter-particle separation scales (l ≡ L/N). Our strategy is to compute the minimum radii where profiles should be converged with respect to either effect, and to use the more stringent of the two criteria.

Specifically, we impose a target accuracy of σ res = 5% in the circular velocity profiles. The effect of force resolution was numerically investigated by Mansfield & Avestruz (2021) and encapsulated in the fitting function σ res (r) = exp[− (Ah/r) β ] with h = /0.357, A = 0.172, and β = −0.522 for Gadget simulations. We invert this function and evaluate it for σ res = 0.05,

where we have rounded to r = 4 for convenience. To derive a comparable limit for the radius where we expect two-body relaxation effects to fall below 5%, we rely on the Power et al. (2003) parametrization of κ P03 ≡ t relax /t orb , the ratio of the relaxation to orbiting time-scales. Analogous to equation (A1), Ludlow et al. (2019) provide a mapping between κ P03 and the accuracy of convergence (their equation 21), which we again invert to find

with (a, b, c, d) = (−0.4, −0.6, −0.55, 0.95). For σ res = 0.05, this formula gives κ P03 = 0.37, meaning that the orbiting timescale should roughly be three times the relaxation time-scale in order to avoid resolution effects of more than 5%. The radius where this criterion is satisfied, r conv , can be computed via an implicit relation that depends on the density profiles (and thus concentrations) of halos as well as on the simulation parameters (Power et al. 2003) . We approximate this calculation using equation (18) of Ludlow et al. (2019) , which is a function of only the inter-particle spacing,

where C = 2.44 and l(z) = l/(1 + z) is the evolving inter-particle spacing in physical units. We have verified that this approximation provides a good description of the full, implicit calculation for our simulation parameters and the halo masses in question. Inserting our target accuracy of κ P03 = 0.37, we find r conv = 0.133 Ω 1/3 m l(z), or r conv = 0.086 l(z) for the WMAP7 simulations, r conv = 0.091 l(z) for the Planck simulations, and r conv = 0.133 l(z) for the self-similar simulations.

For a given simulation, the minimum radii for both two-body relaxation and force resolution now correspond to a fixed comoving scale, meaning that we can directly compare their strength. For the WMAP7 simulations, two-body relaxation dominates over force resolution when r conv > 4 or /l < 1/47 (1/44 in the Planck cosmology). Thus, two-body relaxation is the limiting factor in the 125 h −1 Mpc and smaller boxes, and force resolution is the limiting factor in the 250 h −1 Mpc and larger boxes. In the self-similar simulations, two-body relaxation dominates if /l < 1/30; given the very small softenings of /l < 1/195 in the n = −1 and n = −1.5 simulations, the results are entirely limited by two-body relaxation in those boxes (r conv ≈ 7 × 4 ).

We now numerically verify the criteria derived in the previous section by directly comparing haloes in different simulation at fixed peak height. Given that the simulations systematically vary in mass and force resolution, non-convergence would manifest as disagreements in profiles that should be indistinguishable. Fig. A1 shows representative examples of our convergence tests, which include both the regimes where errors are dominated by force resolution (large boxes in ΛCDM) and by two-body relaxation (small boxes and self-similar simulations). In the left three columns, we consider the WMAP7 cosmology at z = 0 and z = 2. The coloured lines show the mean (left column) and median (second and third columns) profiles from individual simulations as coloured lines. The black lines indicate the means and medians of the entire sample compiled from all simulations, to which we compare the individual simulations in the bottom panels. The shaded areas show our bootstrap uncertainties (Section 2.6). We restrict the profiles to our fiducial limits, namely, only halos with N min ≥ 500 and only radial bins with Self-similar, n = −1.5 1.9 < ν < 2.0 Median z = 1.0 z = 1.5 z = 2.2 z = 3.0 z = 4.2 z = 5.6 n = −1.5 10 −1 10 0 10 1 r/R 200m Figure A1 . Representative examples of convergence tests. Each set of panels shows density profiles (top) and the fractional difference (bottom) to the profile of the combined sample (black lines). The coloured lines refer to different simulation boxes of a ΛCDM cosmology (left three columns) and to different redshifts in the n = −1.5 self-similar simulation (right column). Each set of boxes and redshifts corresponds to different resolutions for haloes at fixed peak height, and their agreement thus demonstrates convergence. As in the other figures, we do not show noisy profile bins with fewer than 80 individual profiles.

r min ≥ max(4 , r conv ). Generally speaking, the total profiles are converged to within 5% once the uncertainties are taken into account. This finding holds for both the mean and median profiles and at all redshifts.

The convergence of the total density profiles reflects our N-body simulations, but it does not imply that our algorithm to split particles into orbiting and infalling is also insensitive to resolution effects. We consider this more challenging test in the second and third rows of Fig. A1 , which show that the split profiles are converged to roughly 10% or better. The only notable exception occurs in the mean orbiting profiles at large radii, which fall to extremely low densities (e.g., 10 −4 ρ m in the left column of Fig. A1 ). Here, individual outlier profiles can drastically affect the mean, for example, if a halo is tidally ripped apart. We should thus be cautious when interpreting the mean orbiting profiles at large radii.

For the self-similar simulations, convergence can be tested by comparing different redshifts of a single simulation. Given that both the power spectrum and time are independent of physical scales, the profiles at fixed peak height must match on average, except for resolution effects Colombi et al. 1996; Jain & Bertschinger 1998; Widrow et al. 2009; Joyce et al. 2021; Leroy et al. 2021; Maleubre et al. 2022) . The right column of Fig. A1 shows profiles in a representative peak height bin for different redshifts in the n = −1.5 simulation. Once again, we find excellent convergence, even in the orbiting and infalling profiles.

So far, we have compared halo samples that were selected with our fiducial resolution limits, N min = 500 and r min ≥ max(4 , r conv ). We visually justify these choices in the left two columns of Fig. A2 , where we directly compare our fiducial sample to others with stricter and less strict cuts. We find no significant differences between the Mass/force resolution Self-sim., n = − 10 −1 10 0 10 1 r/R 200m Figure A2 . Same as in Fig. A1 , but instead of different simulation boxes we now compare samples defined by different resolution cuts (left two columns), time resolutions (third column), and neighbour criteria (right column). In the left two columns, we observe no statistically significant differences beyond the expected 5% convergence between our fiducial sample (N ≥ 500, r > 4 , κ P03 = 0.37) and a version with stricter criteria (purple; N ≥ 2000, r > 8 , κ P03 = 0.74), which indicates that our cuts in mass and force resolution are sufficient. This statement holds for both ΛCDM (left column) and self-similar simulations (second column). The median profiles are even better converged than the mean. Samples with less stringent resolution criteria (r > 1 and κ P03 = 0.1 in red, and N ≥ 100 in yellow) lead to visible deviations. The third column compares profiles from our test simulation with three different time resolutions (snapshot spacings). We detect no significant differences between the high-resolution version (dark blue) and that with the usual 100 snapshots (light blue). In the right column, we show the impact of our criterion to exclude contributions from nearby neighbours on low-mass haloes. As we decrease the allowed M 200m,all /M 200m,bnd ratio, the outer profile keeps decreasing and approaches the profile shape found for more massive haloes. Our chosen limit of 1.5 represents a compromise.

fiducial and stricter samples, which leads us to conclude that our limits are sufficiently strict. Furthermore, we have tested that the scatter around the mean and median values is not increased in our sample compared to higher resolution samples. On the other hand, less stringent cuts do have an impact on the mean. For example, including haloes with as few as 100 particles leads to moderate, although statistically significant, deviations. Similarly, including those parts of profiles that suffer from poor mass and force resolution can lead to the characteristic down-turn of the profiles where they are underestimated due to a reduced forces and two-body scattering (left column of Fig. A2 and Appendix A1).

Given that our profile splitting algorithm relies on particle orbits, it is subject to yet another potential resolution issue, namely the finite number of snapshots per orbit. For example, pericentres can be missed for orbits that are unresolved in time, which could lead to particles being erroneously classified as infalling (Fig. 4) . We test the convergence with time resolution in the third column of Fig. A2 . This test can only be performed with our smaller test simulation, where we can subsample 200 snapshots. We find no significant differences between the full version and that subsampled to about 100 snapshots at any mass or redshift. The slight differences in the total profiles (top row of Fig. A2 ) are due to the slightly different halo

Galactic Dynamics: Second Edition

Galaxy Formation and Evolution

These results give us confidence that the 100 snapshots of the other Erebos simulations are sufficient (see D17 for similar conclusions)

The pericentre detection algorithm presented in Section 2.4 has a number of free parameters that could change the resulting split profiles. One such parameter is the radius at which particles are first tracked, which represents the largest radius at which pericentres can occur. This radius is set to 2 R 200m (D20). We find that increasing this radius does lead to a small number of additional pericentres and thus to a small shift of matter into the orbiting term. However, by visual inspection, we find that most of the additional pericentres appear to be spurious and due to irregular particle motions at large radii that cannot be clearly associated with a pericentre. Nevertheless, the uncertainty in how far out we allow pericentres slightly influences the shape of the orbiting profile at large radii, where the overdensity falls well below unity. Given that there simply is no 'right answer' for what constitutes a pericentre, we conclude that these shapes inherently depend on our definitions to some extent.Other parameters of the algorithm include the minimum angular distance a particle must have traversed at its first pericentre (φ < 0.5) and the angular distance at which a pericentre is triggered in tangential orbits without clear radial velocity switches (φ = −0.8). These values were determined by visual inspection of numerous orbits, but varying them within a reasonable range (φ > ∼ 0 as a minimum and φ ≈ −1 as a certain pericentre) does not noticeably change the profiles.

In Section 2.7, we introduced a limit on the ratio of total and bound mass, M 200m,all /M 200m,bnd < 1.5, to weed out haloes that are strongly affected by close neighbours. The impact of this cut on the mean profiles of low-mass haloes is shown in the right column of Fig. A2 . Here, we consider haloes with high accretion rates where the trend is most visible, which is somewhat by construction as those haloes will typically have grown recently because they are 'accreting' material from another halo. Cutting out such haloes reduces the density in the outer regions, which mostly affects the infalling profile but also has an effect on the orbiting term. The figure highlights that our choice of M 200m,all /M 200m,bnd < 1.5 is a compromise between maintaining a complete halo sample and keeping the profiles free from the effects of nearby neighbours. Importantly, median profiles are much less affected than mean profiles because the effect is driven by relatively few extreme cases. Moreover, the convergence between simulation boxes (Fig. A1) is not affected by the neighbour criterion.Nearby haloes also appear to be at least partly responsible for the break in the slope that occurs around R 200m in some infalling profiles.No such feature appears when haloes with nearby neighbours (or rather, even small fractions of unbound mass) are strictly excluded (yellow lines in Fig. A1 ), but it becomes gradually more prominent as the cut is relaxed. Similarly, the orbiting profiles of haloes with unbound material reach farther out, indicating that particles can be tidally removed by interactions. However, the exact change in the profile shapes due to neighbours is complex, and we defer a more thorough investigation to future work.