Accurate modelling of the onetotwo halo transition has long been difficult to achieve. We demonstrate that physically motivated halo definitions that respect the bimodal phasespace distribution of dark matter particles near halos resolves this difficulty. Specifically the two phasespace components are overlapping and correspond to 1 particles it orbiting rm the halo and 2 particles it infalling rm into the halo for the first time. Motivated by this decomposition Garca R. Garca et. al. MNRAS 521 2464 2023 advocated for defining haloes as the collection of particles orbiting their selfgenerated potential. This definition identifies the traditional onehalo term of the halomass correlation function with the distribution of orbiting particles around a halo while the twohalo term governs the distribution of infalling particles. We use dark matter simulations to demonstrate that the distribution of orbiting particles is finite and can be characterised by a single physical scale rrm h which we refer to as the it halo radius. rm The twohalo term is described using a simple yet accurate empirical model based on the Zeldovich correlation function. We further demonstrate that the halo radius imprints itself on the distribution of infalling particles at small scales. Our final model for the halomass correlation function is accurate at the approx 2 level for r in 0.150 h1 Mpc. The Fourier transform of our best fit model describes the halomass power spectrum with comparable accuracy for kin 0.06 6.0 h Mpc1.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Dynamicsbased halo model for large scale structure
{'Dynamics-based halo model for large scale structure': "Edgar M. Salazar, 1, ∗ Eduardo Rozo, 1 Rafael Garc'ıa, 1 Nickolas Kokron, 2, 3 Susmita Adhikari, 4 Benedikt Diemer, 5 and Calvin Osinga 5  \n1 Department of Physics, University of Arizona, Tucson, AZ 85721, USA  \n2 Department of Astrophysical Sciences, Princeton University, Princeton, NJ, 08540, USA  \n3 School of Natural Sciences, Institute for Advanced Study, Princeton, NJ, 08540, USA  \n4 Department of Physics, Indian Institute of Science Education and Research, Homi Bhaba Road, Pashan, Pune 411008, India  \n5 Department of Astronomy, University of Maryland, College Park, MD 20742, USA  \n(Dated: September 11, 2024)  \nAccurate modelling of the one-to-two halo transition has long been difficult to achieve. We demonstrate that physically motivated halo definitions that respect the bimodal phase-space distribution of dark matter particles near halos resolves this difficulty. Specifically, the two phase-space components are overlapping and correspond to: 1) particles orbiting the halo; and 2) particles infalling into the halo for the first time. Motivated by this decomposition, Garc'ıa [R. Garc'ıa et. al., MNRAS 521, 2464 (2023)] advocated for defining haloes as the collection of particles orbiting their self-generated potential. This definition identifies the traditional one-halo term of the halo-mass correlation function with the distribution of orbiting particles around a halo, while the two-halo term governs the distribution of infalling particles. We use dark matter simulations to demonstrate that the distribution of orbiting particles is finite and can be characterised by a single physical scale r h , which we refer to as the halo radius. The two-halo term is described using a simple yet accurate empirical model based on the Zel'dovich correlation function. We further demonstrate that the halo radius imprints itself on the distribution of infalling particles at small scales. Our final model for the halo-mass correlation function is accurate at the ≈ 2% level for r ∈ [0 . 1 , 50] h -1 Mpc . The Fourier transform of our best fit model describes the halo-mass power spectrum with comparable accuracy for k ∈ [0 . 06 , 6 . 0] h Mpc -1 .", 'I. INTRODUCTION': "The halo model is a well-established tool for describing the N-point functions characterising the large scale structure of the Universe [1, 2]. In it, two point clustering statistics are expressed as a sum of two terms: 1) a one-halo term that describes pairs of points that fall within a single halo; and 2) a two-halo term that describes pairs of points that fall within different haloes. The conceptual simplicity of the model has made it an indispensable tool for describing clustering statistics. However, the model is famously inaccurate in the translinear regime, exhibiting up to 20% biases [3]. Consequently, its utility as a tool of precision cosmology is limited [2].  \nIn the standard halo model, the one halo term is typically described using an NFW [4] or Einasto [5] profile, whereas the two halo term is described as the product of a halo bias times the linear matter correlation function ξ lin ( r ) . To improve the fit of the model in the translinear regime, van den Bosch et al. [6] developed accurate fiting functions for a scale-dependent  \nbias. While accurate, these fitting functions provide little insight into the physics that lead to said scale dependencies.  \nFortunately, our understanding of halo structure has progressed significantly over the past decade. The first step towards our current understanding was made by Diemer and Kravtsov [7], who demonstrated the existence of a 'splashback' feature in the density profile of dark matter haloes. This feature is associated with the apocentric radius of particles in their first halo orbit, and can be used as a well understood and physically-motivated halo boundary [8-10]. Interestingly, if one uses the splashback radius to define the halo boundary, the resulting halo mass function is more universal than that derived using traditional halo definitions [11].  \nAlthough physically motivated, the fact that the apocenters of splashback particles are distributed over a range of radii implies that the definition of the splashback boundary is not unique [12]. Garc'ıa et al. [13] demonstrated this feature is imprinted in the halo-mass correlation function, which can be used to uniquely define a radius-mass relation for dark matter haloes. They further found that doing so enabled them to construct a percent-level accurate - though somewhat ad hoc description of the halo-mass correlation function.  \nIndependently, Fong and Han [14] proposed an alternative  \nhalo boundary definition called the depletion radius. It is nominally defined in terms of the maximum of the mass-flow rate I m ( r ) = 4 πr 2 ρ ( r ) v r ( r ) for the inner depletion radius [see also 15]. However, the infall velocity profile v r ( r ) is never negative at low masses. The authors note that at high mass the halo depletion radius coincides with the minimum of the ratio ξ hm ( r ) /ξ mm ( r ) , and adopt this as their definition of the inner depletion radius for low mass haloes. The inner depletion radius roughly coincides with the splashback boundary, and was used in Zhou and Han [16] to construct an accurate model of the halo bias b ( r, M ) ≡ ξ hm ( r | M ) /ξ mm ( r ) on small and intermediate scales ( r ≲ 20 h -1 Mpc ).  \nBuilding on the insights from these works, the last two years have seen a new paradigm for describing the structure of dark matter haloes begin to emerge. [17] proposed that particles near a halo can be separated into orbiting or infalling particles according to whether or not they had experienced their first pericentric passage. In a follow-up paper, [18] used this decomposition to arrive at fitting functions that can accurately ( ≈ 5% ) describe the density profile of dark matter haloes.  \nThe proposal by [17] was significantly strengthened by [19]. It has long been known that the phase space distribution of particles around dark matter haloes is inherently bimodal [20]. [19] demonstrated that: 1) the two components of this phase space distribution could be cleanly separated using the accretion history of the particles; and 2) that the resulting populations of dark matter correspond to the orbiting/infall split proposed by Diemer. They further proposed that dark matter haloes ought to be defined as the collection of all particles orbiting their self-generated potential. We refer to haloes defined in this way as dynamical haloes . 1 In this proposal there is no halo boundary; the orbiting and infalling populations of particles overlap in phase space. Moreover: 1) the halo mass function of dynamical haloes is Press-Schechter with a slowly moving barrier; and 2) the halo bias of dynamical haloes can be accurately predicted using the peak-background split formalism [21].  \nThe proposal by Garc'ıa et al. [19] has important consequences for the halo model. Specifically, because dynamical haloes are defined as the collection of all orbiting particles, it follows that the one-halo term describes the distribution of orbiting particles in a halo, while the two-halo term describes the distribution of infalling particles. Moreover, because the orbiting and infalling particle distributions overlap in phase space, the one- and two-halo terms necessarily overlap in configuration space.  \nIn this work, we construct an updated halo model based on the definition of dynamical haloes proposed in Garc'ıa et al.  \n[19], focusing specifically on the halo-mass correlation function. We construct an accurate fitting formula for the orbiting and infall components, and arrive at a model of the halo-mass correlation function that is accurate to ≈ 2% (at 1 σ ) on scales r ∈ [0 . 1 , 50] h -1 Mpc . As detailed below, our results provide new insights into the structure of dark matter haloes and their environments.  \nThis paper is organised as follows. Section II presents the simulation we use and how the halo catalogue is constructed. In sections III and IV we characterise the orbiting and infall profiles of dark matter haloes respectively. Section V presents our complete model for the halo-mass correlation function and characterises its accuracy. We summarise our work and discuss some of the advantages and disadvantages of our model relative to alternatives in the literature in Section VI.", 'II. SIMULATION DATA AND HALO CATALOGUE': "This work relies on the Cold Dark Matter (CDM) simulation described in Banerjee et al. [22]. The simulation was run using GADGET-2 [23], and contains 1024 3 particles in a 1 h -1 Gpc box. The softening length is 15 h -1 kpc . The cosmological parameters are Ω m = 0 . 3 , Ω Λ = 0 . 70 , n s = 0 . 96 , h = 0 . 7 , and A s = 2 × 10 -9 , corresponding to σ 8 = 0 . 85 . The initial conditions are set at z = 99 using N-GenIC [24]. In this paper, we focus exclusively on describing the halomass correlation function at z = 0 . An initial halo catalogue is obtained using the ROCKSTAR halo finding algorithm [25]. For this initial catalogue we adopt a standard spherical overdensity definition with an overdensity threshold of 200 with respect to mean. Starting from the most massive ROCKSTAR halo in the catalogue, we classify every dark matter particle within 5 h -1 Mpc of each halo as either orbiting or infalling. The halo is then redefined as the collection of all orbiting particles, and these particles are not allowed to be members of any other halo. We then move on to the next largest halo, and iterate as needed. The end result is a halo catalogue where every halo is comprised exclusively of orbiting particles. We refer to these haloes as dynamical haloes. Additional details can be found in Garc'ıa et al. [19].  \nFigure 1 illustrates how dynamical haloes differ from traditional haloes [see also Fig. 3 in 17]. The figure shows the mass distribution and phase space density of the largest halo in the simulation, with an orbiting mass M orb = 2 . 35 × 10 15 h -1 M ⊙ . The top row shows the distribution of dark matter particles around the halo, and the bottom row shows the corresponding phase space distribution. The left column shows all the particles within a 6 h -1 Mpc box around the halo; the middle column, the orbiting particles only; and the right column, the infalling particles. We see that: 1) the phase space distribution of the halo has two distinct but overlapping components; and 2) the traditional R 200 halo boundary does  \nFIG. 1. Dark matter particle distribution in the vicinity of the most massive halo in the simulation (left), separated into orbiting (middle) and infalling (right) populations. Top: Configuration space y -x projection. The circular dashed line denotes a standard halo radius definition R 200 m for this halo. Bottom: Phase space v r -r projection. The vertical dashed line corresponds to R 200 m .  \n<!-- image -->  \nnot properly separate these two components. By contrast, the orbiting/infall classification proposed by Garc'ıa et al. [19] adequately captures the dynamical structure of the halo.  \nWe measure the halo-mass correlation function ξ hm ( r ) of dynamical haloes binned by orbiting mass. When mass is measured in units of h -1 M ⊙ , the bin edges are log 10 M ∈ [ 13.40, 13.55, 13.70, 13.85, 14.00, 14.15, 14.30, 14.45, 14.65, 15.00 ] . The minimum halo mass is set by the requirement that haloes contain at least 300 orbiting particles. For reference, this is roughly equivalent to M 200m ≈ 2 × 10 13 h -1 M ⊙ haloes. Figure 2 shows ξ hm ( r ) as measured in the simulation. We also show the orbiting (dashed line) and infalling components (dotted line) computed using the orbiting/infall split of Garc'ıa et al. [19]. Both components add up to ξ hm ( r ) . This decomposition is clearly evocative of the halo model, and is the primary motivation for this work. In what follows, we describe both the orbiting and infall contributions to the halo-mass correlation function in terms of simple and, where possible, physically motivated fitting functions.", 'III. THE ORBITING PROFILE': "The left panel in Figure 3 shows the mean orbiting matter density profile for dynamical haloes in bins of M orb . The covariance matrix for the data is obtained by jackknifing the box,  \nFIG. 2. The halo-mass correlation function of dynamical haloes as measured in our simulation. The grey region corresponds to scales smaller than six times the force softening, and are unresolved. The dashed and dotted lines are the orbiting and infalling profiles respectively as measured in the simulations based on the orbiting/infall split proposed by Garc'ıa et al. [19]. By construction, these two contributions exactly add to the full halo-mass correlation function. Errors come from jackknifing the box into 64 regions.  \n<!-- image -->  \n<!-- image -->  \nFIG. 3. Left: Orbiting density profile for stacked haloes in bins of orbiting mass. The errors bars are jackknife. The best fits for each profile are shown as solid lines. Right: Fractional error between data and model: (top-left) each mass bin is fit individually; (bottom-left) r h , α ∞ and a are power laws of mass; (top-right) r h and α ∞ are power laws of mass, and a is mass independent; (bottom-right) r h and a are power laws of mass, and α ∞ is mass independent. In all cases the fit is performed is restricted to scales r ≥ 90 h -1 kpc .  \n<!-- image -->  \nand the error bars correspond to the square-root of the diagonal elements of the covariance matrix. Note that this error does not account for any uncertainty associated with particle classification. Consequently, we expect our error bars to be somewhat underestimated. The grey region corresponds to a scale cut of 90 h -1 kpc , or six times the force resolution of the simulation [26]. We see that the profiles exhibit a slowly steepening slope at small radii, and a sharp exponential decay at large radii. We can achieve an accurate fit to the data with the functional form  \nρ orb ( r | M orb ) = A ( x a ) -α ( x ) exp [ -x 2 2 ] (1)  \nwhere we have defined the dimensionless radius x = r/r h , and r h is the characteristic length scale associated with the exponential truncation of the orbiting profile. The 'slope' α ( r ) of the profile varies slowly with r . We choose to parameterise this dependence so that α ( r ) extrapolates to constant values as r → 0 and r →∞ . We write  \nα ( x ) = α ∞ x a + x (2)  \nThis functional form has α smoothly transitioning from α = 0 to α = α ∞ at a scale r = ar h . The value α = 0 at r = 0 is arbitrary since these scales are not resolved in the simulation. We chose this convention since it keeps the profile finite at r = 0 while providing an accurate fit to the data over the scales we resolve.  \nThere are four parameters in our model, { A,r h , a, α ∞ } . However, the normalisation A is not independent: it is set by the condition that the integral of the orbiting profile ρ orb ( r | M orb ) be equal to M orb . Incorporating this constraint  \nrequires some care due to the impact of resolution effects. We describe the necessary corrections in Appendix A.  \nThe parameter r h governs the spatial extent of the halo, so we refer to it as the halo radius. Note that it is the only physical scale associated with this halo profile. In particular, while the running of the slope introduces a new scale ar h , we find that the parameter a can be modelled as mass independent, thereby establishing r h as the sole characteristic length scale in our profile. Finally, the parameter α ∞ is the asymptotic slope of the power-law component of the halo profile. We find α ∞ is close to isothermal (i.e. -2 ), though it slowly decreases with increasing mass. Consequently, our halo profiles are not self-similar: small haloes are not just scaled versions of larger haloes.  \nOur choice of parameterisation appears to violate the basic finding of Diemer [18], who argued the orbiting profile of haloes depends on two separate length scales, a scale radius r s , and a truncation radius r t . However, the underlying phenomenology is the same. The fact that r s /r t is not constant reflects the lack of self-similarity in the orbiting profile of haloes. Indeed, the scale radius r s is defined in terms of the logarithmic slope of the density profile, so a change in α ∞ in our model corresponds to a change in r s in the Diemer [18] model. The physical degrees of freedom are the same in both cases, but they are parameterized differently.  \nThe likelihood of the orbiting density profile for haloes in the k th mass bin is  \nln L k ( ρ orb ( r | M orb ) | θ ) ∝ -1 2 D ⊤ k ˜ C -1 ρ orb D k -1 2 det( ˜ C ρ orb ) (3)  \nwhere θ = { r h , a, α ∞ } is the vector of model parameters,  \nD k = ρ data orb -ρ model orb , and ˜ C ρ orb = C ρ orb +diag( δ 2 k ρ 2 orb ) (4)  \nThe last term in the above covariance matrix characterises the systematic uncertainty in the model. That is, δ k is the percent error we must add to account for the observed differences between our model and the data.  \nThe best fit model to the orbiting profile for each individual halo mass bin is shown as solid lines in the left panel of Figure 3. The residuals from the fit are shown in the top-left residuals sub-panel. The horizontal dotted lines in the residual panels correspond to a 2% difference between the data and the best fit model.  \nFigure 4 shows how each of our model parameters scales with mass. Unsurprisingly, r h increases with increasing mass: the more massive a halo, the larger it is. Note that while the shape parameter a controls the detailed shape of the halo profile at small scales, the scale r = ar h is unresolved in our simulations, being roughly three times the softening length. Consequently, we do not ascribe it physical significance; it is simply a means of parameterising the running of the power-law index of the halo profile over the scales we resolve.  \nTo construct a model for the density profile as a function of mass and radius, we parameterise the halo-mass dependence of each of our model parameters as a power law in mass. We then simultaneously fit all mass bins for the amplitudes and slopes of these power-law relations. We will refer to this parameterisation as the smooth model . The power laws for our model parameters are given by  \nr h ( M orb ) = r h ,p ( M orb M p ) r h ,s (5)  \nα ∞ ( M orb ) = α ∞ ,p ( M orb M p ) α ∞ ,s , (6)  \na ( M orb ) = a p ( M orb M p ) a s . (7)  \nThe likelihood of the smooth orbiting density profile is then  \nln L ( { ρ orb ( r | M orb ) } N k =1 | θ ) = ∑ k ln L k (8)  \nwhere θ = { r h ,p , r h ,s , α ∞ ,p , α ∞ ,s , a p , a s , δ } is the vector of model parameters. We use a single parameter δ to characterise any required excess variance across all mass bins. We set a flat prior ( -∞ , 0] for α ∞ ,s , and [0 , ∞ ) on the rest. Moreover, we assume no correlation between mass bins. We chose a pivot mass M p = 10 14 h -1 M ⊙ , which roughly decorrelates the amplitude and slope parameters in our fits.  \nThe best fit smooth model obtained when fitting all bins simultaneously is nearly indistinguishable from the individual best  \nfits, and is therefore not shown in Figure 3. The corresponding residuals are shown in the bottom-left sub-panel. The band about zero characterises the error in our fits, and includes the best fit value for the added statistical uncertainty parameter δ . Evidently, the smooth model remains accurate to better than 2% across all masses. Indeed, the best fit value for δ , which characterises the accuracy of our model, is δ = 1 . 0% . The smooth model predictions for each of the orbiting profile parameters r h , α ∞ , and a are shown as black solid lines in Figure 3. Note these lines are not fits to the data points in those panels, but rather fits to the complete set of stacked orbiting halo profiles.  \nTraditional halos are typically defined using an overdensity criterion, which leads to a radius-mass relation with a slope of 1/3. By contrast, we find that the slope of the r h -M orb relation is shallower than 1 / 3 at 3 σ : r h , s = 0 . 226 ± 0 . 003 .  \nTurning to the power laws for α ∞ and a as a function of mass are both very close to zero. For this reason, we also test the accuracy of our model for two other cases: 1) α ∞ is a power law in mass and a is mass independent; and 2) a is a power law in mass and α ∞ is mass independent. The residuals from the best fit model in each of these two cases are shown in the top-right and bottom-right sub-panels of Figure 3 respectively. The corresponding curves for these parameters as a function of mass are also shown in Figure 4.  \nFrom the right-most residual panels in Figure 3 it is clear that making α ∞ mass independent severely degrades the quality of the fits. Conversely, modelling the parameter a as mass independent hardly impacts the quality of our best fit model. Given that the parameter a corresponds to an unresolved scale, and that making this parameter mass-dependent results in negligible improvement in the quality of our fits, we have chosen to model a as mass-independent.  \nIn summary: the orbiting profile of dynamical haloes can be accurately ( ≈ 1% ) described as an exponentially truncated power-law (equation 1). The halo radius r h characterises the exponential truncation, and is the only dimension-full parameter in our model. As expected, r h increases with mass, but the slope of the relation is significantly shallower than the traditional value of 1 / 3 . The slope of the power-law part of the density profile varies slowly with radius as per equation 2. This introduces 2 new parameters: a , and α ∞ . Of these, α ∞ scales with mass as a power-law, whereas a can be modelled as mass independent. Thus, the total number of free parameters in our model is three: two parameters r h and α ∞ that depend on mass, and a population parameter a . The maximum likelihood values for the parameters in our fiducial (constant a ) smooth model are shown in Table I. The corresponding posterior distributions are shown as the red contours in Figure 10.  \nFigure 5 demonstrates that the infall correlation function is feature-rich. Traditionally, one sets the two halo term of ξ hm to ξ 2 -halo hm = bξ lin , where ξ lin is the linear correlation function. More generally, on large scales one can set  \n<!-- image -->  \n<!-- image -->  \nFIG. 4. Best-fit orbiting profile model parameters as a function of halo mass. The data points correspond to the individual profile fit maximum likelihood values, and the error bars are the 68% confidence interval. The solid lines are not a fit to these data points, but rather the best fit smooth model to the orbiting profile for three smooth parameterisations, as labelled.  \n<!-- image -->", 'IV. THE INFALL PROFILE': 'We wish to characterise the infall contribution to - or, equivalently, the two-halo term of - the halo-mass correlation function. Since our orbiting/infall split extends only to a radius of 5 h -1 Mpc away from a halo, we assume all particles at distances larger than this are infalling. This is already an excellent approximation for all haloes at r ≈ 5 h -1 Mpc . For scales less than 5 h -1 Mpc , the infall correlation function is defined as ρ inf ( r ) / ¯ ρ m -1 . On scales larger than this we set ξ inf = ξ hm . We measure the halo-mass correlation function up to a maximum radius of 50 h -1 Mpc . On larger scales, the measurements become noisy due to the cosmic variance in the simulation box. The left panel of Figure 5 shows the product r 1 . 5 ξ inf as a function of radius for each of our halo mass bins. The vertical axis is multiplied by r 1 . 5 to reduce the dynamical range of the y -axis and accentuate the non-power-law features in the correlation function. All errors are jackknife. Our goal is to fit this data.', 'A. The Large Scale Limit of ξ inf': "ξ inf = bξ large . (9)  \nfor an appropriate model of ξ large . Indeed, the ratio ξ hm /ξ mm is very nearly constant for r ∈ [10 , 50] h -1 Mpc [19]. Consequently, we seek a function ξ large that provides an accurate description of the matter correlation function on the scales resolved by this work.  \nWhen defining ξ large we could rely on existing approxima-  \nions for ξ mm , for instance HALOFIT [28-30] or emulators [31-34]. However, we wish to make our model as interpretable as possible. We had anticipated using Effective Field Theory (EFT) models for ξ large [e.g. 27, 35, 36, and others], but found these models ill-suited for our purposes. Specifically, figure 6 compares the z = 0 mass correlation measured using the 15,000 Planck simulations of the Quijote simulations suite [37] with the best-fit EFT model based on one-loop Lagrangian Perturbation Theory (LPT). We see that on scales r ≲ 20 h -1 Mpc , the LPT model: 1) deviates strongly from the simulation data; and 2) has obvious features not present in the simulation data. If we were to choose a one-loop LPT description of ξ large , our small scale parameterization for ξ inf would have to 'undo' these features.  \nFor these reasons, we rely instead on the empirical approximation ξ mm ≈ Bξ ZA where B = 0 . 948 and ξ ZA is the Zeldovich correlation function. As shown in Figure 6, this approximation is accurate to better than 1% over the scales r ∈ [20 , 50] h -1 Mpc . Moreover, the failure of this approximation at small scales is both featureless and less severe than that of the one-loop LPT model, making it a more useful model for our purposes. Consequently, we define  \nξ large ≡ Bξ ZA . (10)  \nwhere B = 0 . 948 is held fixed. The inequality B < 1 is due to the impact of non-linear growth on ξ mm , as evidenced by the fact that: 1) the prefactor B is a function of redshift, with B → 1 as redshift increases; and 2) on scales r ∈ [30 , 50] h -1 Mpc , the one-loop LPT predictions for ξ mm (with a counterterm) agrees with the approximation Bξ ZA with better than 1% accuracy. We caution that on very large scales one must recover the perturbative limit ξ mm → ξ ZA as r → ∞ , or B = 1 . Consequently, we anticipate that our model does not currently extrapolate to linear scales. We will address this problem in future work.  \n<!-- image -->  \nFIG. 5. Left: Infalling halo-mass correlation function for haloes in bins of orbiting mass. The errors bars are jackknife. The best fits for each profile are shown as solid lines. Right: Fractional error between data and model when fitting each bin individually (top), or when performing a joint fit to all mass bins using our smooth parameterisation (bottom). In all cases the fit is performed on scales above the resolution limit of the simulation ( r ≥ 90 h -1 kpc , shown in grey).  \n<!-- image -->", 'B. The Small Scale Limit of ξ inf': "On small ( r ≲ 10 h -1 Mpc ) scales ξ hm deviates from both ξ mm and ξ large in a scale dependent way. This leads us to define a scale dependent bias β ( r ) such that ξ inf = bβ ( r ) ξ large . The function β ( r ) only has an impact on small scales, and is therefore distinct from the higher order parameters appearing in perturbative expansions (e.g. b 2 , b s 2 , b 3 ). We find empirically that β ( r ) increases with decreasing radius, which suggests β should take the form 1 + ( r inf /r ) γ . To avoid divergences at r = 0 we adopt a cored power-law instead, and hypothesise that the core radius should scale with the halo radius r h . These assumptions lead us to write  \nβ ( r ) = 1 + [ r inf µr h + r ] γ . (11)  \nBecause r inf characterises the deviations from linear bias due non-linear growth, we anticipate r inf is related to the nonlinear scale k NL defined via ∆ 2 lin ( k NL ) = 1 . Indeed, we find r inf is independent of halo mass, with k NL r inf ≈ 2 / 3 .  \nOur model is not yet complete, however. Specifically, we find that the orbiting profile imprints itself onto the distribution of infalling matter at small scales. This is shown in Figure 7, where we plot the ratio ξ inf ( r | M ref ) /ξ inf ( r | M ) where M ref is our most massive bin. This ratio is then plotted as a function of r/r h , where r h is the halo radius of haloes of mass M . Evidently, the halo radius r h is imprinted onto the distribution of infalling particles, further cementing r h as the characteristic  \nradial scale for the halo.  \nTo fit the ratios in Figure 7 we define x ≡ r/r h , and then fit the ratio using the function  \nf ( x ) = C + ηxe -x (12)  \nThe fits are shown as solid lines in Figure 7. The vertical displacement parameter can be set to C = 1 by fitting to the ratio b ref ξ ref inf /bξ inf , where b is the large scale bias parameter discussed above. This means that we only need one additional parameter ( η ) to fully characterise this feature in the infall correlation function.  \nOur final model for the infall correlation function is  \nξ inf ( r ) = β ( r ) 1 + ηxe -x bξ large ( r ) . (13)  \nwhere β ( r ) is given by Eq. 11, and x ≡ r/r h . The parameter η characterises the amplitude of the impact of the orbiting particles on the infall profile.  \nOur final model for ξ inf has five parameters: { b, η, r inf , γ, µ } , each of which has a simple physical interpretation: b is the large scale halo bias; η characterises the amplitude of the 'dip' in the infall term due to the gravitational potential generated by the orbiting mass of the halo; r inf and γ characterise the power-law behaviour of the small scale dependent bias; and µ characterises the scale at which β ( r ) reaches its maximum value. Note that in this section we do not consider r h to be a model parameter that needs to be fit for. Rather, we adopt  \nFIG. 6. Top : The mass-mass correlation function measured for 15,000 realisations of Quijote's fiducial cosmology at z = 0 . The points with error bars correspond to the ensemble average, and the error bars show the cosmic variance in a single box, as opposed to the error on the mean. The lines correspond to the Zel'dovich correlation function, the correlation function Bξ ZA where B is a constant multiplicative offset chosen to match the non-linear correlation function at r ≈ 50 h -1 Mpc scales, and the 1-loop LPT prediction calculated using the code VELOCILEPTORS [27]. Bottom : Ratio of these various correlation functions to the mean correlation function measured in Quijote.  \n<!-- image -->  \nthe best fit value for r h obtained from fitting the orbiting profiles.", 'C. Infall Fits': "Our final model for ξ inf is set by Eq. (13) where ξ large is defined via Eq. (10). We fit the simulation data to our model for ξ inf using the likelihood in Eq. (3), where we simply exchange ρ orb for ξ inf . Because of the way our measurements were made, the covariance matrix of the extended infalling correlation function, C ξ inf , is a block-diagonal matrix, as we had no easy way to estimate the covariance between the  \nFIG. 7. Ratio of the infalling halo-mass correlation function of the largest mass bin to all other mass bins. A 'bump' appears at r = r h as we decrease in mass, becoming more prominent for the lowest masses. The solid lines are a fit to the data points using Eq. (12).  \n<!-- image -->  \nr ≤ 5 h -1 Mpc and r ≥ 5 h -1 Mpc for ξ inf . We set flat priors [0 , ∞ ) on all the parameters. Our best fit model is shown as the solid lines in the left panel of Figure 5. The top-right panel of the same figure shows the fractional residuals between the data and our model. We see that our model is accurate to 5% or better.  \nFigure 8 shows the maximum likelihood estimates for each of our model parameters as a function of M orb . The error bars represent the 68% confidence regions obtained from the parameter posteriors. Remarkably, both r inf and µ are massindependent. The fact that r inf is mass independent supports the notion that it is related to the non-linear scale k NL . Likewise, the fact that µ is mass independent emphasises the fact that r h is the unique radial scale associated with haloes of a given orbiting mass.  \nWe now turn to parameterise a smooth model for ξ inf as a function of mass. We parameterise the mass dependence of b ( M orb ) and γ ( M orb ) with power laws. 2 The mass dependence of η ( M orb ) is more complicated, and is modelled using an error function. 3 Both r inf and µ are modelled as being mass independent. Our updated model parameters for this smooth  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFIG. 8. Infalling profile model parameters. The data points correspond to the individual profile fit maximum likelihood values, and the error bars are the 68% confidence interval. The solid lines are not a fit to the data points shown above, but rather the best fit smooth model to the infalling profile.  \n<!-- image -->  \nmodel are given by  \nb = b p ( M orb M p ) b s , (14)  \nγ = γ p ( M orb M p ) γ s , (15)  \nη = η 0 2   1 -erf   log 10 [ M orb M p ] -η m η σ     (16)  \nWe fit for the smooth model parameters and again choose a pivot mass M p = 10 14 h -1 M ⊙ . The likelihood of the smooth infalling correlation function is given by Eq. (8) with the appropriate exchange ρ orb → ξ inf . Likewise, the vector of model parameters is now θ = { b p , b s , γ p , γ s , η 0 , η m , η σ , r inf , µ } , i.e. we describe our 9 ξ inf curves with one parameter per curve. We assume no correlation between mass bins when performing our simultaneous fit. The maximum likelihood parameter values for our smooth ξ inf model relations are shown in Table I. Our smooth mass model is nearly indistinguishable from the single-bin fits in the left panel of Figure 5, so we do not show them there. The percent residuals for our best fit smooth model are shown in the right-most panel of Figure 5. The corresponding curves for the single-bin parameters as a function of mass are shown as solid lines in Figure 8. The posterior distribution of our infall parameters is shown as blue contours in Figure 10. The accuracy of our best fit model is characterised using the excess variance parameter δ , for which we recover δ = 2 . 1% . The single most discrepant point in our model is ≈ 5% offset from our model fit.  \nIn summary: our infall model assumes ξ inf ∝ bξ mm , where ξ mm is described using the empirical approximation ξ mm ≈ Bξ ZA on scales r ∈ [10 , 50] h -1 Mpc . To accurately describe smaller scales we introduce a scale-dependent bias β ( r ) which we model as cored power-law. This introduces three new parameters: the amplitude and slope of the power-law ( r inf and γ ), and the core radius. The core radius is a constant multiple µ of the halo radius. Finally, we found the orbiting density profile imprints itself in the infall profile, which necessitates introducing one additional free parameter ( η ) that  \ncharacterises the amplitude of this effect. The total number of parameters required to describe ξ hm is three mass-dependent parameters ( b , γ , and η ), plus two population parameters µ and r inf that apply to all halos. Our model describes the infall profile on all scales with ≈ 2% accuracy.", 'A. Results': "Following the halo definition advocated by Garc'ıa et al. [19], one can identify the one- and two- halo terms of the halomass correlation function as those arising from the orbiting and infalling population respectively. Note that because all particles in a simulation are classified, this decomposition inherently occurs in density space. That is, we ought to write ρ tot = ρ orb + ρ inf . The correlation function takes the form ξ hm ( r ) = ρ tot / ¯ ρ m -1 , so there is some ambiguity as to whether the '-1' goes with the orbiting term, the infall term, or both. In configuration space, it is most natural to assign the '-1' to the infall term, as we have done so far. With this convention, our model for the halo-mass correlation function can be written as  \nξ hm ( r ) = ρ orb ( r ) ¯ ρ m + ξ inf ( r ) . (17)  \nThe left panel of Figure 9 shows the halo-mass correlation function for each of our mass bins as data points with error bars. All errors come from jackknifing the box. The solid lines correspond to the sum of the orbiting and infall terms predicted from out best bin-by-bin fits of the orbiting and infall profiles discussed in sections III and IV. The fractional errors from this model are shown in the bottom-left panel of the same figure, and are typically ≲ 2% , though the worst fit points exhibit residuals as large as ≈ 8% .  \nCritically, however, the model curves in the left panel of Figure 9 are not a fit to these data. These are the model curves de-  \nFIG. 9. Halo-mass correlation function in bins of orbiting mass. The errors bars are jackknife. Left: The solid lines are the sum of the best fit orbiting and infalling profiles, and are not a fit to the points in this figure. Right: Our best fit smooth model to the data points in the Figure. Bottom: Fractional error between data and model. In both cases the fit is performed above the shaded region corresponding to 90 h -1 kpc . The coloured band on the residuals plot on the left shows the jackknife error bars of the data, while the plot on right adds in quadrature the best-fit scatter parameter δ that characterises the accuracy of the fit.  \n<!-- image -->  \nfined by the fits to the orbiting and infalling components from sections III and IV, which in turn rely on the orbiting/infall split adopted in Garc'ıa et al. [19]. It is therefore possible that errors in particle classification compromise the accuracy of our model.  \nTo test for this possibility, we use our full model (orbiting+infalling) to fit the halo-mass correlation function directly. Note that while our functional form for ξ hm model was informed by the orbiting and infalling profiles, the fit is completely agnostic to how or even if the particles are classified. The resulting model fit is shown in the right hand panel of Figure 9, and the corresponding residuals are shown in the bottom-middle column. The accuracy of our model is characterised be the excess variance needed to fit the data, for which we find δ = 1 . 62% . The worst-fit data point differs from the model by ≈ 5% .", 'B. The Impact of Particle Classification': "Figure 10 compares the smooth model parameter posteriors obtained using two different analysis: 1) fitting ξ hm directly; and 2) fitting ρ orb and ξ inf independently as per sections III and IV. Evidently, the posterior distributions do not agree [see also 18]. These results suggest the presence of systematic errors in the orbiting/infall decomposition of Garc'ıa et al. [19].  \nTo further test this hypothesis we have evaluated the impact of small, plausible modifications to the Garc'ıa et al. [19] algorithm for particle classification. In that work, particles are classified by applying a cut in the space of accretion time a acc vs. mean radial velocity ⟨ v r ⟩ [see 19, for details]. The fiducial cut is obtained by minimizing the number of particles in a narrow band around a straight line in this space, as illustrated in Figure 11. We have slightly modified our orbiting/infall cut  \nTABLE I. Maximum likelihood estimates of the smooth model parameters for the orbiting, infalling and ξ hm profiles. The confidence intervals ( ≈ 68% ) have been set to the largest absolute value of the two bounds.  \nby moving the 'cut' line up or down by the width of the band. We anticipate that this range of variation is roughly comparable to the inherent systematic uncertainty in our orbiting/infall split method.  \nFigure 12 compares the orbiting and infall correlation function contributions to ξ hm from: 1) the best-fit model to the full halo-mass correlation function data; and 2) the three orbiting/infall particle splits (fiducial and ± 1 σ ). We can see that systematic uncertainties in particle classification hardly impact the orbiting profiles. By contrast, ξ inf is clearly sensitive to the details of the orbiting/infall split. Reassuringly, our best-fit smooth infall model to ξ hm is within the range of values spanned by the three orbiting/infall splits we considered except for our highest mass bin. There is a clear trend with mass, such that high-mass haloes are better fit by an orbiting/infall split closer to the bottom of the band shown in Figure 11, whereas low-mass haloes favor a split closer to the top of the band. Based on these results, we conclude that the differences in Figure 10 are due to systematic errors in the particle classification.", 'C. The Halo-Mass Correlation Function of Rockstar Haloes': "Throughout this work we focused on fitting the halo-mass correlation function of dynamical haloes. However, the basic insight that particles in the vicinity of a halo split into orbiting and infalling particles is agnostic to the halo definition. The question arises as to whether our model for the halo-mass correlation function can accurately describe traditional halo catalogues as well. Here, we test whether our model can accurately describe the halo-mass correlation function in the original ROCKSTAR halo catalog used in Garc'ıa et al. [19] as the basis for the dynamical halo catalogue employed in this work.  \nThe two catalogues differ in two ways: 1) our fiducial catalogue defines halo masses in terms M orb rather than M 200 m , which scatters haloes between bins; and 2) we re-percolated the ROCKSTAR algorithm by demanding dark matter particles can only orbit one halo at a time. This condition removes a small fraction of the original ROCKSTAR haloes from the catalogue. For details, we refer the reader do Garc'ıa et al. [19].  \nWe assign every ROCKSTAR halo a mass M 200m , and bin the haloes accordingly. We then fit the resulting halo-mass correlation function using our orbiting+infall model. Unlike for dynamical haloes, however, we do not require that the orbiting profile integrates to the halo mass M 200m , which adds an additional parameter to our fit. Other than that, our model for ROCKSTAR haloes is identical to that employed in our description of dynamical haloes, except that all mass-dependent quantities are assumed to vary with M 200m rather than M orb .  \nOur results are summarized in Figure 13. The left-most panel shows the halo-mass correlation function of ROCKSTAR haloes binned by M 200m , along with our best fit model. The residuals are shown in the corresponding bottom-left panel. Evidently, our model can accurately describe the halo-mass correlation function of ROCKSTAR haloes. The excess variance parameter δ characterising the accuracy of our model fits is δ = 1 . 0% .  \nFigure 13 also compares the performance of our model to the fitting functions proposed by Diemer [18]. The mass density of orbiting particles is modelled as an exponential ρ ( r ) = ρ s e S ( r ) , where:  \nS ( r ) = -2 α [( r r s ) α -1 ] -1 β [ ( r r t ) β -( r s r t ) β ] (18)  \nand { ρ s , α, β, r s , r t } are fit parameters. As discussed in Diemer [18], we also fix α = 0 . 18 and fit for the other four parameters.  \nFIG. 10. 68% and 95% posterior distributions for the parameters our smooth model for the halo-mass correlation function ξ hm . The grey-solid contours come from fitting ξ hm , whereas the red-dashed and blue-dotted contours show the fits to the orbiting and infall components when particles are separated using the criterion of Garc'ıa et al. [19].  \n<!-- image -->  \nTheir infalling model is a power law with a smooth transition given by:  \nρ ( r ) = ¯ ρ m   δ 1 √ ( δ 1 /δ max ) 2 +( r/R ) 2 s -1   (19)  \nwhere { δ 1 , δ max , s } are fit parameters, R is a pivot radius chosen to be R = R 200 m , and ¯ ρ m is the average mean density of  \nthe Universe.  \nWe fit the total density field ρ = ρ orb + ρ inf for each mass bin individually and impose the same priors as in Diemer [18] (see their Table 1). The rightmost panel in Figure 13 shows our fits to the total matter density profile ρ ( r ) . The resulting fits are of comparable quality to ours over the radial range for which the Diemer [18] model was defined ( r/R 200 m ∈ [0 . 1 , 10] ). The  \nFIG. 11. Criteria for splitting particles into orbiting and infalling. The solid cyan line generates the fiducial split outlined in Garc'ıa et al. [19]. The band is at ± 0 . 05 from the fiducial line. For definitions of the truncation mass M t and radius r t see Garc'ıa et al. [13, 19]. We adopt the width of these bands as a rough estimate of the systematic uncertainty associated with orbiting/infall split methodology of Garc'ıa et al. [19].  \n<!-- image -->  \nFIG. 12. Orbiting density profile ρ orb and infall correlation function ξ inf as determined using three different orbiting/infall splits. Each column corresponds to a different mass bin, with mass increasing from left to right. The black data points correspond to our fiducial split; the triangular points to a 1 σ upwards displacement from the fiducial line; and the star points, to a 1 σ downwards displacement. The dotted line is the best fit model recovered from fitting ξ hm directly, and is agnostic to the particle classification scheme.  \n<!-- image -->  \ncentral panel shows this same data, but now with r 2 ξ hm in the y-axis. In this space the fitting formula from Diemer [18] fails at large radii. This is expected. Diemer [18] characterised the halo density profile, not the halo-mass correlation function. Since ρ ( r ) = ¯ ρ m (1 + ξ hm ) , one can achieve high accuracy in  \nρ ( r ) even if the fractional accuracy in ξ hm is low so long as | ξ hm | ≪ 1 .  \nFIG. 13. Halo-mass correlation function for ROCKSTAR haloes binned by M 200m mass. The errors bars are jackknife. The horizontal dotted lines in the residual panels correspond to 2% and 5%. Left: The best-fit model from this work. Our model uses 9 free parameters to describe all 9 curves. Middle and Right panel: Best fit models using the density profiles proposed by Diemer [18]. The latter model was designed to accurately describe the total density profile (right), but not necessarily the halo-mass correlation function (center). Our best fit model has 9 free parameters shared across all 9 curves, while the Diemer [18] has 7 free parameters per curve. Future work will dramatically reduce the number of free parameters in the Diemer [18] model by characterizing the evolution of the model parameters with mass.  \n<!-- image -->", 'D. The Halo-Mass Power Spectrum': 'We have constructed an accurate model for the halo-mass two-point function in configuration space. We briefly explore the accuracy of our model in Fourier space. Since our model is accurate over the range of scales r ∈ [0 . 1 , 50] h -1 Mpc , we test the accuracy of our model in Fourier space on scales k ∈ [0 . 06 , 30] h Mpc -1 (i.e., we identify k ≈ π/r ). We measure the halo-mass power spectrum by pixelizing the box using a grid with N grid = 512 grid points per axis. The halo and matter fields at each pixel are calculated using a cloud-in-cell algorithm. We then take the Fourier transform of the corresponding density fluctuations, and average over narrow shells in k -space. The grid spacing limits the accuracy of this algorithm to scales larger than the Nyquist frequency k Ny = πN grid /L box . For our chosen grid and simulation volume, this corresponds to a relatively small k value of k Ny ≈ 1 . 6 h -1 Mpc .  \nWe overcome this limitation by using the algorithm described in [38]. Briefly, the box is folded up on itself, thereby halving the Nyquist frequency with each fold. For our analysis, we chose to perform five folds of the box, which increases the Nyquist frequency to k Ny ≈ 51 . 5 ( h -1 Mpc) -1 .  \nThe solid lines in the left panel in Figure 14 show the measured power spectrum for all mass bins up to twice the Nyquist frequency of the folded box. The corresponding dashed lines show the Fourier transform of our best-fit model in configuration space. Scales smaller than the Nyquist frequency are shown as a dark grey band, while scales between half the Nyquist frequency and the Nyquist frequency are shown in light grey.  \nFigure 14 demonstrates that there is excellent agreement between the simulation data and the Fourier transform of our model on scales k ∈ [0 . 06 , 6 . 0] h Mpc -1 . The top right panel of Figure 14 shows the ratio of the data to our model. We average this ratio over comparatively broad bins in k -space to produce the bottom-right panel. The error bars are the error on the mean in each of these broad bins. The blue band averages the ratios recovered in each of our different mass bins using inverse variance weights. The width of the band is defined as the standard deviation of ratios as measured using all 9 mass bins. We see that our model for P hm becomes highly biased at high k in a mass-dependent way: while our model accurately describes high mass bins, our model overestimates the power at high k for our low mass bins. It is possible that this overestimate is due to the resolution limit of the simulation, a possibility that we will address in future work. For now, we quantify the accuracy of our model for k ∈ [0 . 06 , 0 . 6] h Mpc -1 by  \n<!-- image -->  \nFIG. 14. Halo-mass power spectrum of dynamical haloes. Left: The solid line is the power spectrum measured up to twice the Nyquist frequency k Ny of the box after five folds (see text for details). The shaded regions represent scales smaller than half the Nyquist frequency (light grey) or the Nyquist frequency (dark grey). Each curve corresponds to a different mass bin as labelled. The dashed line is the Fourier transform of our best fit model for ξ hm ( r ) as per Table I. Top-Right: The ratio between our best fit model and the measured power spectrum. Bottom-Right: The ratio between the data and model fit above, averaging across relatively broad bins in k -space. The blue band is the inverse-variance weighed mean offset, averaged across all mass bins. The width of the band is the standard deviation in this ratio across all mass bins. For reference, the horizontal dotted and dashed lines correspond to 2% and 5% errors respectively.  \n<!-- image -->  \nfitting for the excess variance required for our model to give an adequate description of the data. When we do so, we find the best fit excess variance is 1 . 2% at 1 σ .', 'VI. SUMMARY AND DISCUSSION': "We have proposed a new halo model prescription for the halo-mass correlation function based on the orbiting/infall dichotomy established in Garc'ıa et al. [19]. The existence of this dichotomy implies that the 1-halo/2-halo split in our model is better physically motivated than that employed in standard halo model approaches. Our full model is given by Eq. (17). The density profile ρ orb is given by Eq. (1), while the infall correlation function is given by Eq. (13). The function ξ large ( r ) is defined as per Eq. (10).  \nThe orbiting profile ρ orb is as a running power law with index α ( r ) . This running power law is exponentially truncated at large scales. The model has three free parameters: r h , α ∞ , and a . The parameter r h governs the spatial extent of the halo, and is therefore dubbed the halo radius. The parameter α ∞ is the large-scale limit of α ( r ) , while the parameter a governs the running of α . We find a can be modeled as being mass-independent, making r h the only radial scale associated with our halo profile. Unsurprisingly, r h scales with mass: the more massive a halo, the bigger it is. In addition, α ∞ is close to isothermal, but varies slowly as a function of mass. Note that one could use this variation in slope to define a 'scale radius' as was done in Diemer [18], so either parameterization  \nis valid.  \nThe infall correlation model assumes ξ inf → bξ large where b is the linear halo bias, and ξ large is a suitably defined model for ξ mm on large scales. After considering several options, we settled for the empirical approximation ξ mm ≈ Bξ ZA where B = 0 . 948 at z = 0 . This approximation is accurate to better than 1% on scales r ∈ [30 , 50] h -1 Mpc , and extrapolates reasonably well to small scales. We further modify the basic model ξ inf = bξ large by incorporating: 1) a newly discovered scaling of the infall correlation function (see Figure 7); and 2) a scale-dependent bias β ( r ) that boosts the clustering amplitude around haloes on r ≲ 10 h -1 Mpc scales. The latter is modeled as a cored power law, where the core scale can be associated with the halo radius r h .  \nOur final description of the infall profile for a single mass bin has three mass-dependent parameters ( b , η , and γ ) and two population parameters r inf and µ . The parameter b is the linear halo bias. The parameter η accounts for the newly discovered scaling of ξ inf shown in Figure 7. Finally, the parameters r inf and γ are the amplitude and slope of the cored power law defining the scale-dependent bias. The core radius is set to µr h , where r h is the halo radius. We expect that the amplitude of the small-scale dependent bias r inf must be related to nonlinear structure growth, so we anticipate k NL r inf ≈ 1 , where k NL is the non-linear scale defined by ∆ 2 lin ( k NL ) = 1 . Empirically, we recovered k NL r inf ≈ 2 / 3 .  \nIt is worth commenting on the differences between the model presented in this work and that of Garc'ıa et al. [13] and Zhou and Han [16]. Both of these works defined haloes using phys-  \nally motivated radial apertures. Moreover, both works explicitly incorporated the impact of halo exclusion into their models for ξ hm . A consequence of this type of modelling is that the one- and two-halo terms overlap: the one-halo term extends into the two-halo regime, and vice-versa. This is a necessary ingredient for achieving high-accuracy. In those works, this overlap comes about in an ad-hoc way, whereas in our model it is a natural consequence of our definition. As a specific example, the two-halo term in the model of Garc'ıa et al. [13] does not extend as far into the one-halo regime as our infall model does. This difference is purely cosmetic if one is agnostic to the physical interpretation of the fitting function. However, we have demonstrated that our one and two halo terms have a clear physical interpretation that is lacking in the work of Garc'ıa et al. [13].  \nOne possible source of concern is that we have made no attempt to incorporate halo-exclusion effects. However, we do not believe this is problematic. Although the current data set is expected to be impacted by halo-exclusion effects, we anticipate these will cease to exist in the future. When haloes are defined as collections of orbiting particles, the restriction that haloes cannot overlap vanishes: if a small halo is falling into a larger halo, and the small halo has not yet experienced its first pericentric passage, then that halo will still be considered a parent halo. This is true even if the smaller halo is entirely contained within the radius of the larger halo. That is, a dynamical halo finder will produce halo catalogs with no halo exclusion effects. For this reason, we have chosen to ignore halo exclusion effects in our current model. We will present a halo finder that produces halo catalogues with no halo exclusion effects in future work.  \n- [1] A. Cooray and R. Sheth, ßPhys. Rep. 372 , 1 (2002), astroph/0206508.\n- [2] M. Asgari, A. J. Mead, and C. Heymans, The Open Journal of Astrophysics 6 , 39 (2023), arXiv:2303.08752 [astro-ph.CO].\n- [3] E. Hayashi and S. D. M. White, MNRAS 388 , 2 (2008), arXiv:0709.3933 [astro-ph].\n- [4] J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490 , 493. (1997).\n- [5] J. Einasto, Trudy Astrofizicheskogo Instituta Alma-Ata 5 , 87 (1965).\n- [6] F. C. van den Bosch, S. More, M. Cacciato, H. Mo, and X. Yang, MNRAS 430 , 725 (2013), arXiv:1206.6890 [astro-ph.CO].", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'ACKNOWLEDGMENTS': 'E. Salazar is supported by NSF grant 3047980. In addition to NSF 3047980, E. Rozo is also funded by DOE grant DE-SC0009913. NK acknowledges support from NSF award AST-2108126.', 'Appendix A: Orbiting Mass Correction': "By definition, the orbiting density profile ρ ( r ) integrates to the total orbiting mass of a halo. Suppose the density profile in simulations is well described by our model ρ model ( r ) at scales r > R , whereas ρ data ( r ) differs from ou rmodel on scales below R (e.g. due to resolution effects). By integrating our model over all space we must account for the deficit or excess in mass due to the extrapolation of the model to small scales. We define the orbiting mass correction to be:  \n∆ M orb = ∫ r<R ( ρ model -ρ data ) dV (A1)  \nWe have chosen R to be six times the simulation softening length. Whenever we refer to the 'orbiting mass' of stacked profiles in the text, it is always including the correction to the ensemble mass average as:  \nM orb = ⟨ M orb ⟩ +∆ M orb (A2)  \nThis amounts to at most a 2% correction to the halo mass for the lowest mass haloes.  \n- [7] B. Diemer and A. V. Kravtsov, Astrophys. J. 789 , 1 (2014), arXiv:1401.1216.\n- [8] S. Adhikari, N. Dalal, and R. T. Chamberlain, Journal of Cosmology and Astro-Particle Physics 2014 , 019 (2014), arXiv:1409.4482 [astro-ph.CO].\n- [9] S. More, H. Miyatake, M. Takada, B. Diemer, A. V. Kravtsov, N. K. Dalal, A. More, R. Murata, R. Mandelbaum, E. Rozo, E. S. Rykoff, M. Oguri, and D. N. Spergel, Astrophys. J. 825 , 39 (2016), arXiv:1601.06063 [astro-ph.CO].\n- [10] P. Mansfield, A. V. Kravtsov, and B. Diemer, Astrophys. J. 841 , 34 (2017), arXiv:1612.01531 [astro-ph.CO].\n- [11] B. Diemer, Astrophys. J. 909 , 112 (2021), arXiv:2007.10992", '[astro-ph.CO].': "- [12] B. Diemer, ApJS 231 , 5 (2017), arXiv:1703.09712.\n- [13] R. Garc'ıa, E. Rozo, M. R. Becker, and S. More, MNRAS 505 , 1195 (2021), arXiv:2006.12751 [astro-ph.CO].\n- [14] M. Fong and J. Han, MNRAS 503 , 4250 (2021), arXiv:2008.03477 [astro-ph.CO].\n- [15] A. J. Cuesta, F. Prada, A. Klypin, and M. Moles, MNRAS 389 , 385 (2008), arXiv:0710.5520.\n- [16] Y. Zhou and J. Han, MNRAS 525 , 2489 (2023), arXiv:2303.10886 [astro-ph.CO].\n- [17] B. Diemer, MNRAS 513 , 573 (2022), arXiv:2112.03921 [astroph.CO].\n- [18] B. Diemer, MNRAS 519 , 3292 (2023), arXiv:2205.03420 [astro-ph.CO].\n- [19] R. Garc'ıa, E. Salazar, E. Rozo, S. Adhikari, H. Aung, B. Diemer, D. Nagai, and B. Wolfe, MNRAS 521 , 2464 (2023), arXiv:2207.11827 [astro-ph.CO].\n- [20] J. Diemand and M. Kuhlen, ApJL 680 , L25 (2008), arXiv:0804.4185.\n- [21] V. Desjacques, D. Jeong, and F. Schmidt, ßPhys. Rep. 733 , 1 (2018), arXiv:1611.09787 [astro-ph.CO].\n- [22] A. Banerjee, S. Adhikari, N. Dalal, S. More, and A. Kravtsov, JCAP 2020 , 024 (2020), arXiv:1906.12026 [astro-ph.CO].\n- [23] V. Springel, MNRAS 364 , 1105 (2005), astro-ph/0505010.\n- [24] V. Springel, N-GenIC: Cosmological structure initial conditions, Astrophysics Source Code Library, record ascl:1502.003 (2015), ascl:1502.003.\n- [25] P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, Astrophys. J. 762 , 109 (2013), arXiv:1110.4372 [astro-ph.CO].\n- [26] P. Mansfield and C. Avestruz, MNRAS 500 , 3309 (2021), arXiv:2008.08591 [astro-ph.CO].\n- [27] S.-F. Chen, Z. Vlah, and M. White, JCAP 2020 , 062 (2020), arXiv:2005.00523 [astro-ph.CO].  \n- [28] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Astrophys. J. 761 , 152 (2012), arXiv:1208.2701 [astro-ph.CO].\n- [29] R. E. Smith and R. E. Angulo, MNRAS 486 , 1448 (2019), arXiv:1807.00040 [astro-ph.CO].\n- [30] D. J. Bartlett, B. D. Wandelt, M. Zennaro, P. G. Ferreira, and H. Desmond, arXiv e-prints , arXiv:2402.17492 (2024), arXiv:2402.17492 [astro-ph.CO].\n- [31] K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, Astrophys. J. 780 , 111 (2014), arXiv:1304.7849 [astro-ph.CO].\n- [32] T. Nishimichi, M. Takada, R. Takahashi, K. Osato, M. Shirasaki, T. Oogi, H. Miyatake, M. Oguri, R. Murata, Y. Kobayashi, and N. Yoshida, Astrophys. J. 884 , 29 (2019), arXiv:1811.09504 [astro-ph.CO].\n- [33] R. E. Angulo, M. Zennaro, S. Contreras, G. Aric'o, M. Pellejero-Iba˜nez, and J. Stucker, MNRAS 507 , 5869 (2021), arXiv:2004.06245 [astro-ph.CO].\n- [34] K. R. Moran, K. Heitmann, E. Lawrence, S. Habib, D. Bingham, A. Upadhye, J. Kwan, D. Higdon, and R. Payne, MNRAS 520 , 3443 (2023), arXiv:2207.12345 [astro-ph.CO].\n- [35] R. A. Porto, L. Senatore, and M. Zaldarriaga, Journal of Cosmology and Astroparticle Physics 2014 (05), 022-022.\n- [36] Z. Vlah, M. White, and A. Aviles, Journal of Cosmology and Astroparticle Physics 2015 (09), 014-014.\n- [37] F. Villaescusa-Navarro, C. Hahn, E. Massara, A. Banerjee, A. M. Delgado, D. K. Ramanah, T. Charnock, E. Giusarma, Y. Li, E. Allys, A. Brochard, C. Uhlemann, C.-T. Chiang, S. He, A. Pisani, A. Obuljen, Y. Feng, E. Castorina, G. Contardo, C. D. Kreisch, A. Nicola, J. Alsing, R. Scoccimarro, L. Verde, M. Viel, S. Ho, S. Mallat, B. Wandelt, and D. N. Spergel, ApJS 250 , 2 (2020), arXiv:1909.05273 [astro-ph.CO].\n- [38] S. Colombi, A. Jaffe, D. Novikov, and C. Pichon, MNRAS 393 , 511 (2009), arXiv:0811.0313 [astro-ph]."}