A dynamics-based density profile for dark haloes – II. Fitting function

Diemer, Benedikt

doi:10.1093/mnras/stac3778

ABSTRACT

The density profiles of dark matter haloes are commonly described by fitting functions such as the NFW or Einasto models, but these approximations break down in the transition region where haloes become dominated by newly accreting matter. Here, we present a simple accurate new fitting function that is inspired by the asymptotic shapes of the separate orbiting and infalling halo components. The orbiting term is described as a truncated Einasto profile, ρ_orb ∝ exp [ − 2/α (r/r_s)^α − 1/β (r/r_t)^β], with a five-parameter space of normalization, physically distinct scale and truncation radii, and α and β, which control how rapidly the profiles steepen. The infalling profile is modelled as a power law in overdensity that smoothly transitions to a constant at the halo centre. We show that these formulae fit the averaged total profiles in simulations to about 5 per cent accuracy across almost all of an expansive parameter space in halo mass, redshift, cosmology, and accretion rate. When fixing α = 0.18 and β = 3, the formula becomes a three-parameter model that fits individual haloes better than the Einasto profile on average. By analogy with King profiles, we show that the sharp truncation resembles a cut-off in binding energy.

methods: numerical, dark matter, large-scale structure of Universe

1 INTRODUCTION

The shapes that make up the cosmic web of dark matter are generally too complex to be described mathematically (e.g. Doroshkevich & Shandarin 1978; Bond, Kofman & Pogosyan 1996). However, the densest structures, or haloes, collapse into relatively uniform and roughly spherical shapes. Describing haloes analytically is critical because they contain about half the dark and most of the observable baryons in the Universe. The most common descriptions are spherically averaged density profiles, which conveniently summarize simulation results (Dubinski & Carlberg 1991), serve as the basis of galaxy formation models (Wechsler & Tinker 2018), underlie descriptions of large-scale structure via the halo model (Cooray & Sheth 2002), are used to fit observed profile data (e.g. Courteau et al. 2014), and provide a parameter space in which to compare such observations to theory and simulations (e.g. Umetsu 2020; Eckert et al. 2022).

To facilitate these applications, numerous fitting functions have been proposed, although many were originally intended for galaxies rather than haloes. The three-parameter (Einasto 1965, 1969) profile features a slope that smoothly steepens with radius and thus contains a finite mass. Many later proposals are essentially double power laws that smoothly transition between two slopes around a scale radius, r_s (Dehnen 1993). This group includes the profiles of Jaffe 1983 (with inner slope −2 and outer slope −4), Hernquist 1990 (−1 and −4), Burkert 1995 (0 and −3), Navarro, Frenk & White 1995, 1996, 1997 (hereafter NFW, −1 and −3), and Moore et al. 1999 (−1.5 and −3, see also Ghigna et al. 2000). Numerous studies have found the extra parameter of the Einasto profile to be justified by a superior fit to simulation results (Navarro et al. 2004, 2010; Graham et al. 2006; Merritt et al. 2006; Gao et al. 2008; Stadel et al. 2009; Ludlow et al. 2011; Wang et al. 2020). Naturally, the fit of power-law models can be improved by adding more slope parameters (Zhao 1996; Dekel et al. 2017; Freundlich et al. 2020). Finally, an entire different class of models can be derived by assuming a distribution of particle energies and computing the corresponding density structure (King 1966; Shapiro, Iliev & Raga 1999; Hjorth & Williams 2010; Pontzen & Governato 2013).

All of these models were intended to fit density profiles out to roughly the virial radius, a focus that makes sense given that galaxy formation happens at much smaller radii. For example, one of the key debates has been whether the central slope of the profiles is characteristic of a flat core or a power-law cusp (e.g. de Blok 2010; Teyssier et al. 2013; Di Cintio et al. 2014; Genina et al. 2018). A second reason is that the nature of the profiles changes fundamentally beyond the virial radius, where the profile becomes dominated by particles that are falling into the halo for the first time. In this transition region, the profile shape is the result of a complex interplay between the orbiting and infalling terms, an inherently dynamical distinction that cannot easily be made based on the profiles alone (Fukushige & Makino 2001, Diemand & Kuhlen 2008). As a result, conventional fitting functions trained on the orbiting (or one-halo) term struggle to fit simulated or observed profiles in this regime (Becker & Kravtsov 2011; Oguri & Hamana 2011). While some models for the outer profiles have been proposed (Prada et al. 2006; Betancort-Rijo et al. 2006; Tavio et al. 2008; Baltz, Marshall & Oguri 2009), they fail to capture the detailed profile shapes (Diemer & Kravtsov 2014, hereafter DK14).

Despite these complications, the transition region and the outer profiles have recently gained renewed relevance because they are much less affected by baryons and because they are sensitive to otherwise inaccessible halo properties. For example, the mass accretion rate sets the position of the edge of the orbiting term, which is also known as the splashback radius (DK14, Adhikari, Dalal & Chamberlain 2014; More, Diemer & Kravtsov 2015). Moreover, the outer profiles have become observationally accessible via the density of satellite galaxies and the weak lensing signal around clusters (e.g. More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019; Murata et al. 2020; Bianconi et al. 2021). For their inferences, all of these works have relied on the fitting function of DK14, an Einasto profile multiplied by a power-law cutoff term with variable truncation radius, sharpness, and asymptotic slope (a total of 6 free parameters). This flexible function fits averaged profiles to 5–10 per cent accuracy including the transition region, but its design suffered from the same fundamental issue as previous models: at its truncation, the orbiting term is concealed by the infalling term, and its asymptotic shape was thus unknown. As a result, the DK14 fitting function is hard to interpret physically. First, the large number of parameters causes well-known degeneracies that necessitate informative priors (DK14; Baxter et al. 2017; Chang et al. 2018; Umetsu & Diemer 2017). Second, the steepening term approaches a somewhat arbitrary slope that does not actually reflect the shape of the orbiting term. Some observational works have extrapolated this shape based on DK14 fits (e.g. Baxter et al. 2017; Shin et al. 2021), but it is not clear that the results are physically meaningful. Third, the slope of the DK14 profile is a complex function, which makes it difficult to establish the relationship between the profile parameters and physical features such as the splashback radius.

In this paper, we set out to design a more physically motivated fitting function. For the first time, we can draw inspiration from profiles that have been dynamically split into orbiting and infalling particles (using the algorithm presented in Diemer 2022, hereafter Paper I; see also García et al. 2022). Similar splits have recently been explored in observational data using galaxy properties such as colour (Baxter et al. 2017; Adhikari et al. 2021; O’Donnell, Behroozi & More 2021; Shin et al. 2021; Aung et al. 2022; Dacunha et al. 2022; O’Neil et al. 2022), necessitating a fitting function that captures the shapes of the orbiting and infalling terms. We require that this fitting function should (1) describe the profiles accurately even in the transition region, (2) be flexible enough to apply to individual haloes and to stacks with different selection criteria, (3) work across all halo masses, redshifts, cosmologies, and accretion rates, (4) rely on as few free parameters as possible with clear, physical interpretations, and (5) exhibit minimal parameter degeneracies. We achieve these goals with a 5-parameter Einasto-like form with a truncated exponent. Our aim is not necessarily to achieve more accurate fits than DK14 because there are complex variations in the averaged profiles that make it difficult to systematically improve on a 5–10 per cent fit. Moreover, baryons affect the profiles at roughly this level even at large radii (e.g. Velliscig et al. 2014; Schneider et al. 2019). Thus, we instead focus on creating a parameter space that can be used to meaningfully describe observed and simulated profiles. In this second paper of the series, we present and test the fitting function. In Diemer (in preparation; hereafter Paper III), we analyse the resulting parameter space and connect it to halo physics such as accretion histories.

The paper is structured as follows. In Section 2, we mathematically describe previous fitting functions and the new models. We summarize our numerical methods in Section 3, largely referring the reader to Paper I. We assess the quality of our fits to stacked and individual halo profiles in Section 4, and we compare to previously proposed models in Section 5. Section 6 summarizes our results. We discuss an alternative model variant in Appendix B and derive additional mathematical properties in Appendices C and D. Supplementary figures are provided online on the author’s website at benediktdiemer.com/data. Our fitting function is implemented in the publicly available code colossus (Diemer 2018).

Throughout the paper, we follow the notation of Paper I. We normalize halo radii by R_200m, which encloses an average of 200 times the mean density of the Universe, ρ_m, including all particles (bound or unbound). We denote the mass enclosed within this radius as M_200m. We mostly express mass as peak height, ν ≡ δ_c/σ(M_200m), where δ_c = 1.686 is the critical density for collapse (Gunn & Gott 1972) and σ(M_200m) is the variance of the linear power spectrum on the Lagrangian scale corresponding to M_200m. We define the logarithmic mass accretion rate Γ ≡ Δln M_200m/Δln a, measured over one dynamical time (which we define as the mass-independent halo crossing time, about 5 Gyr at z = 0). More details on these calculations are given in Paper I.

2 PROFILE MODELS

In this section, we introduce new fitting functions for the orbiting and infalling components of density profiles, which can be combined to fit the entire profile. We describe the general logic of constructing models for the orbiting term in Section 2.1, previously proposed models in Sections 2.2 and 2.3, and our new model in Section 2.4. We introduce a slightly different version of the same model in Appendix B. We also discuss previous and new models for the infalling profile in Sections 2.5 and 2.6. The derivatives of the models with respect to their parameters are given in Appendix C.

2.1 Generalized exponential profiles

Our new function for the orbiting term can be thought of as a generalization of the Einasto profile, where the slope smoothly changes with radius. We write such profiles as the exponential of a radial function S(r),

$$\begin{eqnarray*} \rho _{\rm orb} = \rho _{\rm s}e^{S(r)} \, . \end{eqnarray*}$$

(1)

This family of profiles is more intuitively understood by considering the logarithmic slope γ(r),

$$\begin{eqnarray*} \gamma \equiv \frac{{\rm d} \ln \rho }{{\rm d} \ln r} = r \frac{{\rm d} S}{{\rm d} r} \, . \end{eqnarray*}$$

(2)

This expression demonstrates the power of equation (1): we can construct a density profile by defining its logarithmic slope and integrating to get S(r),

$$\begin{eqnarray*} S(r) = \int \frac{\gamma (r)}{r} {\rm d}r \, . \end{eqnarray*}$$

(3)

Clearly, one condition for the fitting function to be useful is that γ/r must integrate to a reasonably simple, analytical expression. For example, exponential cutoffs in the slope such as γ ∝ exp (r/r_t) are problematic because their integrals (if solvable) involve Gamma functions or other complex mathematical expressions.

2.2 The Einasto profile

We begin by reviewing the Einasto profile in light of the discussion above. Its defining feature is a slope of γ = −2(r/r_s)^α. The parameter α determines how rapidly the profile steepens, and the meaning of the scale radius is that γ(r_s) = −2. Moreover, the Einasto profile leads to a core at small radii, γ(0) = 0. The slope integrates to S(r) = −(2/α)(r/r_s)^α + C. We can choose to set C = 0, in which case ρ_s describes the density at r = 0, but the profile is more commonly written with C = 2/α,

$$\begin{eqnarray*} S(r) = -\frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha - 1 \right] \, , \end{eqnarray*}$$

(4)

such that S(r_s) = 0 and thus ρ(r_s) = ρ_s. In practice, we find that the C ≠ 0 version of the Einasto profile behaves better in least-squares fits because it does not allow a degeneracy between ρ_s, r_s, and α (although the latter two can still be degenerate at fixed ρ_s).

2.3 The DK14 profile

We briefly review the DK14 model because we will later use its performance as a benchmark for our new model. The model is

$$\begin{eqnarray*} \rho _{\rm DK14} = \rho _{\rm Einasto} \left[1 + \left(\frac{r}{r_{\rm t}} \right)^\beta \right]^{-\frac{\gamma }{\beta }} \, , \end{eqnarray*}$$

(5)

where the steepening term suppresses the density at r ≳ r_t as a power law with slope −γ, while β governs the sharpness of the truncation. Combined with the Einasto parameters, this function has six free parameters, but DK14 suggested fixing (β, γ) to (4, 8) in fits to mass-selected samples and to (6, 4) when the sample is also selected by mass accretion rate (they did not investigate fits to individual haloes). The main disadvantage of equation (5) is that the β and γ parameters do not have a clear physical meaning because the asymptotic slope of the overall orbiting profile depends on r_s, r_t, α, and γ in a complex fashion. Setting β and γ to fixed values reduces the profile to a four-parameter fit with a meaningful truncation radius parameter r_t, but it no longer allows for a varying shape of the truncation term.

2.4 The new model: truncated exponential

In Paper I, we confirmed that the orbiting profile has two characteristic radial scales, which can be captured by the scale and the truncation radii. Moreover, we found that the sharpness of the truncation varies with halo properties and that the profiles decline with an exponential-like, accelerating slope (rather than a constant, power-law slope). Inspired by these observations, we construct a model by adding a truncation term into the slope,

$$\begin{eqnarray*} \gamma (r) = -2 \left(\frac{r}{r_{\rm s}} \right)^\alpha - \left(\frac{r}{r_{\rm t}} \right)^\beta \, . \end{eqnarray*}$$

(6)

While similar in spirit, this model significantly differs from the ‘truncated Sersic’ and ‘broken exponential’ models that are popular in the observational literature (e.g. Peng et al. 2010; Erwin 2015). We find the corresponding profile by integrating γ/r,

$$\begin{eqnarray*} S(r) = -\frac{2}{\alpha } \left(\frac{r}{r_{\rm s}} \right)^\alpha -\frac{1}{\beta } \left(\frac{r}{r_{\rm t}} \right)^\beta + C \, . \end{eqnarray*}$$

(7)

For simplicity, we could set C = 0, in which case ρ_s would become the density at the centre. A more elegant form is obtained by setting the integration constant such that ρ(r_s) = ρ_s,

$$\begin{eqnarray*} C = \frac{2}{\alpha } + \frac{1}{\beta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \, , \end{eqnarray*}$$

(8)

and thus

$$\begin{eqnarray*} S(r) = -\frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha - 1 \right] -\frac{1}{\beta } \left[ \left(\frac{r}{r_{\rm t}} \right)^\beta - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \right] \, . \end{eqnarray*}$$

(9)

We use this parametrization throughout because the density at r_s is numerically constrained, whereas the central density represents an extrapolation. However, we give derivatives with respect to the free parameters for both variants in Appendix C2. In Fig. 1, we visually explore how the free parameters affect the shape of our model. Throughout the rest of the paper, we show that equation (9) describes the orbiting profiles with great accuracy.

Figure 1.

Effect of the free parameters in our model for the orbiting density profile (top) and its logarithmic slope (bottom). Each column shows variations of a single parameter; ρ_s is omitted as it corresponds to a simple normalization. All profiles are normalized by some radius R and integrate to the same M(R). The black lines show the same fiducial profile in each column (r_s = 0.2 R, α = 0.15, r_t = R, β = 4), and the dotted vertical lines highlight the fiducial r_s and r_t. The horizontal dashed lines mark the mean density in the top panels and a slope of −2 in the bottom panels. All slopes converge to γ(0) = 0 at the halo centre, but we would need to extend the plot to much smaller radii to see this effect.

Open in new tab Download slide

The model does, however, have two minor shortcomings. First, the truncation term ‘breaks’ the meaning of the scale radius because the slope at r_s is now γ(r_s) = −2 − (r_s/r_t)^β instead of −2. The difference is small in most cases since r_t > r_s and generally β ≫ 1, but it can manifest itself for extreme parameter values (red lines in the right two columns of Fig. 1). We note that the γ(r_s) = −2 condition will likely be violated regardless as soon as an infalling profile is added. None the less, in Appendix B, we present a model variant that enforces γ(r_s) = −2 at the expense of an additional term in S(r). Both models give the same fit and best-fitting parameters for virtually all profiles, but the variant can be preferable in cases where r_s and α are strongly degenerate. For either model, the scale radii (and thus concentrations) are directly comparable to those from Einasto fits, but we expect small systematic shifts similar in scale to the differences between NFW and Einasto concentrations (Dutton & Macciò 2014) or differences due to the fitting procedure (Dooley et al. 2014).

The second shortcoming is that the integrals of ρ, namely the enclosed mass M(r), the projected density Σ(R), and the lensing signal ΔΣ, cannot be computed analytically. Similar issues hamper even the simpler Einasto profile, where M(r) is a relatively complicated analytical expression (Cardone, Piedipalumbo & Tortora 2005; Retana-Montenegro et al. 2012) and Σ(R) can only be approximated (Dhar & Williams 2010; Dhar 2021). We will investigate similar approximate solutions for the new model in future work, and we have provided fast numerical solution in the colossus code (Diemer 2018).

2.5 The infalling profile: previous models

Before shells of dark matter begin to cross, the overdensity due to non-linearly infalling matter is expected to roughly follow a power law with a slope of −3/2 (Bertschinger 1985). At large radii, the statistical contribution from large-scale structure comes to dominate, but this transition happens at radii larger than the 10 R_200m we consider (Paper I). Thus, it is not surprising that DK14 found the profiles to be well described by a simple power law in overdensity,

$$\begin{eqnarray*} \rho (r) = \rho _{\rm m}(z) \left[ \delta _{\rm 1}\left(\frac{r}{R} \right)^{-s} + 1 \right] \, , \end{eqnarray*}$$

(10)

where δ₁ represents the normalization at radius R and −s the slope. We will set R = R_200m throughout, but any other pivot radius could be chosen instead. One issue with equation (10) is that it can reach arbitrarily high values at small r, which is clearly unphysical. We can prevent this artefact by introducing a maximum overdensity (e.g. Diemer 2018),

$$\begin{eqnarray*} \rho (r) = \rho _{\rm m}\left(\delta _{\rm 1}\left[ \frac{\delta _{\rm 1}}{\delta _{\rm max}} + \left(\frac{r}{R} \right)^{s} \right]^{-1} + 1 \right) \, . \end{eqnarray*}$$

(11)

This expression approaches ρ → ρ_m(δ_max + 1) as r → 0. The exact value of δ_max does not matter because it is reached at radii where the infalling profile dominates by orders of magnitude.

2.6 New infalling model: power law with smooth transition

Inspecting the infalling profiles in Paper I, we found that they do indeed seem to approach a constant density at small radii, but the transition to this value happens less sharply than suggested by equation (11). We thus introduce a transition smoothness parameter, ζ,

$$\begin{eqnarray*} \rho (r) &= \rho _{\rm m}\left(\delta _{\rm 1}\left[ \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right)^{\frac{1}{\zeta }} + \left(\frac{r}{R} \right)^{\frac{s}{\zeta }} \right]^{-\zeta } + 1 \right) \equiv \rho _{\rm m}\left[ \delta _{\rm 1}Q(r)^{-\zeta } + 1 \right] \, , \nonumber \\ \end{eqnarray*}$$

(12)

where we have defined

$$\begin{eqnarray*} Q(r) \equiv \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right)^{\frac{1}{\zeta }} + \left(\frac{r}{R}\right)^{\frac{s}{\zeta }} \, . \end{eqnarray*}$$

(13)

The shape of this model as a function of its parameters is visualized in Fig. 2. We note that the function’s logarithmic slope,

$$\begin{eqnarray*} \gamma (r) = - \left(1 - \frac{\rho _{\rm m}}{\rho (r)} \right) \frac{s}{Q(r)} \left(\frac{r}{R} \right)^{\frac{s}{\zeta }} \, , \end{eqnarray*}$$

(14)

remains at a value of −s only across a very narrow radial range for most parameter combinations, if at all. The full freedom of this function allows for very accurate fits, but also for a number of pathological cases where not all parameters are well constrained by the simulated profiles. We find that fixing ζ = 0.5 is a compromise that works for the majority of profiles. The resulting function,

$$\begin{eqnarray*} \rho (r) = \rho _{\rm m}\left(\frac{\delta _{\rm 1}}{\sqrt{(\delta _{\rm 1}/ \delta _{\rm max})^{2} + (r / R)^{2s}}} + 1 \right) \, , \end{eqnarray*}$$

(15)

is the infalling profile that we use throughout this paper.

Figure 2.

Same as Fig. 1 but for the new infalling profile model, which smoothly transitions from a power law to an asymptotic central density. The parameters are the normalization at radius R, the slope of the power law, the maximum central overdensity, and the transition smoothness. The δ_max parameter can affect the overdensity at R if it is sufficiently small.

Open in new tab Download slide

2.7 Allowed ranges of parameters

We impose an allowed range on each parameter to avoid extreme unphysical values. These flat priors (in logarithmic space) are listed in Table 1. They are designed to be uninformative, with the exception that they enforce r_s < r_t. Some parameters are fixed when fitting individual halo profiles (Table 1). We discuss the reasons behind the chosen values in detail in Appendix A2. In Paper III, we will also present the distributions of best-fitting parameters for individual and stacked profiles, which could be used as a prior when fitting to observational data. The individual fits can fill out the allowed range of the parameters because not all individual profiles constrain all parameters equally well, whereas the ranges are generous for the averaged profiles.

Table 1.

Open in new tab

Summary of the free parameters of the fitting functions for the orbiting and infalling terms, including their allowed ranges allowed in fits. If the lower and upper bounds are the same, the parameter is held fixed. The last column lists possible modifications in fits to individual profiles. If ‘free’, the same limits apply as for averaged profiles. See Section 2.7 and Appendix A2 for details.

Symbol	Description	Lower	Upper	Ind.
Orbiting term
ρ_s/ρ_m	Overdensity at scale radius	10	10⁷	free
r_s/R_200m	Scale radius	0.01	0.45	free
r_t/R_200m	Truncation radius	0.5	3	<10
α	Radial evolution of slope	0.03	0.4	0.18
β	Sharpness of truncation	0.1	10	3
Infalling term
δ₁	Normalization at R_200m	1	100	free
s	Power-law slope	0.01	4	free
δ_max	Central overdensity	10	2000	free
ζ	Smoothness of transition	0.5	0.5	0.5
Model variations
ρ₀/ρ_m	Overdensity at r = 0	10	10²⁰	–
η	Term in Model B (App. B)	0.1	0.1	0.1

Symbol	Description	Lower	Upper	Ind.
Orbiting term
ρ_s/ρ_m	Overdensity at scale radius	10	10⁷	free
r_s/R_200m	Scale radius	0.01	0.45	free
r_t/R_200m	Truncation radius	0.5	3	<10
α	Radial evolution of slope	0.03	0.4	0.18
β	Sharpness of truncation	0.1	10	3
Infalling term
δ₁	Normalization at R_200m	1	100	free
s	Power-law slope	0.01	4	free
δ_max	Central overdensity	10	2000	free
ζ	Smoothness of transition	0.5	0.5	0.5
Model variations
ρ₀/ρ_m	Overdensity at r = 0	10	10²⁰	–
η	Term in Model B (App. B)	0.1	0.1	0.1

Table 1.

Open in new tab

Summary of the free parameters of the fitting functions for the orbiting and infalling terms, including their allowed ranges allowed in fits. If the lower and upper bounds are the same, the parameter is held fixed. The last column lists possible modifications in fits to individual profiles. If ‘free’, the same limits apply as for averaged profiles. See Section 2.7 and Appendix A2 for details.

Symbol	Description	Lower	Upper	Ind.
Orbiting term
ρ_s/ρ_m	Overdensity at scale radius	10	10⁷	free
r_s/R_200m	Scale radius	0.01	0.45	free
r_t/R_200m	Truncation radius	0.5	3	<10
α	Radial evolution of slope	0.03	0.4	0.18
β	Sharpness of truncation	0.1	10	3
Infalling term
δ₁	Normalization at R_200m	1	100	free
s	Power-law slope	0.01	4	free
δ_max	Central overdensity	10	2000	free
ζ	Smoothness of transition	0.5	0.5	0.5
Model variations
ρ₀/ρ_m	Overdensity at r = 0	10	10²⁰	–
η	Term in Model B (App. B)	0.1	0.1	0.1

Symbol	Description	Lower	Upper	Ind.
Orbiting term
ρ_s/ρ_m	Overdensity at scale radius	10	10⁷	free
r_s/R_200m	Scale radius	0.01	0.45	free
r_t/R_200m	Truncation radius	0.5	3	<10
α	Radial evolution of slope	0.03	0.4	0.18
β	Sharpness of truncation	0.1	10	3
Infalling term
δ₁	Normalization at R_200m	1	100	free
s	Power-law slope	0.01	4	free
δ_max	Central overdensity	10	2000	free
ζ	Smoothness of transition	0.5	0.5	0.5
Model variations
ρ₀/ρ_m	Overdensity at r = 0	10	10²⁰	–
η	Term in Model B (App. B)	0.1	0.1	0.1

3 SIMULATIONS AND METHODS

In this section, we describe our simulations and algorithms. In the interest of brevity, we do not excessively duplicate information from Paper I. We briefly summarize our simulations in Section 3.1 and the profile data in Section 3.2, largely referring the reader to Paper I. We devote more detail to our fitting procedure in Section 3.3 and Appendix A.

3.1 N-body simulations

Our analysis is based on the Erebos suite of dissipationless N-body simulations (Diemer & Kravtsov 2014, 2015), which consist of 14 simulations of 1024³ dark matter particles (see table 1 in Paper I). The suite covers different box sizes and resolutions, as well as two ΛCDM and four self-similar cosmologies. The first ΛCDM cosmology is that of the Bolshoi simulation (Klypin, Trujillo-Gomez & Primack 2011), consistent with WMAP7 (Komatsu et al. 2011), namely a flat ΛCDM cosmology with Ω_m = 0.27, Ω_b = 0.0469, σ₈ = 0.82, and n_s = 0.95. The second is a Planck-like cosmology (Planck Collaboration XVI 2014; Ω_m = 0.32, Ω_b = 0.0491, h = 0.67, σ₈ = 0.834, and n_s = 0.9624).

We also consider four self-similar Einstein-de Sitter universes with power-law initial spectra of slopes n = −1, −1.5, −2, and −2.5. Here, length and time are scale free, meaning that the density profiles are independent of redshift when rescaled by a meaningful radius. The self-similar simulations highlight the impact of the initial power spectrum (Efstathiou et al. 1988; Knollmann, Power & Knebe 2008) and allow us to test a wide range of cosmologies with few simulations because ΛCDM can be seen as an interpolation between different power-spectrum slopes. For example, redshift trends in ΛCDM profiles are mostly trends in n (Paper I), meaning that a good fit of our function to profiles from self-similar cosmologies implies a wide range in redshift.

The power spectra for the ΛCDM simulations were generated using camb (Lewis, Challinor & Lasenby 2000). They were translated into initial conditions using 2LPTic (Crocce, Pueblas & Scoccimarro 2006). All simulations were run with gadget2 (Springel 2005). We use the phase–space halo finder rockstar (Behroozi, Wechsler & Wu 2013a) to identify haloes and subhaloes. We construct merger trees by connecting haloes across time using the consistent-trees code (Behroozi et al. 2013b). The resulting halo catalogues are described in detail in Diemer (2020).

3.2 Dynamically split density profiles

The key new ingredient on which we base our models is the splitting of dark matter particles into infalling and orbiting, the transition between which is defined to occur at a particle’s first pericentre (though alternative definitions exist, e.g. Sugiura et al. 2020; García et al. 2022). The pericentres are reliably detected by a novel algorithm that follows the orbit of each particle in each halo. This algorithm was implemented in the sparta framework (Diemer 2017, 2020) and applied to all 14 Erebos simulations (Paper I). Despite a few ambiguous particle orbits (e.g. due to the limited time resolution of snapshots), the resulting split profiles were shown to be robust.

We apply the same resolution cuts as in Paper I. We consider only host haloes with at least 500 particles within R_200m, and we cut out those parts of profiles that lie within |$r_{\rm min} = {\rm max}[4\epsilon , 0.133\ \Omega _{\rm m}^{1/3} l(z)]$|⁠, where ϵ and l(z) are the comoving force resolution and interparticle spacing of the simulation. These cuts limit the effects of suppressed centripetal forces and two-body scattering (see Appendix A1 of Paper I for details; Power et al. 2003; Ludlow, Schaye & Bower 2019; Mansfield & Avestruz 2021). We also impose a limit on the fraction of a halo’s mass that is unbound, M_{200m, all}/M_{200m, bnd} < 1.5, which excludes haloes that are being disrupted by interactions with neighbours (Paper I). However, we find that the effects of neighbours can still significantly distort the mean orbiting profiles of low-mass haloes near the truncation. Thus, we exclude profile bins where ρ < 0.1 ρ_m from the fits to mean orbiting profiles, although we do plot those bins in the following figures.

We fit both individual halo profiles and averaged (mean and median) profiles. The latter are constructed by considering halo samples defined by a range of peak height (or mass), a redshift, and possibly mass accretion rate. Each sample combines haloes from different Erebos simulation boxes. In the self-similar simulations, there is no physical time, so that each sample combines profiles from different redshifts (Paper I). We calculate a bootstrap uncertainty of the averaged profiles by randomly subsampling them 500 times.

When selecting individual haloes, we again apply a lower radial limit of r_min (Section 3.2) to avoid poorly resolved parts of the profiles, and we also omit bins that are expected to contain fewer than 5 particles based on the mean profile of the given sample (Paper I). Additionally, we only fit well-resolved haloes with N_200m ≥ 5000 because their larger resolved radial range leads to less noisy distributions of the best-fitting parameters. In total, we have fitted the profiles of about 378 000 haloes in the WMAP7, Planck, and self-similar cosmologies. The WMAP7 sample shrinks with redshift from about 45 000 at z = 0 to about 1300 at z = 6. Similarly, the self-similar samples shrink from about 120 000 haloes in the n = −1 simulation to about 13 000 in the n = −2.5 box.

3.3 Fitting procedure

Fitting density profiles is notoriously fickle, and the chosen routine can systematically influence the results (e.g. O’Neil et al. 2021). We have designed a robust routine, which we describe in detail in Appendix A1. We minimize the logarithmic residual between data and fitting functions using a Cauchy loss function, which severely reduces influence of outliers. For averaged profiles, the residual is compared to the bootstrap error plus a 5 per cent systematic uncertainty added to all bins. For individual haloes, we use the sum of the Poisson error due to (potentially) low number of particles and an arbitrary 25 per cent systematic error. Each fit is performed multiple times from a number of initial guesses to ensure that the true minimum is found (see Appendix A1 for details).

4 RESULTS

In Sections 4.1 and 4.2, we fit averaged and individual profiles, respectively, and analyse the quality of the fits. Given the large parameter space of mass, redshift, cosmology, and accretion rate, we show only representative examples and refer the reader to a collection of additional online figures (Section 1).

4.1 Fits to averaged profiles

Figs 3 and 4 show median and mean profiles for a wide variety of halo samples. We omit intermediate bins to avoid crowding the figures; their profiles fall between those shown. As in Paper I, all profile plots are split into total, orbiting, and infalling profiles (from top to bottom), with smaller panels showing the logarithmic slope and the fractional deviation of the fitting function from the profiles. We show the latter on a symmetric log scale, with the region between |$\pm 5{{\ \rm per\ cent}}$| being linear. The dotted lines highlight differences of |$5{{\ \rm per\ cent}}$|⁠, |$10{{\ \rm per\ cent}}$|⁠, and |$20{{\ \rm per\ cent}}$| to guide the eye.

$Fits to the selected median profiles of haloes binned by mass, redshift, accretion rate, and power spectrum slope (from left to right). Large panels show the total, orbiting, and infalling profiles (solid lines, from top to bottom), as well as fits (dashed lines). The smaller bottom panels show the logarithmic slope and the relative difference between simulation and fit, plotted on a symmetric log scale where the range between $\pm 5{{\ \rm per\ cent}}$ is linear; the dotted gray lines mark 5 per cent, 10 per cent, and 20 per cent. Shaded areas highlight the statistical uncertainty. The first column shows mass-selected samples in the WMAP7 cosmology. The second column shows the redshift evolution of a representative peak height bin (2 < ν < 3). In the third column, this bin is additionally split by accretion rate. The fourth column compares profiles at fixed ν and Γ in self-similar cosmologies. The total profiles are fit to 5 per cent accuracy, with the only exception of the wiggles in low-Γ profiles that are caused by particles on their second orbit. The orbiting profiles are fit to $\sim 10{{\ \rm per\ cent}}$ except where they sharply drop and even tiny differences in density lead to large relative differences. The infalling profiles are similarly fit to between 10 per cent and 20 per cent accuracy, depending on the halo sample.$

Figure 3.

Fits to the selected median profiles of haloes binned by mass, redshift, accretion rate, and power spectrum slope (from left to right). Large panels show the total, orbiting, and infalling profiles (solid lines, from top to bottom), as well as fits (dashed lines). The smaller bottom panels show the logarithmic slope and the relative difference between simulation and fit, plotted on a symmetric log scale where the range between |$\pm 5{{\ \rm per\ cent}}$| is linear; the dotted gray lines mark 5 per cent, 10 per cent, and 20 per cent. Shaded areas highlight the statistical uncertainty. The first column shows mass-selected samples in the WMAP7 cosmology. The second column shows the redshift evolution of a representative peak height bin (2 < ν < 3). In the third column, this bin is additionally split by accretion rate. The fourth column compares profiles at fixed ν and Γ in self-similar cosmologies. The total profiles are fit to 5 per cent accuracy, with the only exception of the wiggles in low-Γ profiles that are caused by particles on their second orbit. The orbiting profiles are fit to |$\sim 10{{\ \rm per\ cent}}$| except where they sharply drop and even tiny differences in density lead to large relative differences. The infalling profiles are similarly fit to between 10 per cent and 20 per cent accuracy, depending on the halo sample.

Open in new tab Download slide

$Same as Fig. 3 but for mean rather than median profiles. The total mean profiles are fit as well as their median counterparts. Similarly, the mean orbiting profiles are fit to about $10{{\ \rm per\ cent}}$ accuracy except at the truncation. The fit ignores bins where the density falls below ρ = 0.1 ρm; the corresponding deviations are driven by the influence of neighbouring haloes and thus not particularly meaningful. The mean infalling profiles are fit less accurately than the medians, largely because they show strong effects of the R200m boundary that is used to define subhaloes. However, this issue has no impact on the total profiles.$

Figure 4.

Same as Fig. 3 but for mean rather than median profiles. The total mean profiles are fit as well as their median counterparts. Similarly, the mean orbiting profiles are fit to about |$10{{\ \rm per\ cent}}$| accuracy except at the truncation. The fit ignores bins where the density falls below ρ = 0.1 ρ_m; the corresponding deviations are driven by the influence of neighbouring haloes and thus not particularly meaningful. The mean infalling profiles are fit less accurately than the medians, largely because they show strong effects of the R_200m boundary that is used to define subhaloes. However, this issue has no impact on the total profiles.

Open in new tab Download slide

In this section, we focus on fits with the full parameter freedom, that is five parameters for the orbiting term (Section 2.4) and three for the infalling term (Section 2.6). The orbiting and infalling profiles are fit separately with the respective models, and the results serve as initial conditions for a combined fit to the total profiles (in which we vary all parameters except for δ_max; Section 2.7). Thus, the total fits differ from the sum of the orbiting and infalling fits. In Paper III, we will calibrate some parameters to produce a more predictive profile model.

4.1.1 The total profile

We begin by evaluating fits to the median and mean total profiles, shown in the top rows of Figs 3 and 4. These fits do not make use of the separate orbiting and infalling components. In the left columns, we split halo samples from the WMAP7 cosmology by mass, or rather by peak height, ν. The corresponding mass bins range from M_200m ≈ 1.4 × 10¹⁰ M_⊙ (ν = 0.5) to M_200m > 8 × 10¹⁴M_⊙ (ν > 3). In the second column, we investigate the redshift evolution of the 2 < ν < 3 bin, which corresponds to 1.4 × 10¹⁴ < M_200m/M_⊙ < 8 × 10¹⁴ at z = 0 and to 7 × 10⁸ < M_200m/M_⊙ < 4.2 × 10¹⁰ at z = 6. The total profiles are fit excellently, to within 5 per cent or better. The slope panels demonstrate that the fits faithfully reproduce the complicated shape evolution of the profiles. Moreover, we do not notice any significant differences in the fit quality to the median and mean total profiles. The results are essentially the same for the Planck cosmology or for other bins in peak height.

DK14 and Paper I showed that the profile shapes depend more fundamentally on the mass accretion rate, Γ, than on mass, redshift, or cosmology. Many of the apparent trends with mass are actually trends in Γ, which are convolved with the mass-dependent distribution of accretion rates (larger haloes accrete more actively when structure forms hierarchically). In the third columns of Figs 3 and 4, we split the 2 < ν < 3 bin by accretion rates; the trends are the same for other peak heights (Paper I). An accretion rate of Γ ≲ 0.5 corresponds to pure pseudo-evolution due to the evolving overdensity threshold (Diemer, More & Kravtsov 2013), meaning that the sample shown in black has undergone essentially no actual changes to the density profile over the past dynamical time. Conversely, haloes with Γ > 4 (yellow) have grown rapidly. Fitting the resulting profiles is more challenging than for mass-selected samples because they exhibit a greater diversity of shapes. Nonetheless, our model fits the Γ-selected profiles almost as well as their ν-selected counterparts, namely to about 5 per cent within the statistical uncertainties in most samples and to about 10 per cent for the most extreme (low or high) accretion rates. One exception are the ‘wiggles’ near the second caustic in low-Γ profiles, where the fit can deviate by up to 20 per cent. These caustics form at the second orbit of recently accreted particles (Adhikari et al. 2014), but modelling this feature in our function would introduce significant complexity for a modest return in accuracy.

Finally, the right columns of Figs 3 and 4 show profiles from the self-similar simulations. We have selected haloes by both ν and Γ to isolate the effect of the slope of the linear power spectrum, n. Even though extreme values such as n = −1 lead to profiles that differ significantly from ΛCDM, the total profiles are fit to 5 per cent or better (subject to the trends with Γ discussed above). This match gives us confidence that our model can fit profiles in a wide range of cosmologies.

4.1.2 The orbiting term

We now turn to the orbiting profiles (rows 4–6 in Figs 3 and 4). The exact quality of the fit is somewhat academic because the orbiting and infalling terms cannot be exactly distinguished in observations, but it was a stated goal for our model to capture the physical nature of the orbiting term. At its truncation, however, the orbiting term becomes slightly sensitive to the numerical definition of particles’ pericentres (Paper I). Moreover, particles can be tidally stripped by interactions with nearby neighbours. While we try to suppress this effect by excluding haloes with a large unbound component (Section 3.2), a few ‘orbiting’ particles can be dragged to arbitrarily large radii. Given that this artefact is caused by relatively few haloes, the resulting differences are much more apparent in the mean than in the median profiles. These issues caution us not to overinterpret fitting errors at very low densities, which can reach ρ = 10⁻⁵ρ_m in some mean profiles. We thus ignore densities below ρ < 0.1 ρ_m in fits to the mean orbiting profiles.

The ν-selected orbiting profiles are fit to 10 per cent accuracy out to roughly the truncation radius, where the relative differences can become arbitrarily large. As for the total profiles, we find a somewhat degraded fit quality in Γ-selected samples, chiefly because the second-caustic wiggles are more pronounced in the orbiting than in the total profiles where the infalling term partially smooths them out. Overall, the fit quality is a little better for the median profiles than for the means.

4.1.3 The infalling term

Fitting the infalling term to small radii is an entirely new challenge because the shape of this contribution at r ≲ 0.5 R_200m was hitherto concealed by the orbiting component. Key insights from Paper I were that the infalling term approaches a constant value at small radii, continuously steepens with radius towards the truncation radius, and flattens again beyond r_t. However, the detailed profile shapes are diverse and depend on ν, Γ, and n, as evidenced by rows 7–9 of Figs 3 and 4.

Our model fits most median infalling profiles to about 10 per cent and some Γ-selected samples to 20 per cent accuracy, with no real trend with redshift, mass, or cosmology. Some of the deviations seen at small radii are not statistically significant considering the uncertainties (shaded areas). The mean infalling profiles are more difficult to fit and can deviate by more than 20 per cent for low-mass samples. One reason is that the mean profiles are noticeably affected by nearby neighbours, and thus by the definition of subhaloes as residing within R_200m (Paper I). This arbitrary boundary leads to sharp features around R_200m (e.g. blue lines in the left column of Fig. 4), which cannot (and should not) be fit by our simple power-law expression.

4.1.4 Summary

We find that the new truncated-exponential model provides 5–10 per cent fits to the total mean and median profiles selected by mass and/or accretion rate. Any features that are systemically fit poorly are understood to be somewhat unphysical (e.g. due to neighbouring haloes), too detailed for a fitting model (second caustic), or susceptible to noise (the exact shape at the orbiting truncation, where the density plummets to very low values). We note that baryonic effects enter at about the same level of accuracy, meaning that a more accurate fit would need to rely on hydrodynamical simulations. Fitting the orbiting and infalling terms separately poses a greater challenge, but the fitting function still captures their salient features and provides good accuracy over the radial ranges where the profiles are well constrained.

We have also experimented with fitting the averaged profiles of subhaloes. Here, we select by bound (as opposed to total) mass. The sparta algorithm keeps identifying first pericentres after a halo becomes a subhalo, but the orbiting term will increasingly include host material that the subhalo has drifted through. When binning by the resulting peak height, the truncated exponential model fits the median orbiting profiles reasonably well, as they do exhibit a clear truncation. The mean profiles, however, appear fairly irregular due to some subhaloes with strong host contributions. Similarly, the ‘infalling’ term now contains mostly host material and follows an entirely different shape from that of host haloes. Given these difficulties, we leave an investigation of subhalo profiles to future work.

4.2 Fits to individual haloes

One requirement for our fitting function was that it should fit the total profiles of individual haloes, as well as their dynamically split components. For each halo in our sample of 378 000, we separately fit equation (9) to the orbiting term (fixing α = 0.18 and β = 3) and equation (15) to the infalling term (Section 3.3). We then use the best-fitting parameters as the initial guess in a combined fit to the total profile. As discussed previously, individual halo profiles generally do not contain sufficient information to simultaneously constrain r_s and α, or to reliably determine r_t (necessitating other ways to measure the splashback radius, such as the sparta algorithm). Thus, we experiment with different levels of parameter freedom in the combined profiles. In one case, we let all parameters except δ_max float. This procedure naturally leads to the best fit in a statistical sense, but we observe numerous degeneracies and unphysical behaviour where the infalling profile adjusts to make up for a poor fit to the orbiting term. To avoid these issues, we also run fits where we vary only ρ_s and r_s in the combined fit, meaning that r_t, δ₁, δ_max, and s are fixed to their values from the separate fits. Given the restricted freedom of the combined fits, they are generally very similar to the sum of the separate fits. The goal of the restricted fits is not to find the formally most optimal fit but to obtain meaningful best-fitting parameters. We compare the fit quality of both types of fit to the NFW and Einasto models in Section 5.2. In Paper III, we show that the totally free fits do recover the parameters of the separate orbiting and infalling fits on average, albeit with large scatter.

Fig. 5 shows a few representative example fits. Each panel shows the orbiting, infalling, and total profiles of an individual halo as solid lines; the logarithmic slopes are omitted because they are noisy for individual haloes. Our model (dashed lines) fits the profiles well within statistical error and stochastic variations, except for the bottom-right-hand panel which was selected to show a poor fit. The top row shows haloes with increasing accretion rates. As discussed in Paper I, their peak height (mass) has very little influence on the profile shapes, but the increasing Γ clearly manifests itself in decreasing values of r_t/R_200m and a sharper truncation. The second row of Fig. 5 shows a number of examples from other redshifts and cosmologies. First, at high redshift, the average mass accretion rate is higher, which can lead to very small truncation radii such as in the left bottom panel, where r_t reaches the minimum imposed in our fits, 0.55 R_200m. However, we have not observed examples where the true r_t is so much smaller than this limit that it would lead to a bad fit. The high accretion rate also leads to a prominent infalling profile, with δ₁ ≈ 45 at R_200m; this example shows that our generous upper limit of 100 is necessary. The second and third panels demonstrate that our functions for the orbiting and infalling profiles can fit extreme cases, which are often found in the self-similar simulations with very different power-spectrum slopes. For example, the infalling profile in the second panel reaches the maximum slope of s = 4, and infalling matter contributes equally to orbiting matter all the way to r ≈ 0.2 R_200m. Conversely, the third panel shows a halo with flat infalling and orbiting profiles, which lead to r_s reaching the maximum of 0.45 R_200m. Finally, the bottom right-hand panel was selected to show a poor fit, in this case because of a major merger that deposited unusually large amount of orbiting material at large radii.

Figure 5.

Representative examples of fits to individual halo profiles. Solid lines show the orbiting, infalling, and total profiles (in dark-blue, purple, and light-blue, respectively). Shaded contours highlight the combination of statistical and systematic error that is used in the fits (equation A3). Thin lines without error contours indicate unreliable profile bins smaller than the resolution limit. The dashed dark-blue and purple lines show fits of our models for the orbiting and infalling terms; the dashed light-blue lines correspond to the sum of these fits rather than the combined fit. The dotted gray lines show a three-parameter Einasto fit to the orbiting term for comparison. The legends in each panel list properties of the halo in question and the most informative best-fitting parameters. The top row shows haloes from the WMAP7 cosmology at z = 0 with ascending Γ, which leads to a decreasing r_t/R_200m. The normalization, slope, and asymptotic values of the infalling profiles differ dramatically between individual haloes. The bottom row shows a few extreme cases. First, a halo at z = 4 has a very high accretion rate, Γ ≈ 8, which leads to a truncation radius at the lower limit of r_t ≥ 0.55 R_200m as well as a high normalization of the infalling profile. The second example is a halo in the n = −1 self-similar cosmology, which often produces very steep infalling profiles that can reach our maximum of s = 4. Conversely, the n = −2.5 cosmology commonly produces shallow, power-law like total and infalling profiles with low concentrations (third example). The final example was deliberately picked to represent a poor fit, in this case because of a large amount of orbiting material at large radii. While an Einasto profile with the same number of free parameters (α instead of r_t) can fit low-Γ profiles, it cannot capture the sharp truncation (dotted lines).

Open in new tab Download slide

We conclude that the new truncated-exponential model captures individual profiles well, including the full range of halo properties and cosmologies investigated. In Section 5.2, we furthermore show that the model outperforms the Einasto function without a truncation. In Paper III, we analyse the resulting parameter distributions in detail.

5 COMPARISON TO OTHER MODELS

We have convinced ourselves that the new fitting function successfully reproduces simulated profiles, but we should also ask how unique that success is and whether it could be achieved with fewer free parameters. In Section 5.1, we compare models for the orbiting term based on averaged profiles, where we can discern the detailed shapes that lead to the success or failure of certain fitting functions. In Section 5.2, we evaluate the fit quality for individual haloes by comparing to the NFW and Einasto forms. In Section 5.3, we compare the surface density profiles and lensing signal to other models and show that the features of the three-dimensional profiles do persist in projection. In Section 5.4, we consider models based on the distribution function and discuss the dynamical properties of the new model.

5.1 Models for the orbiting term

Fig. 6 shows a comparison of different fits to the total and orbiting profiles of two representative halo samples. The infalling part of the total profile is fit with equation (15) in all cases and is thus omitted from the figure. The slope panels for the orbiting term reach down to extremely steep slopes to highlight the predictions for the asymptotic shape of this component that are made by different fitting functions (without the benefit of knowing the split profiles, of course).

$Comparison of different fitting functions for the orbiting term. The profiles in rows 1 and 4 are rescaled by r2 to compress the dynamical range. The following rows show the logarithmic slope and fractional differences, as in Figs 3 and 4. The columns show two representative orbiting profiles (thick black lines), namely the median profiles of haloes with 2 < ν < 3 (left) and the mean profiles of a sub-sample additionally selected to have 2 < Γ < 4 (right). The statistical uncertainty on the average profiles is smaller than the thickness of the black lines. The gray lines show NFW (dashed) and Einasto (dot–dashed) profiles, which are not designed to fit the sharp truncation of the orbiting term; those profiles are omitted from the difference panels for clarity. The coloured lines show the profiles of Baltz et al. (2009, BMO, dashed maroon), DK14 with the slopes β and γ as free parameters (yellow dot–dashed), DK14 with the slopes fixed to their suggested values for ν- and Γ-selected samples (orange dot–dashed), and the new model (light blue). While the full-freedom DK14 fits well, it has an extra free parameter compared to the new function. Moreover, it is clear that its slope asymptotically approaches a certain log-linear evolution with radius due to the underlying Einasto profile. This shape fundamentally differs from the observed truncation, which is reproduced by the new model. For this figure, we fit to the entire radial range of the mean profiles, including possibly unphysical features at low density that are overfit by the free DK14 model (right column).$

Figure 6.

Comparison of different fitting functions for the orbiting term. The profiles in rows 1 and 4 are rescaled by r² to compress the dynamical range. The following rows show the logarithmic slope and fractional differences, as in Figs 3 and 4. The columns show two representative orbiting profiles (thick black lines), namely the median profiles of haloes with 2 < ν < 3 (left) and the mean profiles of a sub-sample additionally selected to have 2 < Γ < 4 (right). The statistical uncertainty on the average profiles is smaller than the thickness of the black lines. The gray lines show NFW (dashed) and Einasto (dot–dashed) profiles, which are not designed to fit the sharp truncation of the orbiting term; those profiles are omitted from the difference panels for clarity. The coloured lines show the profiles of Baltz et al. (2009, BMO, dashed maroon), DK14 with the slopes β and γ as free parameters (yellow dot–dashed), DK14 with the slopes fixed to their suggested values for ν- and Γ-selected samples (orange dot–dashed), and the new model (light blue). While the full-freedom DK14 fits well, it has an extra free parameter compared to the new function. Moreover, it is clear that its slope asymptotically approaches a certain log-linear evolution with radius due to the underlying Einasto profile. This shape fundamentally differs from the observed truncation, which is reproduced by the new model. For this figure, we fit to the entire radial range of the mean profiles, including possibly unphysical features at low density that are overfit by the free DK14 model (right column).

Open in new tab Download slide

Most previously proposed fitting functions encounter at least one of two fundamental issues. First, Paper I showed that the slope of the orbiting term reaches arbitrarily steep values (as low as −13 in Fig. 6). This kind of cutoff cannot be captured by power laws because they inevitably asymptote to a fixed slope at large radii. For example, the NFW profile approaches a slope of −3 (dashed gray lines), too shallow for the simulated profiles at r ≳ r_t. This issue persists for profiles with other slope transitions such as Burkert (1995) or superNFW (Lilley, Evans & Sanders 2018), models with steeper outer slopes (Hernquist 1990), or even more complicated combinations of power-law slopes such as generalized power-law models (Zhao 1996; An & Zhao 2013; Di Cintio et al. 2014; Dekel et al. 2017; Freundlich et al. 2020). To improve the fit of NFW profiles around the transition region, Baltz et al. (2009, BMO) introduced a multiplicative steepening, |$\rho \propto \rho _{\rm NFW} \times [r_{\rm t}^2/(r^2 + r_{\rm t}^2)]^n$| (maroon dashed lines in Fig. 6). This form can fit the total profiles to about 25 per cent accuracy, but it fails to capture the orbiting term because it approaches a fixed slope of −3 − 2n. Here n ≈ 11, a value much larger than the n ≈ 2 usually employed to fit weak lensing signals (Oguri & Hamana 2011). An approach similar to BMO was taken by Tavio et al. (2008); we do not include their function in Fig. 6 because it was already shown not to reach sufficiently steep slopes (fig. 15 in DK14).

The second fundamental problem encountered by many fitting functions is that the orbiting profiles have two physical scales, which can be expressed as r_s and r_t. These scales are set by different epochs in a halo’s accretion history, namely, by the formation time, which correlates with c, and r_t, which correlates with the recent accretion history (Wechsler et al. 2002; Tasitsiomi et al. 2004; Zhao et al. 2009; Ludlow et al. 2013; Lucie-Smith, Adhikari & Wechsler 2022b; Shin & Diemer 2022). The need for an extra variable to describe the profiles can be independently discovered using machine learning (Lucie-Smith et al. 2022a). Single-scale functions cannot address this problem and typically smooth over the sharp truncation, though they may work if the fit does not extend beyond roughly R_vir. The Einasto model exemplifies this behaviour: while its slope can reach arbitrarily steep values in principle, the steepening occurs at the same pace at all radii (gray dot–dashed lines in Fig. 6). Some models have introduced a second radial scale to account for variations at small scales, e.g. in the coreNFW (Read, Agertz & Collins 2016) or coreEinasto (Lazar et al. 2020) models, but these modifications obviously do not change the fit at large radii. Springel & White (1999) suggested ρ ∝ exp [ − (r − R_200c)/r_s] to model tidal truncation, but this term does not introduce a new scale and thus cannot capture variations in r_t/r_s (see also Fielder et al. 2020).

To account for the second radial scale and the strong steepening, DK14 introduced a flexible truncation term (Section 2.3). Their results were based on the same simulations as this paper, but they did not split profiles into orbiting and infalling, meaning that their functional form (equation 5) cannot be expected to reproduce the truncation shape at low densities. Fig. 6 shows two versions of the DK14 model. In the first, the sharpness of the transition, β, and the asymptotic slope of the steepening term, γ, are free parameters (hereafter ‘DK14-6’, yellow dot–dashed lines). In the second, they are fixed to the optimal values recommended by DK14, namely β = 4 and γ = 8 for the ν-selected sample in the left column and β = 6 and γ = 4 for Γ-selected sample on the right (hereafter ‘DK14-4’, orange dot–dashed lines). We will keep in mind that the former version has one more free parameter than the new model, and the latter one fewer.

Fig. 6 demonstrates that the DK14 profile approaches unphysical slopes at large radii: regardless of the power-law slope γ, the exponential term from the Einasto profile eventually takes over and leads to a gradually steepening slope that does not match the simulated profiles. While this transition occurs at steep slopes (e.g. about −10 in the left column of Fig. 6) and does not cause noticeable fit errors, it does lead to a conceptually wrong prediction for the shape of the orbiting term. For the purposes of this test, we have included even radii with low-mean densities in the fit (unlike in our fiducial procedure, see Section 3.2).

We can get a quantitative sense of the new and DK14 fits by comparing their χ²/N_dof values. The conclusions depend somewhat on whether samples are selected by ν only or also Γ, mean and median, and whether we consider the DK14-6 or DK14-4 variants. We first consider fits to only the orbiting term, which test whether the fitting function physically describes the correct underlying profile. As expected, the new form outperforms DK14-4 for virtually all mean and median samples, and often by a sizeable margin, with Δ_χ ≡ Δlog₁₀(χ²/N_dof) of up to 2. Compared to the DK14-6 form that has one more free parameter, the median ν-selected profiles tend to be better fit by the new form, and the mean profiles by DK14-6. Most Γ-selected samples are better fit by DK14-6, but at the expense of an extreme range of best-fitting β and γ. The new model better fits particularly sharp truncations in median samples with Γ ≳ 2. All differences in fit quality are strongly reduced when fitting the total rather than the orbiting profiles, but DK14-4 still tends to perform slightly worse and DK14-6 slightly better and in most fits. Some of the latter model’s success, however, arises from overfitting profiles that carry signatures of neighbouring haloes, which leave a power-law tail that the DK14 fit latches onto. Physically, it is not clear that the profile model should fit such profiles.

In summary, we confirm the expected trend that models with more free parameters achieve a lower χ²/N_dof. However, giving the DK14 model its full freedom of six parameters often leads to meaningless values of β and γ. Our new truncated-exponential model achieves more physical asymptotic profiles and a comparable fit quality with fewer well-defined parameters.

5.2 Individual profiles

While the examples of individual haloes in Fig. 5 are encouraging, we need a quantitative metric to decide whether the additional complexity of the new model is warranted compared to simpler functions. Without an absolute measure of fit quality in the absence of a well-motivated uncertainty on the profiles (Appendix A1), our best option is to compare the relative χ²/N_dof of different fits. To this end, we repeat the procedure described in Section 4.2 with three models: NFW, Einasto with free α, and Einasto with α = 0.18. In Fig. 7, we compare these fits for all 378 000 haloes in our individual sample. While there are mild trends when we break up the sample by cosmology and redshift, the overall picture remains unchanged.

$Comparison of the relative fit quality of equation (9) and the NFW and Einasto models, for all 378 000 haloes in our individual halo sample. The functions have been applied to the orbiting term (left) and to the total profile (right). All models including the same fit to the infalling profile. The histograms show the logarithmic ratio of χ2/Ndof between each model and our new model (with three free parameters, ρs, rs, and rt). The points with error bars mark the mean ratio and 68 per cent interval (the median ratios are very similar). When fitting the orbiting term, the new model strongly outperforms the three-parameter Einasto profile (green), highlighting that freedom in α is not sufficient to fit the truncation. When fitting the total profiles but keeping α and rt fixed based on the separate fit (right-hand panel), the fit quality is improved by 10 per cent on average (and a median of $5{{\ \rm per\ cent}}$). This result holds even when fitting only to the total term with no knowledge of the split profiles (yellow). Fixing α in the Einasto fits (light blue) degrades the orbiting fit but also reduces the number of degrees of freedom so that the χ2/Ndof difference in the total fits remains similar. The two-parameter NFW model compares less well than the Einasto profile (dark blue).$

Figure 7.

Comparison of the relative fit quality of equation (9) and the NFW and Einasto models, for all 378 000 haloes in our individual halo sample. The functions have been applied to the orbiting term (left) and to the total profile (right). All models including the same fit to the infalling profile. The histograms show the logarithmic ratio of χ²/N_dof between each model and our new model (with three free parameters, ρ_s, r_s, and r_t). The points with error bars mark the mean ratio and 68 per cent interval (the median ratios are very similar). When fitting the orbiting term, the new model strongly outperforms the three-parameter Einasto profile (green), highlighting that freedom in α is not sufficient to fit the truncation. When fitting the total profiles but keeping α and r_t fixed based on the separate fit (right-hand panel), the fit quality is improved by 10 per cent on average (and a median of |$5{{\ \rm per\ cent}}$|⁠). This result holds even when fitting only to the total term with no knowledge of the split profiles (yellow). Fixing α in the Einasto fits (light blue) degrades the orbiting fit but also reduces the number of degrees of freedom so that the χ²/N_dof difference in the total fits remains similar. The two-parameter NFW model compares less well than the Einasto profile (dark blue).

Open in new tab Download slide

The most ‘fair’ comparison is between the new model and an Einasto profile with free α because both have three parameters. When fitting only the orbiting profile, the new model has a lower χ²/N_dof by a median factor of four (green histogram in the left-hand panel of Fig. 7). This dramatic improvement highlights that the truncation cannot generally be fit with Einasto profiles. This observation is also apparent in Fig. 5, where low-Γ haloes are well described by the Einasto form (gray dotted lines) but where the increasingly sharp truncation leads to poor fits at large radii.

In realistic applications, however, we fit to total profiles. This case is shown in the right-hand panel of Fig. 7. We test two variations of the fits. First, we vary only ρ_s and r_s in the combined fit, keeping α, r_t, and the parameters of the infalling profile to their values from the separate fits (green). The new model now improves χ²/N_dof by about 10 per cent on average because exact shape of the truncation is less important in total profiles. Moreover, in most real-world applications, we do not have access to the separate orbiting profile. This result holds when we ignore the results of the separate orbiting fit and let all free parameters of the orbiting term vary (yellow). This relatively modest difference reminds us that it is difficult to extract the truncation radius from total individual profiles, but it still represents a measurable improvement over the Einasto model. We also compare the new model to the two-parameter Einasto fit with α = 0.18 and to the NFW profile (light and dark blue histograms in Fig. 7). The orbiting profile is described much less accurately by these forms, with mean χ²/N_dof increases of about 7 and 12, respectively. Once again, the difference is much smaller when fitting the entire profile, with mean increases of 17 per cent and 34 per cent. We have confirmed that these results remain unchanged when directly optimizing χ² instead of the Cauchy loss function ln (1 + χ²) (Appendix A1).

Overall, Fig. 7 provides strong justification for our new fitting function because it improves the fit quality at a fixed number of free parameters. While this improvement is most noticeable when fitting the orbiting term separately, we also record a statistically significant improvement when fitting the total profiles without any knowledge of the orbiting term. We conclude that the new function captures both averaged and individual profiles well. In Paper III, we analyse the connection between the properties of individual haloes and their best-fitting parameters.

5.3 Projected profiles and lensing signal

While we have demonstrated that the new profile model is preferred by simulated halo profiles, the three-dimensional density profile is generally not accessible in observations. In general, we measure either the surface density or the lensing signal. The surface (or projected) density, Σ(r), can be measured via the satellite distribution, for example (More et al. 2016). The lensing signal corresponds to the excess surface mass density, |$\Delta \Sigma \equiv \overline{\Sigma }(\lt r) - \Sigma (r)$|⁠, the difference between the averaged-surface density within a given radius and the surface density at that radius. In this section, we investigate how different the new profile is from other models in projection, and whether variations in the truncation radius can, in principle, be extracted from observational data.

In Fig. 8, we show projected profile quantities for a typical cluster halo with |$M_{\rm 200m}= 5 \times 10^{14}\ h^{-1}\, \mathrm{M}_\odot$| and c_200m = 5 at z = 0. We set α = 0.18 and β = 0.3, but we vary r_t between values corresponding to low- and high-accretion rates (about 0.7 < r_t/R_200m < 1.2, Paper III). The surface density profiles of these models can easily be distinguished by eye, and all truncated exponential profiles visibly differ from the NFW and Einasto forms. The differences are smaller in the lensing signal, but they persist across a wide radial range. This finding highlights that the entire profile must be fitted in order to extract the full information content (e.g. Xhakaj et al. 2020, 2022).

$The projected (surface) density profile (top) and lensing signal (excess surface mass density, bottom) for characteristic profiles of a cluster halo with $M_{\rm 200m}= 5 \times 10^{14}\, h^{-1}\, \mathrm{M}_\odot$ and c200m = 5. We compare the NFW, Einasto, and truncated exponential profiles, where we force the truncation radius of the latter to the values listed in the legend and renormalize the profile to the same total mass. The differences in the transition region clearly propagate to the surface density profiles, where the different values of rt are easily distinguishable by eye. The differences are smaller but propagate to a much larger radial range in ΔΣ because Σ(r) is subtracted from the averaged surface density.$

Figure 8.

The projected (surface) density profile (top) and lensing signal (excess surface mass density, bottom) for characteristic profiles of a cluster halo with |$M_{\rm 200m}= 5 \times 10^{14}\, h^{-1}\, \mathrm{M}_\odot$| and c_200m = 5. We compare the NFW, Einasto, and truncated exponential profiles, where we force the truncation radius of the latter to the values listed in the legend and renormalize the profile to the same total mass. The differences in the transition region clearly propagate to the surface density profiles, where the different values of r_t are easily distinguishable by eye. The differences are smaller but propagate to a much larger radial range in ΔΣ because Σ(r) is subtracted from the averaged surface density.

Open in new tab Download slide

5.4 Dynamical models

Assuming that haloes are in dynamical equilibrium and that the particle velocities are isotropic (which is not exactly realistic, e.g. Hansen & Moore 2006), a model for the density profile implies a particular form of the velocity dispersion, σ², and the phase–space distribution function, f(ε). The mathematical expressions for these quantities are given in Appendix D. The velocity dispersion can be measured in principle, e.g. from the motions of satellite galaxies or stars (Mamon, Biviano & Murante 2010; Okumura et al. 2018; Adhikari et al. 2019; Hamabata, Oguri & Nishimichi 2019; Tomooka et al. 2020; Bose & Loeb 2021; Aung et al. 2021, 2022). We do not expect the dynamics to be fundamentally altered by the presence of baryons (Callingham et al. 2020).

In this section, we check whether the new model makes reasonable predictions for the dynamical quantities and we question whether other energy-based models could make similar predictions. The former is not automatically guaranteed since not all profiles have positive continuous distribution functions (Baes & Camps 2021). We define ψ(r) to be the (positive) gravitational potential of the halo and ε ≡ ψ(r) − v²/2 to be the relative binding energy of particles. Particles with ε < 0 would be unbound and thus lead to f(ε < 0) = 0. Conversely, the maximum binding energy would correspond to particles that reside at the centre of the halo with velocity v = 0, meaning that ε = ψ(0).

Fig. 9 shows radial profiles of density, enclosed mass, potential, and velocity dispersion, as well as the distribution function. We compare our new model (dark blue) to the NFW and Einasto forms (gray). The profiles are normalized to have the same total mass at the largest radius shown, and they are expressed in dimensionless units rescaled by R, some scale density ρ₀, and the gravitational constant G. The results thus depend only on the relative scales r_s/R and r_t/R, as well as the steepening parameters α and β. We choose a representative set of parameters for all models, namely r_s/R = 0.2, α = 0.18, r_t/R = 1, and β = 4. We observe that M(r) flattens as a consequence of the truncation in ρ, and that the potential approaches ψ ∝ 1/r at a smaller radius than for the Einasto profile.

Figure 9.

Dynamical properties of our new model (dark blue), a King (1966) profile (light blue, with a comparable cutoff radius), as well as the NFW and Einasto forms (gray). The panels show the density ρ, enclosed mass M, potential ψ, isotropic velocity dispersion σ², and distribution function f(ε) (from top to bottom). The latter falls to zero for unbound particles with ε < 0; conversely, ε = ψ(0) corresponds to maximally bound particles at the centre of the halo. The profiles are expressed in dimensionless units and normalized to reach the same total mass M = 1 at the right edge of the plot. Compared to the conventional models, the sharp density cutoff in the new model leads to a mass distribution that suddenly flattens, a rapidly dropping velocity dispersion, and a sharp drop in the distribution function at the energy corresponding to orbits that reach their apocentre near the truncation radius.

Open in new tab Download slide

The density cutoff is also mirrored in energy space, where the velocity dispersion falls sharply. The corresponding distribution function is similar to the NFW and Einasto profile at high-binding energies, but it drops steeply at the binding energy corresponding to particles at the truncation radius. This is the desired behaviour for a profile that attempts to model the edge of the orbital distribution (e.g. King 1966; Drakos, Taylor & Benson 2017; Amorisco 2021). In other words, the splashback radius (the apocentre of the most recently accreted particles) is governed by a soft limit in the distribution of particle kinetic energies. Conversely, the NFW and Einasto profiles have small, but finite, support out to the smallest binding energies (see also e.g. Cardone et al. 2005; Mamon & Łokas 2005; Beraldo e Silva et al. 2014; Baes & Ciotti 2019). For large values of β such as 10 (i.e. very sharp truncation in density space), the distribution function can become non-monotonic near its cutoff.

These observations raise the questions of whether the truncation of the orbiting term could be described by a truncation of f(ε) at some maximum energy. One well-known example of such models is the King (1966) family of profiles, which corresponds to an isothermal sphere that is sharply cut off at a certain binding energy (see also Michie 1963). This model is defined by the distribution function

$$\begin{eqnarray*} f_{\rm king}(\varepsilon) = f_0 \left(e^{(\varepsilon - \varepsilon _{\rm t}) / \sigma ^2} - 1 \right) \, , \end{eqnarray*}$$

(16)

where f₀ is a normalization, σ the velocity dispersion of the isothermal sphere, and ε_t the energy where the object is truncated (f → 0). We use the dimensionless variables of Drakos et al. (2017) to numerically integrate equation (16) to find the potential, density, and velocity dispersion profiles. Given our fundamental scale R and the mass normalization, the profile has only one free parameter Z_t ≡ ε_t/[σ²ψ(0)]. The light-blue lines in Fig. 9 show a King profile with Z_t = 0.2, a cutoff energy chosen to roughly match the real-space truncation radius of r_t = 0.2 R. The position and sharpness of the truncation can be adjusted using Z_t (e.g. fig. 4.8 in Binney & Tremaine 2008), but Fig. 9 immediately explains why the King model cannot describe realistic halo profiles: the density approaches a large fixed-density core at r ≪ R. Moreover, the truncation in binding energy leads to a sharper real-space truncation than we observe in simulations. Our new model allows for lower binding energies (bottom panel) and thus for a smoother truncation.

Of course, the King model is only one example of dynamics-based profiles. Another recently suggested model is DARKexp (Hjorth & Williams 2010), but this model cannot fit the truncation because its slope approaches γ = −4 at large radii, regardless of the ϕ₀ parameter that controls the distribution of particle energies (Williams & Hjorth 2010; Williams, Hjorth & Wojtak 2010). Similarly, Pontzen & Governato (2013) derive a distribution function by maximizing the entropy, but it is not clear whether this model produces a truncation at the right radius (their fig. 7; see also Wagner 2020).

6 CONCLUSIONS

We have presented a new fitting function for halo-density profiles, which is composed of models for the orbiting (one-halo) and infalling terms. The orbiting term is modelled as an extension of the Einasto profile, ρ_orb ∝ exp [ − 2/α (r/r_s)^α − 1/β (r/r_t)^β]. The infalling term is a power law in density that smoothly approaches a maximum value at small radii. These functions are implemented in the publicly available colossus code, including numerical routines for their mass, surface density, and ΔΣ integrals. Our main conclusions are as follows.

The new model fits mean and median total profiles to roughly |$\pm 5{{\ \rm per\ cent}}$| across a vast range of radius, halo mass, redshift, and cosmology. When haloes are also selected by accretion rate, the fit quality can degrade to |$\pm 20{{\ \rm per\ cent}}$| near the transition between the orbiting and infalling profiles.
The new model accurately fits the orbiting and infalling terms separately, even down to densities well below the cosmic mean. The exact fit quality depends on the chosen mass, accretion rate, and whether mean or median profiles are considered.
By fixing α = 0.18 and β = 3, the orbiting model becomes a 3-parameter fit that captures the profiles of individual haloes more accurately than Einasto profiles (on average).
The sharp truncation of the orbiting term corresponds to a truncation in the binding energy of particles, which is not replicated in more extended profile models.
Different truncation radii lead to clear differences in the projected surface density profiles and in the lensing signal ΔΣ.
We introduce an augmented ‘Model B’ that fixes the logarithmic slope at the scale radius to be −2. This formulation is almost equivalent to the fiducial model but alleviates some rare parameter degeneracies at the cost of slightly increased complexity.

In Paper III, we will analyse the best-fitting parameters of averaged and individual halo profiles. We have left a number of theoretical and practical questions unexplored. For example, we intend to provide analytical approximations for the mass and projected density of the new fitting function. Another urgent question is the relationship between the profile parameters and definitions of the halo boundary based on the orbiting population of subhaloes (Aung et al. 2021; García et al. 2021), or other definitions based on the interplay between the infalling and orbiting components (Fong & Han 2021). Finally, we intend to apply the new model to observational data.

ACKNOWLEDGEMENTS

I am grateful to Nicole Drakos for help with computing King profiles and to Maarten Baes for verifying the dynamical properties of the new model. I thank Andrew Hearin and Keiichi Umetsu for comments on a draft, and Susmita Adhikari, Han Aung, Barun Dhar, Rafael Garcia, Daisuke Nagai, Eduardo Rozo, and Angus Wright 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).

DATA AVAILABILITY

The sparta code that was used to extract the dynamically split density profiles from our simulations 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 on the author’s website 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.

REFERENCES

Adhikari

S.

,

Dalal

N.

,

Chamberlain

R. T.

,

2014

,

JCAP

,

11

,

19

10.1088/1475-7516/2014/11/019

Crossref

Search ADS

Adhikari

S.

,

Dalal

N.

,

More

S.

,

Wetzel

A.

,

2019

,

ApJ

,

878

,

9

10.3847/1538-4357/ab1a39

Crossref

Search ADS

Adhikari

S.

et al. ,

2021

,

ApJ

,

923

,

37

10.3847/1538-4357/ac0bbc

Crossref

Search ADS

Amorisco

N. C.

,

2021

,

preprint

(

arXiv

)

An

J.

,

Zhao

H.

,

2013

,

MNRAS

,

428

,

2805

10.1093/mnras/sts175

Crossref

Search ADS

Aung

H.

,

Nagai

D.

,

Rozo

E.

,

García

R.

,

2021

,

MNRAS

,

502

,

1041

10.1093/mnras/staa3994

Crossref

Search ADS

Aung

H.

,

Nagai

D.

,

Rozo

E.

,

Wolfe

B.

,

Adhikari

S.

,

2022

,

preprint

(

arXiv

)

Baes

M.

,

Camps

P.

,

2021

,

MNRAS

,

503

,

2955

10.1093/mnras/stab634

Crossref

Search ADS

Baes

M.

,

Ciotti

L.

,

2019

,

A&A

,

626

,

A110

10.1051/0004-6361/201935320

Crossref

Search ADS

Baltz

E. A.

,

Marshall

P.

,

Oguri

M.

,

2009

,

J. Cosmology Astropart. Phys.

,

1

,

15

10.1088/1475-7516/2009/01/015

Crossref

Search ADS

Baxter

E.

et al. ,

2017

,

ApJ

,

841

,

18

10.3847/1538-4357/aa6ff0

Crossref

Search ADS

Becker

M. R.

,

Kravtsov

A. V.

,

2011

,

ApJ

,

740

,

25

10.1088/0004-637X/740/1/25

Crossref

Search ADS

Behroozi

P. S.

,

Wechsler

R. H.

,

Wu

H.-Y.

,

2013a

,

ApJ

,

762

,

109

10.1088/0004-637X/762/2/109

Crossref

Search ADS

Behroozi

P. S.

,

Wechsler

R. H.

,

Wu

H.-Y.

,

Busha

M. T.

,

Klypin

A. A.

,

Primack

J. R.

,

2013b

,

ApJ

,

763

,

18

10.1088/0004-637X/763/1/18

Crossref

Search ADS

Beraldo e Silva

L.

,

Lima

M.

,

Sodré

L.

,

Perez

J.

,

2014

,

Phys. Rev. D

,

90

,

123004

10.1103/PhysRevD.90.123004

Crossref

Search ADS

Bertschinger

E.

,

1985

,

ApJS

,

58

,

39

10.1086/191028

Crossref

Search ADS

Betancort-Rijo

J. E.

,

Sanchez-Conde

M. A.

,

Prada

F.

,

Patiri

S. G.

,

2006

,

ApJ

,

649

,

579

10.1086/506311

Crossref

Search ADS

Bianconi

M.

,

Buscicchio

R.

,

Smith

G. P.

,

McGee

S. L.

,

Haines

C. P.

,

Finoguenov

A.

,

Babul

A.

,

2021

,

ApJ

,

911

,

136

10.3847/1538-4357/abebd7

Crossref

Search ADS

Binney

J.

,

1982

,

MNRAS

,

200

,

951

10.1093/mnras/200.4.951

Crossref

Search ADS

Binney

J.

,

Tremaine

S.

,

2008

,

Galactic Dynamics

, 2nd edn.

Princeton Univ. Press

,

Princeton

Bond

J. R.

,

Kofman

L.

,

Pogosyan

D.

,

1996

,

Nature

,

380

,

603

10.1038/380603a0

Crossref

Search ADS

Bose

S.

,

Loeb

A.

,

2021

,

ApJ

,

912

,

114

10.3847/1538-4357/abec77

Crossref

Search ADS

Brown

S. T.

,

McCarthy

I. G.

,

Diemer

B.

,

Font

A. S.

,

Stafford

S. G.

,

Pfiefer

S.

,

2020

,

MNRAS

,

495

,

4994

10.1093/mnras/staa1491

Crossref

Search ADS

Burkert

A.

,

1995

,

ApJ

,

447

,

L25

10.1086/309560

Crossref

Search ADS

Callingham

T. M.

,

Cautun

M.

,

Deason

A. J.

,

Frenk

C. S.

,

Grand

R. J. J.

,

Marinacci

F.

,

Pakmor

R.

,

2020

,

MNRAS

,

495

,

12

10.1093/mnras/staa1089

Crossref

Search ADS

Cardone

V. F.

,

Piedipalumbo

E.

,

Tortora

C.

,

2005

,

MNRAS

,

358

,

1325

10.1111/j.1365-2966.2005.08834.x

Crossref

Search ADS

Chang

C.

et al. ,

2018

,

ApJ

,

864

,

83

10.3847/1538-4357/aad5e7

Crossref

Search ADS

Cooray

A.

,

Sheth

R.

,

2002

,

Phys. Rep.

,

372

,

1

10.1016/S0370-1573(02)00276-4

Crossref

Search ADS

Courteau

S.

et al. ,

2014

,

Rev. Mod. Phys.

,

86

,

47

10.1103/RevModPhys.86.47

Crossref

Search ADS

Crocce

M.

,

Pueblas

S.

,

Scoccimarro

R.

,

2006

,

MNRAS

,

373

,

369

10.1111/j.1365-2966.2006.11040.x

Crossref

Search ADS

Dacunha

T.

,

Belyakov

M.

,

Adhikari

S.

,

Shin

T.-h.

,

Goldstein

S.

,

Jain

B.

,

2022

,

MNRAS

,

512

,

4378

10.1093/mnras/stac392

Crossref

Search ADS

de Blok

W. J. G.

,

2010

,

Adv. Astron.

,

2010

,

789293

10.1155/2010/789293

Crossref

Search ADS

Dehnen

W.

,

1993

,

MNRAS

,

265

,

250

10.1093/mnras/265.1.250

Crossref

Search ADS

Dekel

A.

,

Ishai

G.

,

Dutton

A. A.

,

Maccio

A. V.

,

2017

,

MNRAS

,

468

,

1005

10.1093/mnras/stx486

Crossref

Search ADS

Dhar

B. K.

,

2021

,

MNRAS

,

504

,

4583

10.1093/mnras/stab1029

Crossref

Search ADS

Dhar

B. K.

,

Williams

L. L. R.

,

2010

,

MNRAS

,

405

,

340

10.1111/j.1365-2966.2010.16446.x

Crossref

Search ADS

Di Cintio

A.

,

Brook

C. B.

,

Macciò

A. V.

,

Stinson

G. S.

,

Knebe

A.

,

Dutton

A. A.

,

Wadsley

J.

,

2014

,

MNRAS

,

437

,

415

10.1093/mnras/stt1891

Crossref

Search ADS

Diemand

J.

,

Kuhlen

M.

,

2008

,

ApJ

,

680

,

L25

10.1086/589688

Crossref

Search ADS

Diemer

B.

,

2017

,

ApJS

,

231

,

5

10.3847/1538-4365/aa799c

Crossref

Search ADS

Diemer

B.

,

2018

,

ApJS

,

239

,

35

10.3847/1538-4365/aaee8c

Crossref

Search ADS

Diemer

B.

,

2020

,

ApJS

,

251

,

17

10.3847/1538-4365/abbf51

Crossref

Search ADS

Diemer

B.

,

2022

,

MNRAS

,

513

,

573

(Paper I)

10.1093/mnras/stac878

Crossref

Search ADS

Diemer

B.

,

Kravtsov

A. V.

,

2014

,

ApJ

,

789

,

1

(DK14)

10.1088/0004-637X/789/1/1

Crossref

Search ADS

Diemer

B.

,

Kravtsov

A. V.

,

2015

,

ApJ

,

799

,

108

10.1088/0004-637X/799/1/108

Crossref

Search ADS

Diemer

B.

,

More

S.

,

Kravtsov

A. V.

,

2013

,

ApJ

,

766

,

25

10.1088/0004-637X/766/1/25

Crossref

Search ADS

Dooley

G. A.

,

Griffen

B. F.

,

Zukin

P.

,

Ji

A. P.

,

Vogelsberger

M.

,

Hernquist

L. E.

,

Frebel

A.

,

2014

,

ApJ

,

786

,

50

10.1088/0004-637X/786/1/50

Crossref

Search ADS

Doroshkevich

A. G.

,

Shandarin

S. F.

,

1978

,

Soviet Ast.

,

22

,

653

Drakos

N. E.

,

Taylor

J. E.

,

Benson

A. J.

,

2017

,

MNRAS

,

468

,

2345

10.1093/mnras/stx652

Crossref

Search ADS

Dubinski

J.

,

Carlberg

R. G.

,

1991

,

ApJ

,

378

,

496

10.1086/170451

Crossref

Search ADS

Dutton

A. A.

,

Macciò

A. V.

,

2014

,

MNRAS

,

441

,

3359

10.1093/mnras/stu742

Crossref

Search ADS

Eckert

D.

,

Ettori

S.

,

Robertson

A.

,

Massey

R.

,

Pointecouteau

E.

,

Harvey

D.

,

McCarthy

I. G.

,

2022

,

A&A

,

666

,

A41

10.1051/0004-6361/202243205

Crossref

Search ADS

Eddington

A. S.

,

1916

,

MNRAS

,

76

,

572

10.1093/mnras/76.7.572

Crossref

Search ADS

Efstathiou

G.

,

Frenk

C. S.

,

White

S. D. M.

,

Davis

M.

,

1988

,

MNRAS

,

235

,

715

10.1093/mnras/235.3.715

Crossref

Search ADS

Einasto

J.

,

1965

,

Tr. Astrofiz. Inst. Alma-Ata

,

5

,

87

Einasto

J.

,

1969

,

Astrophys.

,

5

,

67

10.1007/BF01013353

Crossref

Search ADS

Erwin

P.

,

2015

,

ApJ

,

799

,

226

10.1088/0004-637X/799/2/226

Crossref

Search ADS

Fielder

C. E.

,

Mao

Y.-Y.

,

Zentner

A. R.

,

Newman

J. A.

,

Wu

H.-Y.

,

Wechsler

R. H.

,

2020

,

MNRAS

,

499

,

2426

10.1093/mnras/staa2851

Crossref

Search ADS

Fong

M.

,

Han

J.

,

2021

,

MNRAS

,

503

,

4250

10.1093/mnras/stab259

Crossref

Search ADS

Freundlich

J.

et al. ,

2020

,

MNRAS

,

499

,

2912

10.1093/mnras/staa2790

Crossref

Search ADS

Fukushige

T.

,

Makino

J.

,

2001

,

ApJ

,

557

,

533

10.1086/321666

Crossref

Search ADS

Gao

L.

,

Navarro

J. F.

,

Cole

S.

,

Frenk

C. S.

,

White

S. D. M.

,

Springel

V.

,

Jenkins

A.

,

Neto

A. F.

,

2008

,

MNRAS

,

387

,

536

10.1111/j.1365-2966.2008.13277.x

Crossref

Search ADS

García

R.

,

Rozo

E.

,

Becker

M. R.

,

More

S.

,

2021

,

MNRAS

,

505

,

1195

10.1093/mnras/stab1317

Crossref

Search ADS

García

R.

,

Salazar

E.

,

Rozo

E.

,

Adhikari

S.

,

Aung

H.

,

Diemer

B.

,

Nagai

D.

,

Wolfe

B.

,

2022

,

preprint

(

arXiv

)

Genina

A.

et al. ,

2018

,

MNRAS

,

474

,

1398

10.1093/mnras/stx2855

Crossref

Search ADS

Ghigna

S.

,

Moore

B.

,

Governato

F.

,

Lake

G.

,

Quinn

T.

,

Stadel

J.

,

2000

,

ApJ

,

544

,

616

10.1086/317221

Crossref

Search ADS

Graham

A. W.

,

Merritt

D.

,

Moore

B.

,

Diemand

J.

,

Terzić

B.

,

2006

,

AJ

,

132

,

2701

10.1086/508990

Crossref

Search ADS

Gunn

J. E.

,

Gott III

J. R.

,

1972

,

ApJ

,

176

,

1

10.1086/151605

Crossref

Search ADS

Hamabata

A.

,

Oguri

M.

,

Nishimichi

T.

,

2019

,

MNRAS

,

489

,

1344

10.1093/mnras/stz2227

Crossref

Search ADS

Hansen

S. H.

,

Moore

B.

,

2006

,

New A

,

11

,

333

10.1016/j.newast.2005.09.001

Crossref

Search ADS

Harris

C. R.

et al. ,

2020

,

Nature

,

585

,

357

10.1038/s41586-020-2649-2

Crossref

Search ADS

PubMed

Hernquist

L.

,

1990

,

ApJ

,

356

,

359

10.1086/168845

Crossref

Search ADS

Hjorth

J.

,

Williams

L. L. R.

,

2010

,

ApJ

,

722

,

851

10.1088/0004-637X/722/1/851

Crossref

Search ADS

Hunter

J. D.

,

2007

,

Comput. Sci. Eng.

,

9

,

90

10.1109/MCSE.2007.55

Crossref

Search ADS

Jaffe

W.

,

1983

,

MNRAS

,

202

,

995

10.1093/mnras/202.4.995

Crossref

Search ADS

King

I. R.

,

1966

,

AJ

,

71

,

64

10.1086/109857

Crossref

Search ADS

Klypin

A. A.

,

Trujillo-Gomez

S.

,

Primack

J.

,

2011

,

ApJ

,

740

,

102

10.1088/0004-637X/740/2/102

Crossref

Search ADS

Klypin

A.

,

Yepes

G.

,

Gottlöber

S.

,

Prada

F.

,

Heß

S.

,

2016

,

MNRAS

,

457

,

4340

10.1093/mnras/stw248

Crossref

Search ADS

Knollmann

S. R.

,

Power

C.

,

Knebe

A.

,

2008

,

MNRAS

,

385

,

545

10.1111/j.1365-2966.2008.12857.x

Crossref

Search ADS

Komatsu

E.

et al. ,

2011

,

ApJS

,

192

,

18

10.1088/0067-0049/192/2/18

Crossref

Search ADS

Lazar

A.

et al. ,

2020

,

MNRAS

,

497

,

2393

10.1093/mnras/staa2101

Crossref

Search ADS

Lewis

A.

,

Challinor

A.

,

Lasenby

A.

,

2000

,

ApJ

,

538

,

473

10.1086/309179

Crossref

Search ADS

Lilley

E. J.

,

Evans

N. W.

,

Sanders

J. L.

,

2018

,

MNRAS

,

476

,

2086

10.1093/mnras/sty295

Crossref

Search ADS

Łokas

E. L.

,

Mamon

G. A.

,

2001

,

MNRAS

,

321

,

155

10.1046/j.1365-8711.2001.04007.x

Crossref

Search ADS

Lucie-Smith

L.

,

Peiris

H. V.

,

Pontzen

A.

,

Nord

B.

,

Thiyagalingam

J.

,

Piras

D.

,

2022a

,

Phys. Rev. D

,

105

,

103533

10.1103/PhysRevD.105.103533

Crossref

Search ADS

Lucie-Smith

L.

,

Adhikari

S.

,

Wechsler

R. H.

,

2022b

,

MNRAS

,

515

,

2164

10.1093/mnras/stac1833

Crossref

Search ADS

Ludlow

A. D.

,

Angulo

R. E.

,

2017

,

MNRAS

,

465

,

L84

10.1093/mnrasl/slw216

Crossref

Search ADS

Ludlow

A. D.

,

Navarro

J. F.

,

White

S. D. M.

,

Boylan-Kolchin

M.

,

Springel

V.

,

Jenkins

A.

,

Frenk

C. S.

,

2011

,

MNRAS

,

415

,

3895

10.1111/j.1365-2966.2011.19008.x

Crossref

Search ADS

Ludlow

A. D.

et al. ,

2013

,

MNRAS

,

432

,

1103

10.1093/mnras/stt526

Crossref

Search ADS

Ludlow

A. D.

,

Bose

S.

,

Angulo

R. E.

,

Wang

L.

,

Hellwing

W. A.

,

Navarro

J. F.

,

Cole

S.

,

Frenk

C. S.

,

2016

,

MNRAS

,

460

,

1214

10.1093/mnras/stw1046

Crossref

Search ADS

Ludlow

A. D.

,

Schaye

J.

,

Bower

R.

,

2019

,

MNRAS

,

488

,

3663

10.1093/mnras/stz1821

Crossref

Search ADS

Mamon

G. A.

,

Łokas

E. L.

,

2005

,

MNRAS

,

362

,

95

10.1111/j.1365-2966.2005.09225.x

Crossref

Search ADS

Mamon

G. A.

,

Biviano

A.

,

Murante

G.

,

2010

,

A&A

,

520

,

A30

10.1051/0004-6361/200913948

Crossref

Search ADS

Mansfield

P.

,

Avestruz

C.

,

2021

,

MNRAS

,

500

,

3309

10.1093/mnras/staa3388

Crossref

Search ADS

Merritt

D.

,

Graham

A. W.

,

Moore

B.

,

Diemand

J.

,

Terzić

B.

,

2006

,

AJ

,

132

,

2685

10.1086/508988

Crossref

Search ADS

Michie

R. W.

,

1963

,

MNRAS

,

125

,

127

10.1093/mnras/125.2.127

Crossref

Search ADS

Moore

B.

,

Quinn

T.

,

Governato

F.

,

Stadel

J.

,

Lake

G.

,

1999

,

MNRAS

,

310

,

1147

10.1046/j.1365-8711.1999.03039.x

Crossref

Search ADS

More

S.

,

Diemer

B.

,

Kravtsov

A. V.

,

2015

,

ApJ

,

810

,

36

10.1088/0004-637X/810/1/36

Crossref

Search ADS

More

S.

et al. ,

2016

,

ApJ

,

825

,

39

10.3847/0004-637X/825/1/39

Crossref

Search ADS

Murata

R.

,

Sunayama

T.

,

Oguri

M.

,

More

S.

,

Nishizawa

A. J.

,

Nishimichi

T.

,

Osato

K.

,

2020

,

PASJ

,

72

,

64

10.1093/pasj/psaa041

Crossref

Search ADS

Navarro

J. F.

,

Frenk

C. S.

,

White

S. D. M.

,

1995

,

MNRAS

,

275

,

720

Crossref

Search ADS

Navarro

J. F.

,

Frenk

C. S.

,

White

S. D. M.

,

1996

,

ApJ

,

462

,

563

10.1086/177173

Crossref

Search ADS

Navarro

J. F.

,

Frenk

C. S.

,

White

S. D. M.

,

1997

,

ApJ

,

490

,

493

10.1086/304888

Crossref

Search ADS

Navarro

J. F.

et al. ,

2004

,

MNRAS

,

349

,

1039

10.1111/j.1365-2966.2004.07586.x

Crossref

Search ADS

Navarro

J. F.

et al. ,

2010

,

MNRAS

,

402

,

21

10.1111/j.1365-2966.2009.15878.x

Crossref

Search ADS

O’Donnell

C.

,

Behroozi

P.

,

More

S.

,

2021

,

MNRAS

,

509

,

3285

10.1093/mnras/stab3170

Crossref

Search ADS

O’Neil

S.

,

Barnes

D. J.

,

Vogelsberger

M.

,

Diemer

B.

,

2021

,

MNRAS

,

504

,

4649

10.1093/mnras/stab1221

Crossref

Search ADS

O’Neil

S.

,

Borrow

J.

,

Vogelsberger

M.

,

Diemer

B.

,

2022

,

MNRAS

,

513

,

835

10.1093/mnras/stac850

Crossref

Search ADS

Oguri

M.

,

Hamana

T.

,

2011

,

MNRAS

,

414

,

1851

10.1111/j.1365-2966.2011.18481.x

Crossref

Search ADS

Okumura

T.

,

Nishimichi

T.

,

Umetsu

K.

,

Osato

K.

,

2018

,

Phys. Rev. D

,

98

,

023523

10.1103/PhysRevD.98.023523

Crossref

Search ADS

Peng

C. Y.

,

Ho

L. C.

,

Impey

C. D.

,

Rix

H.-W.

,

2010

,

AJ

,

139

,

2097

10.1088/0004-6256/139/6/2097

Crossref

Search ADS

Planck Collaboration XVI

,

2014

,

A&A

,

571

,

A16

10.1051/0004-6361/201321591

Crossref

Search ADS

Pontzen

A.

,

Governato

F.

,

2013

,

MNRAS

,

430

,

121

10.1093/mnras/sts529

Crossref

Search ADS

Power

C.

,

Navarro

J. F.

,

Jenkins

A.

,

Frenk

C. S.

,

White

S. D. M.

,

Springel

V.

,

Stadel

J.

,

Quinn

T.

,

2003

,

MNRAS

,

338

,

14

10.1046/j.1365-8711.2003.05925.x

Crossref

Search ADS

Prada

F.

,

Klypin

A. A.

,

Simonneau

E.

,

Betancort-Rijo

J.

,

Patiri

S.

,

Gottlöber

S.

,

Sanchez-Conde

M. A.

,

2006

,

ApJ

,

645

,

1001

10.1086/504456

Crossref

Search ADS

Read

J. I.

,

Agertz

O.

,

Collins

M. L. M.

,

2016

,

MNRAS

,

459

,

2573

10.1093/mnras/stw713

Crossref

Search ADS

Retana-Montenegro

E.

,

van Hese

E.

,

Gentile

G.

,

Baes

M.

,

Frutos-Alfaro

F.

,

2012

,

A&A

,

540

,

A70

10.1051/0004-6361/201118543

Crossref

Search ADS

Ricotti

M.

,

Pontzen

A.

,

Viel

M.

,

2007

,

ApJ

,

663

,

L53

10.1086/520113

Crossref

Search ADS

Schneider

A.

,

Teyssier

R.

,

Stadel

J.

,

Chisari

N. E.

,

Le Brun

A. M. C.

,

Amara

A.

,

Refregier

A.

,

2019

,

J. Cosmology Astropart. Phys.

,

2019

,

020

10.1088/1475-7516/2019/03/020

Crossref

Search ADS

Shapiro

P. R.

,

Iliev

I. T.

,

Raga

A. C.

,

1999

,

MNRAS

,

307

,

203

10.1046/j.1365-8711.1999.02609.x

Crossref

Search ADS

Shin

T.

et al. ,

2019

,

MNRAS

,

487

,

2900

10.1093/mnras/stz1434

Crossref

Search ADS

Shin

T.

et al. ,

2021

,

MNRAS

,

507

,

5758

10.1093/mnras/stab2505

Crossref

Search ADS

Shin

T.-h.

,

Diemer

B.

,

2022

,

preprint

(

arXiv

)

Springel

V.

,

2005

,

MNRAS

,

364

,

1105

10.1111/j.1365-2966.2005.09655.x

Crossref

Search ADS

Springel

V.

,

White

S. D. M.

,

1999

,

MNRAS

,

307

,

162

10.1046/j.1365-8711.1999.02613.x

Crossref

Search ADS

Stadel

J.

,

Potter

D.

,

Moore

B.

,

Diemand

J.

,

Madau

P.

,

Zemp

M.

,

Kuhlen

M.

,

Quilis

V.

,

2009

,

MNRAS

,

398

,

L21

10.1111/j.1745-3933.2009.00699.x

Crossref

Search ADS

Sugiura

H.

,

Nishimichi

T.

,

Rasera

Y.

,

Taruya

A.

,

2020

,

MNRAS

,

493

,

2765

10.1093/mnras/staa413

Crossref

Search ADS

Tasitsiomi

A.

,

Kravtsov

A. V.

,

Gottlöber

S.

,

Klypin

A. A.

,

2004

,

ApJ

,

607

,

125

10.1086/383219

Crossref

Search ADS

Tavio

H.

,

Cuesta

A. J.

,

Prada

F.

,

Klypin

A. A.

,

Sanchez-Conde

M. A.

,

2008

,

preprint

(

arXiv

)

Teyssier

R.

,

Pontzen

A.

,

Dubois

Y.

,

Read

J. I.

,

2013

,

MNRAS

,

429

,

3068

10.1093/mnras/sts563

Crossref

Search ADS

Tomooka

P.

,

Rozo

E.

,

Wagoner

E. L.

,

Aung

H.

,

Nagai

D.

,

Safonova

S.

,

2020

,

MNRAS

,

499

,

1291

10.1093/mnras/staa2841

Crossref

Search ADS

Udrescu

S. M.

,

Dutton

A. A.

,

Macciò

A. V.

,

Buck

T.

,

2019

,

MNRAS

,

482

,

5259

10.1093/mnras/sty3112

Crossref

Search ADS

Umetsu

K.

,

2020

,

A&A Rev.

,

28

,

7

10.1007/s00159-020-00129-w

Crossref

Search ADS

Umetsu

K.

,

Diemer

B.

,

2017

,

ApJ

,

836

,

231

10.3847/1538-4357/aa5c90

Crossref

Search ADS

Velliscig

M.

,

van Daalen

M. P.

,

Schaye

J.

,

McCarthy

I. G.

,

Cacciato

M.

,

Le Brun

A. M. C.

,

Vecchia

C. D.

,

2014

,

MNRAS

,

442

,

2641

10.1093/mnras/stu1044

Crossref

Search ADS

Virtanen

P.

et al. ,

2020

,

Nat. Methods

,

17

,

261

10.1038/s41592-019-0686-2

Crossref

Search ADS

PubMed

Wagner

J.

,

2020

,

Gen. Relativ. Gravit.

,

52

,

61

10.1007/s10714-020-02715-w

Crossref

Search ADS

Wang

J.

,

Bose

S.

,

Frenk

C. S.

,

Gao

L.

,

Jenkins

A.

,

Springel

V.

,

White

S. D. M.

,

2020

,

Nature

,

585

,

39

10.1038/s41586-020-2642-9

Crossref

Search ADS

PubMed

Wechsler

R. H.

,

Tinker

J. L.

,

2018

,

ARA&A

,

56

,

435

10.1146/annurev-astro-081817-051756

Crossref

Search ADS

Wechsler

R. H.

,

Bullock

J. S.

,

Primack

J. R.

,

Kravtsov

A. V.

,

Dekel

A.

,

2002

,

ApJ

,

568

,

52

10.1086/338765

Crossref

Search ADS

Williams

L. L. R.

,

Hjorth

J.

,

2010

,

ApJ

,

722

,

856

10.1088/0004-637X/722/1/856

Crossref

Search ADS

Williams

L. L. R.

,

Hjorth

J.

,

Wojtak

R.

,

2010

,

ApJ

,

725

,

282

10.1088/0004-637X/725/1/282

Crossref

Search ADS

Xhakaj

E.

,

Diemer

B.

,

Leauthaud

A.

,

Wasserman

A.

,

Huang

S.

,

Luo

Y.

,

Adhikari

S.

,

Singh

S.

,

2020

,

MNRAS

,

499

,

3534

10.1093/mnras/staa3046

Crossref

Search ADS

Xhakaj

E.

,

Leauthaud

A.

,

Lange

J.

,

Hearin

A.

,

Diemer

B.

,

Dalal

N.

,

2022

,

MNRAS

,

514

,

2876

10.1093/mnras/stac941

Crossref

Search ADS

Zhao

H.

,

1996

,

MNRAS

,

278

,

488

10.1093/mnras/278.2.488

Crossref

Search ADS

Zhao

D. H.

,

Jing

Y. P.

,

Mo

H. J.

,

Börner

G.

,

2009

,

ApJ

,

707

,

354

10.1088/0004-637X/707/1/354

Crossref

Search ADS

APPENDIX A: DETAILS ON FITTING PROCEDURE

In this Appendix, we give additional detail on the technical aspects of our profile fits. In Section A1, we discuss our fitting procedure and in Section A2, we describe the reasoning behind and meaning of the parameter limits given in Section 2.7 and Table 1.

A1 Fitting procedure

A first important choice is the definition of the residual. The conventional approach of minimizing the square of (ρ_fit − ρ_data)/σ works poorly in the outer profiles because the summed χ² is dominated by the much higher densities at small radii. Some authors have chosen to minimize r²(ρ_fit − ρ_data)/σ instead, but this function still represents an arbitrary weighting with radius. One can even fit the slope profile, but the results depend on how the slope is estimated from noisy data (O’Neil et al. 2021). To obtain a fit of equal quality at all radii, we minimize the square logarithmic differences between simulated and fitted profiles. We find that we obtain the most physically sensible fits if we aggressively reduce the effect of outliers by using a Cauchy loss function, so that the minimized quantity is

$$\begin{eqnarray*} \chi _{\rm cauchy}^2 \equiv \sum _i \ln \left(1 + \chi _{{\rm i}}^2 \right) \quad \mathrm{with} \quad \chi _{{\rm i}} \equiv \frac{\ln \rho _{\rm i,fit} - \ln \rho _{\rm i,data}}{\ln (1 + \sigma _{\rm i}/ \rho _{\rm i,data})} \, , \end{eqnarray*}$$

(A1)

where ρ_{i, fit} is the fit density in bin i, ρ_{i, data} is the density of the simulated profile (individual or averaged), and σ_i is the uncertainty in ρ_{i, data}. Throughout the paper, we quote the usual values of |$\chi ^2 = \sum \chi _{\rm i}^2$|⁠, i.e. the logarithmic differences without taking the loss function into account. The additional logarithm in the Cauchy loss function allows the fitter to ignore features such as wiggles in favour of a better fit to the well-defined inner profile.

The next decision is how to define the uncertainty σ_i of each profile bin. The bootstrap estimate, σ_bs, captures statistical variations in averaged profiles, but it does not account for systematic errors that could arise from numerics, from our profile splitting algorithm, from baryonic effects, and (perhaps most importantly) from substructure. We crudely account for such errors by adding a systematic uncertainty of 5 per cent in quadrature,

$$\begin{eqnarray*} \sigma _{\rm i,mean/median} = \sqrt{\sigma _{\rm i,bs}^2 + \sigma _{\rm sys,av}^2} = \sqrt{\sigma _{\rm i,bs}^2 + (0.05\ \rho _{\rm i,data})^2} \, . \end{eqnarray*}$$

(A2)

The systematic error does have some effect on on the best-fitting profiles and parameters. It essentially represents a trade-off: if its value is small, the fit ‘trusts’ the statistically well-constrained inner profile over the noisier transition region; if its value is large, all bins are weighted more or less equally. Given that the tests in Paper I demonstrated convergence to roughly 5 per cent accuracy, this value is a reasonable approximation to the systematic uncertainty.

We minimize χ² using the ‘trust region reflective’ (trf) algorithm of the least-squares function in scipy. We find this algorithm to converge more reliably than the similar ‘dogbox’ variant. Both methods are faster than a conventional Levenberg–Marquardt solver and allow us to impose the parameter limits listed in Table 1. All parameters are fit in log space to avoid zero or negative values. Despite these optimizations, any steepest-gradient solver is liable to ‘fall into’ the first local minimum in χ² it encounters, which may be far from the global minimum. We avoid false minima with the following procedure. We begin from a fiducial initial guess and run the fitter with a relatively low relative target accuracy of 10⁻³. Using the output as a starting point, we vary the initial guess of each parameter to its extremes, that is, fractions of 0.05 and 0.95 between its lower and upper bounds in log space. If the resulting fit improves on the previous χ², we set the new results as the initial guess. If the total gradient at the solution (‘optimality’) is larger than 10⁻³, we repeat the procedure for more initial guesses of the parameter in question, namely fractions of [0.2, 0.35, 0.5, 0.65, 0.8] of its allowed log interval. Once all parameters have been varied, we accept the initial guess that resulted in the lowest χ² of all attempts and rerun the fit with a higher accuracy of 10⁻⁵ to obtain the final result. At modest cost, this algorithm ensures that virtually all fits converge to a reasonable solution and that the results are independent of the initial guess. The values for χ²/N_dof when fitting the averaged orbiting term hover around unity, though with a large spread that depends on the sample selection and whether we fit mean or median. When fitting the total profiles, χ²/N_dof roughly ranges between 0.01 and 1.

Fitting individual halo profiles can be particularly challenging because substructure and halo-to-halo scatter can lead to significant deviations from any fitting function (e.g. Umetsu & Diemer 2017). We thus slightly adjust the procedure laid out in the previous section for individual halo fits, as well as the parameter ranges (Table 1). Since the simulations do not provide a meaningful estimate of the uncertainty on individual profiles, we assume a Poisson-like statistical uncertainty due to the number of particles as a minimum,

$$\begin{eqnarray*} \sigma _{\rm i,ind} = \rho _{\rm i,data} \sqrt{\frac{1}{N_{\rm i,ptl}} + \sigma _{\rm sys,ind}^2} \, , \end{eqnarray*}$$

(A3)

where N_{i, ptl} is the number of particles in the given bin and σ_{sys, ind} is a fixed systematic error. Since most bins contain N_{i, ptl} ≫ 1 particles, the normalization of χ² is largely determined by the somewhat arbitrary value for σ_{sys, ind}. This systematic error term effectively weighs bins with low particle number against well-resolved bins because all bins with |$N_{\rm i,ptl} \gg 1 / \sigma _{\rm sys,ind}^2$| receive roughly the same weight in the fit. This method prevents the fit from being dominated by a few bins with high N_ptl. We set σ_{sys, ind} = 0.25, meaning that the systematic error dominates for bins with N_{i, ptl} > 16. With this value, the median χ²/N_dof ranges from 0.5 to 2 for the different halo samples. We have verified that the choice of systematic error does not greatly influence the overall distribution of the best-fitting parameters.

A2 Parameter Limits

The normalization, ρ_s, is tightly related to the formation redshift of haloes (Navarro, Frenk & White 1997; Ludlow et al. 2013) and thus cannot vary arbitrarily. None the less, we allow a generous range because ρ_s is well constrained in all fits. If we instead use the density at the centre, ρ₀, by setting C = 0 in Section 2.4, some parameter combinations can extrapolate to extreme central densities and the limits need to be much more flexible (Table 1).

Some profiles steepen so slowly that they exhibit no true scale radius, i.e. γ ≈ −2 over a wide range of radii. In such cases, the fit tends to shift r_s towards large values that can even exceed r_t, which is technically a valid solution because of the symmetry between the r_s and r_t terms in equation (9). After extensive experimentation, we find it necessary to enforce the correct, physical ordering of the terms by explicitly requiring r_t > r_s. In particular, we set 0.01 < r_s/R_200m < 0.45 and 0.5 < r_t/R_200m < 3. The range of r_s limits the concentration to 2.2 < c_200m < 100, but haloes outside of this range would correspond to extreme outliers in the population (e.g. Diemer & Kravtsov 2015); we never observe a preference for r_s ≳ 0.4 R_200m in well-constrained fits. The lower limit corresponds roughly to the lowest radii where our averaged profiles are resolved, meaning that r_s cannot be measured reliably near 0.01 R_200m. This limit could be revisited with high-resolution simulations that probe smaller radii. Similarly, we do not observe a truncation at radii smaller than 0.5 R_200m in any profiles, and we generally find little evidence for r_t > 2 R_200m. If r_t pushes against the upper limit, the profile does not experience a clear truncation, meaning that the values of r_t and β are arbitrary. We have experimented with falling back on an Einasto fit in such cases, but the fit quality is almost always degraded, which demonstrates that almost all orbiting profiles experience some steepening beyond that of the Einasto slope.

For the steepening parameter, we allow 0.03 < α < 0.4, the most extreme values that could reasonably be ascribed to any profiles we have tested. The lower limit corresponds to an almost invariant slope near −2, which is the case for some haloes with high accretion rates (Paper I). The Einasto profile can suffer from a well-known degeneracy between r_s and α if profiles do not extend to small radii or if their slope evolves slowly with radius (e.g. Ricotti, Pontzen & Viel 2007; Udrescu et al. 2019). While virtually all averaged profiles do constrain α to some extent, most individual halo profiles do not. In particular, the average fit quality remains roughly constant within the range 0.1 < α < 0.2, highlighting that individual profiles simply do not contain sufficient information to determine α. We could set it based on an ν–α relation (Gao et al. 2008; Klypin et al. 2016), but these trends with peak height turn out to be an artefact of the Einasto form poorly fitting the true, truncated profiles (Paper III). Instead, we find that α ≈ 0.18 in most averaged profiles, and that it depends on the accretion rate (Paper III) and on the slope of the linear power spectrum (Paper I; Ludlow & Angulo 2017; Brown et al. 2020). While the profiles of individual haloes with very high Γ do indeed show a mild preference for low α, the differences in fit quality are modest at best. Thus, we simply fix α = 0.18 in all fits to individual profiles. This value leads to concentrations that are similar to those measured using NFW profiles (Dutton & Macciò 2014; Ludlow et al. 2016), but other choices have been made in the literature (e.g. α = 0.16; Wang et al. 2020).

We have no intuitive prior for the transition sharpness β, but the profiles shown in Paper I range from gradual to very sharp transitions. Thus, we allow β to fluctuate between 0.1 and 10. While some fits technically prefer even larger values, those correspond to an essentially instantaneous cutoff, where β is not well defined. We might worry that the r_t-β term suffers from a degeneracy similar to that of the r_s-α term, but we find that r_t and β are not intrinsically degenerate and both well-constrained for most averaged profiles (Paper III). In individual haloes, however, the orbiting profiles commonly become unresolved at the radii that most constrain β. We thus fix β = 3, which barely decreases the average fit quality of the total profiles. To accommodate profiles that do not exhibit any noticeable truncation, we allow the truncation radius to shift to very large values where it does not matter, r_t/R_200m ≤ 10.

In the infalling profile, we limit the normalization to 1 < δ₁ < 100. Added to the mean density of the Universe, the lower limit corresponds to ρ_inf(R_200m) = 2ρ_m. Such low overdensities are observed in the self-similar simulation with n = −1. We leave the slope to fluctuate between 0.01 < s < 4, which includes any reasonable (positive) value encountered in simulated profiles. The asymptotic central overdensity δ_max reaches values between 10 and 1000 in the profiles where it is constrained, but we allow 1 < δ_max < 10⁴ to allow for profiles that do not exhibit a measurable maximum. We apply the limits discussed so far in the separate fits to the orbiting and infalling profile components. When we combine them to fit the total profiles, we do not vary δ_max because the infalling profile is subdominant to the orbiting term at r < R_200m, meaning that δ_max is unconstrained. In a situation where only the total profile is known, δ_max can be set to a reasonable value (e.g. 500) without a significant impact on the fit.

APPENDIX B: MODEL VARIANT WITH CORRECTION AT SCALE RADIUS

In Section 2.4, we noted that the truncation term in equation (9) has the undesirable effect of breaking the γ(r_s) = −2 condition. We now construct a model variant that maintains this condition, which we call Model B (in contrast to the fiducial model, which we shall call Model A for brevity). We wish to introduce a new term in the slope that enforces γ(r_s) = −2, but we need to be careful. For example, if we add (r_s/r_t)^β into the slope of Model A, γ(r) is still integrable and adds the term (r_s/r_t)^βln (r) to S(r) from equation (7). While this new model would now satisfy γ(r_s) = −2, γ has been renormalized at all radii with a fixed term that survives all the way to the centre of the halo, γ(0) = (r_s/r_t)^β. This feature is highly undesirable because it ties the central slope to the entirely unrelated properties of the truncation term. Instead, we introduce a multiplicative term that vanishes at r = 0 and approaches (r_s/r_t)^β at r = r_s,

$$\begin{eqnarray*} \gamma (r) = -2 \left(\frac{r}{r_{\rm s}} \right)^\alpha - \left(\frac{r}{r_{\rm t}} \right)^\beta + \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta \, . \end{eqnarray*}$$

(B1)

We integrate γ/r to find

$$\begin{eqnarray*} S(r) = -\frac{2}{\alpha } \left(\frac{r}{r_{\rm s}} \right)^\alpha -\frac{1}{\beta } \left(\frac{r}{r_{\rm t}} \right)^\beta + \frac{1}{\eta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta + C \, . \end{eqnarray*}$$

(B2)

Once again, we can set C = 0 or ensure S(r_s) = 0 by setting

$$\begin{eqnarray*} C = \frac{2}{\alpha } + \left[ \frac{1}{\beta } - \frac{1}{\eta } \right] \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \, , \end{eqnarray*}$$

(B3)

which causes the familiar extra terms from equations (4) and (7),

$$\begin{eqnarray*} S &=& -\frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha - 1 \right] -\frac{1}{\beta } \left[ \left(\frac{r}{r_{\rm t}} \right)^\beta - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \right]\nonumber \\ && + \frac{1}{\eta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[ \left(\frac{r}{r_{\rm s}} \right)^\eta - 1 \right] \, . \end{eqnarray*}$$

(B4)

While this form looks somewhat complicated, it has a straightforward interpretation: its slope is renormalized to match −2 at the scale radius but this normalization factor decays towards r = 0. Here, η is a nuisance parameter that determines how fast this decay happens: the larger η, the faster the slope approaches zero at small radii. On the other hand, η also introduces an undesirable change in the slope at larger radii (where r/r_s ≳ r_s/r_t). However, as long as β ≫ 1 and η < 1, this correction is barely noticeable. We fix η = 0.1, but the Model B fits are insensitive in the range 0.01 < η < 0.2. In some fits to mean orbiting-only profiles, the best-fitting value of α depends slightly on η, but this dependence disappears when fitting the entire profile or median orbiting profiles.

Fig. A1 shows the impact of different parameter values on Model B and highlights differences to Model A as dashed lines. As expected, Model B produces profiles that are indistinguishable from Model A for all but the most extreme parameter values. Differences are apparent for small r_t and/or small β. The example of r_t/R_200m = 0.3 is outside of the allowed parameter range for our fits (r_t ≥ 0.5R_200m), but we do allow very small values of β. Model A fixes the slope at r_t to be γ(r_t) = −2(r_t/r_s)^α − 1, which lies between −3 and −4 for most parameter values. Conversely, Model B fixes the slope at r_s to be γ = −2 regardless of the other parameters, but the slope at r_t becomes γ(r_t) = −2(r_t/r_s)^α − 1 + (r_t/r_s)^{η − β}. The β-dependent correction term is negligible except for the smallest β.

Figure A1.

Same as Fig. 1, but showing the modified Model B with η = 0.1 (solid lines) and a comparison to our fiducial Model A [equation (9), dashed lines]. Wherever the dashed lines are invisible, the models make indistinguishable predictions. The most important differences occur in the right two columns. While Model B fixes the slope at γ(r_s) = −2 regardless of the values of r_t and β, Model A can exhibit slight deviations from this value if r_t and/or β are small.

Open in new tab Download slide

The advantage of Model B is that r_s takes on its usual meaning as the radius where γ = −2, which avoids parameter degeneracies in power law-like profiles where the scale radius is poorly constrained (Paper I). We have, however, verified that models A and B are indistinguishable for virtually all mean and median profile fits. Wherever there are small differences, Model B tends to fit slightly better. The most significant differences in the best-fitting parameters occur for r_s in low-Γ averaged samples and in some individual halo profiles. The values of r_t and β are almost never affected, except in ill-defined fits. In summary, Model B is, fundamentally, the same profile model as Model A, and the best-fitting parameters do not depend on the model except in some pathological cases. While we have exclusively shown Model A fits throughout this paper, we will also investigate best-fitting parameters from Model B in Paper III.

APPENDIX C: DERIVATIVES WITH RESPECT TO PARAMETERS

Least-squares fits are often faster and more reliable if the user provides derivatives of the fitted function with respect to its free parameters θ. For orbiting profiles that follow the form ρ = ρ_se^S(r), those derivatives become

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln \theta } = \theta \frac{{\rm d} S}{{\rm d} \theta } \, . \end{eqnarray*}$$

(C1)

The derivative with respect to the normalization is dln ρ/dln ρ_s = 1 regardless of S(r), and we thus do not repeat it for the models listed below. We have double-checked the following expressions numerically and against Wolfram Alpha.

C1 Orbiting: Einasto

If we set C = 0, the derivatives with respect to the parameters are

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm s}} & =& 2 \left(\frac{r}{r_{\rm s}} \right)^\alpha \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \alpha } & =& \frac{2}{\alpha } \left(\frac{r}{r_{\rm s}} \right)^\alpha \left[1 - \alpha \ln \left(\frac{r}{r_{\rm s}} \right)\right] \, . \end{eqnarray*}$$

(C2)

If we set C = 2/α as in equation (4), the derivative with respect to r_s remains the same but

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln \alpha } & = \frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha \left[1 - \alpha \ln \left(\frac{r}{r_{\rm s}} \right)\right] - 1 \right] \, . \end{eqnarray*}$$

(C3)

C2 Orbiting: New model (Model A)

If C = 0, the derivatives with respect to r_s and α are the same as for the Einasto profile with C = 0 (equation C2), and, analogously,

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm t}} & =& \left(\frac{r}{r_{\rm t}} \right)^\beta \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \beta } & =& \frac{1}{\beta } \left(\frac{r}{r_{\rm t}} \right)^\beta \left[ 1 - \beta \ln \left(\frac{r}{r_{\rm t}} \right)\right] \, . \end{eqnarray*}$$

(C4)

If C is set according to equation (8), we get

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm s}} & =& 2 \left(\frac{r}{r_{\rm s}} \right)^\alpha + \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \alpha } & =& \frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha \left[1 - \alpha \ln \left(\frac{r}{r_{\rm s}} \right)\right] - 1 \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm t}} & =& \left(\frac{r}{r_{\rm t}} \right)^\beta - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \beta } & =& \left(\frac{r}{r_{\rm t}} \right)^\beta \left[ \frac{1}{\beta } - \ln \left(\frac{r}{r_{\rm t}} \right)\right] - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[ \frac{1}{\beta } - \ln \left(\frac{r_{\rm s}}{r_{\rm t}} \right)\right] \, . \nonumber \\ \end{eqnarray*}$$

(C5)

C3 Orbiting: New model (Model B)

When setting C = 0, the parameter derivatives for Model B are

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm s}} & =& 2 \left(\frac{r}{r_{\rm s}} \right)^\alpha + \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta \left[ \frac{\beta }{\eta } - 1 \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \alpha } & =& \frac{2}{\alpha } \left(\frac{r}{r_{\rm s}} \right)^\alpha \left[ 1 - \alpha \ln \left(\frac{r}{r_{\rm s}} \right)\right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm t}} & =& \left(\frac{r}{r_{\rm t}} \right)^\beta - \frac{\beta }{\eta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \beta } & =& \left(\frac{r}{r_{\rm t}} \right)^\beta \left[ \frac{1}{\beta } - \ln \left(\frac{r}{r_{\rm t}} \right)\right] + \frac{\beta }{\eta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta \ln \left(\frac{r_{\rm s}}{r_{\rm t}} \right)\nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \eta } & =& \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left(\frac{r}{r_{\rm s}} \right)^\eta \left[ \ln \left(\frac{r}{r_{\rm s}} \right)- \frac{1}{\eta } \right] \, . \end{eqnarray*}$$

(C6)

When setting C as in equation (B3), we get

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm s}} & =& 2 \left(\frac{r}{r_{\rm s}} \right)^\alpha + \left(\frac{\beta }{\eta } - 1 \right) \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[ \left(\frac{r}{r_{\rm s}} \right)^\eta - 1 \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \alpha } & =& \frac{2}{\alpha } \left[ \left(\frac{r}{r_{\rm s}} \right)^\alpha \left[ 1 - \alpha \ln \left(\frac{r}{r_{\rm s}} \right)\right] - 1 \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln r_{\rm t}} & =& \left(\frac{r}{r_{\rm t}} \right)^\beta - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[ \frac{\beta }{\eta } \left[ \left(\frac{r}{r_{\rm s}} \right)^\eta - 1 \right] + 1 \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \beta } & =& \left(\frac{r}{r_{\rm t}} \right)^\beta \left[ \frac{1}{\beta } - \ln \left(\frac{r}{r_{\rm t}} \right)\right] \nonumber \\ && - \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[\frac{1}{\beta } - \ln \left(\frac{r_{\rm s}}{r_{\rm t}} \right)\left[ \frac{\beta }{\eta } \left[ \left(\frac{r}{r_{\rm s}} \right)^\eta - 1 \right] + 1 \right] \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \eta } & =& \frac{1}{\eta } \left(\frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left[ \left(\frac{r}{r_{\rm s}} \right)^\eta \left[ \eta \ln \left(\frac{r}{r_{\rm s}} \right)- 1 \right] + 1 \right]\, . \end{eqnarray*}$$

(C7)

C4 Infalling: New model

We write the derivatives in terms of Q(r) as defined in equation (13),

$$\begin{eqnarray*} \frac{{\rm d} \ln \rho }{{\rm d} \ln \delta _{\rm 1}} & =& f_{\rm m} \left[ 1 - \frac{1}{Q(r)} \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right)^{1 / \zeta } \right] \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln s} & =& -f_{\rm m} \frac{s}{Q(r)} \left(\frac{r}{R} \right)^{s / \zeta } \ln \left(\frac{r}{R_{\rm 200m}} \right)\nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \delta _{\rm max}} & =& f_{\rm m} \frac{1}{Q(r)} \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right)^{1/\zeta } \nonumber \\ \frac{{\rm d} \ln \rho }{{\rm d} \ln \zeta } & =& f_{\rm m} \frac{\ln (Q)}{\zeta } \left[ \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right)^{1 / \zeta } \ln \left(\frac{\delta _{\rm 1}}{\delta _{\rm max}} \right) + s \left(\frac{r}{R} \right)^{s / \zeta } \ln \left(\frac{r}{R} \right) \right] \, . \nonumber \\ \end{eqnarray*}$$

(C8)

All derivatives share the same f_m ≡ 1 − ρ_m/ρ(r) prefactor, which ensures that they vanish at large radii where ρ asymptotically approaches the mean density. This factor is unity when taking the logarithmic derivative of ρ − ρ_m, i.e. of only the excess density over the mean. The parameter derivatives of the model without transition sharpness (Section 2.5) are the same but with ζ = 1.

APPENDIX D: DYNAMICAL PROPERTIES

In Section 5.4, we investigated the space of binding energy, whose structure is determined by the positive gravitational potential of a spherical mass distribution. We find the potential by numerically integrating

$$\begin{eqnarray*} \psi (r) = 4 \pi G \left[ \frac{M(r)}{r} + \int _r^\infty \rho (r^{\prime }) r^{\prime } {\rm d}r^{\prime } \right] \, . \end{eqnarray*}$$

(D1)

For the Einasto profile, we can numerically reproduce the analytical expression for ψ of Retana-Montenegro et al. (2012). Assuming isotropy, the velocity dispersion is found by integrating the Jeans equation,

$$\begin{eqnarray*} \sigma ^2(r) = \frac{1}{\rho (r)} \int _r^\infty \frac{\rho (r^{\prime }) M(r^{\prime })}{r^{\prime 2}} {\rm d}r^{\prime } \, . \end{eqnarray*}$$

(D2)

Analytical expressions for σ exist for some profiles, but they are generally complicated and tend not to extend to Einasto (and other non-power law) forms (e.g. Łokas & Mamon 2001). The distribution function can be found with the expression of Eddington (1916), but this formula involves derivatives of the form dρ/dψ, which are not easy to compute for arbitrary systems. We follow an easier route by converting to derivatives in r (Binney 1982; Baes & Camps 2021),

$$\begin{eqnarray*} f(\varepsilon) = \frac{1}{\sqrt{8} \pi ^2} \int _{r_\varepsilon }^\infty \frac{g(r)}{\sqrt{\varepsilon - \psi (r)}} {\rm d}r \, , \end{eqnarray*}$$

(D3)

where r_ε is the radius where ψ(r_ε) = ε and

$$\begin{eqnarray*} g(r) \equiv \frac{r^2}{M(r)} \left[ \frac{{\rm d}^2 \rho }{{\rm d} r^2} + \frac{{\rm d} \rho }{{\rm d} r} \left(\frac{2}{r} - \frac{4 \pi \rho (r) r^2}{M(r)} \right) \right] \, . \end{eqnarray*}$$

(D4)

The problem of finding f(ε) has now been reduced to computing dρ/dr and d²ρ/dr². We can generally write

$$\begin{eqnarray*} \frac{{\rm d} \rho }{{\rm d} r} = Q(r) \rho \qquad \mathrm{and} \qquad \frac{{\rm d}^2 \rho }{{\rm d} r^2} = \rho \left(\frac{{\rm d} Q}{{\rm d} r} + Q^2 \right) \, . \end{eqnarray*}$$

(D5)

For the Einasto profile, we compute

$$\begin{eqnarray*} Q = -\frac{2}{r} \left(\frac{r}{r_{\rm s}} \right)^\alpha \qquad \mathrm{and} \qquad \frac{{\rm d} Q}{{\rm d} r} = \frac{Q (\alpha - 1)}{r} \, . \end{eqnarray*}$$

(D6)

For Model A, the expressions become somewhat more complicated,

$$\begin{eqnarray*} Q = -\frac{1}{r} \left[ 2 \left(\frac{r}{r_{\rm s}} \right)^\alpha + \left(\frac{r}{r_{\rm t}} \right)^\beta \right] \end{eqnarray*}$$

(D7)

and

$$\begin{eqnarray*} \frac{{\rm d} Q}{{\rm d} r} = -\frac{1}{r} \left(Q + \frac{1}{r} \left[ 2 \alpha \left(\frac{r}{r_{\rm s}} \right)^\alpha + \left(\frac{r}{r_{\rm t}} \right)^\beta \right] \right) \, . \end{eqnarray*}$$

(D8)

Finding the equivalent expressions for Model B or for the DK14 fitting function would be tedious but straightforward; we omit them because their properties would be almost indistinguishable from Model A. For all profiles, we have checked that we can recover the input density profile from the distribution function via the integral over energy,

$$\begin{eqnarray*} \rho (r) = 4 \pi \int _0^{\psi (r)} f(\varepsilon) \sqrt{2[\psi (r) - \varepsilon ]}\ {\rm d}\varepsilon \, . \end{eqnarray*}$$

(D9)

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
December 2022	4
January 2023	21
February 2023	13
March 2023	13
April 2023	15
May 2023	11
June 2023	3
July 2023	2
August 2023	4
September 2023	6
October 2023	10
November 2023	3
December 2023	6
January 2024	30
February 2024	19
March 2024	44
April 2024	21
May 2024	16
June 2024	23
July 2024	26
August 2024	14
September 2024	26
October 2024	38
November 2024	38
December 2024	27
January 2025	39
February 2025	26
March 2025	46
April 2025	32
May 2025	13

Article Contents

A dynamics-based density profile for dark haloes – II. Fitting function Free

ABSTRACT

1 INTRODUCTION

2 PROFILE MODELS

2.1 Generalized exponential profiles

2.2 The Einasto profile

2.3 The DK14 profile

2.4 The new model: truncated exponential

2.5 The infalling profile: previous models

2.6 New infalling model: power law with smooth transition

2.7 Allowed ranges of parameters

3 SIMULATIONS AND METHODS

3.1 N-body simulations

3.2 Dynamically split density profiles

3.3 Fitting procedure

4 RESULTS

4.1 Fits to averaged profiles

4.1.1 The total profile

4.1.2 The orbiting term

4.1.3 The infalling term

4.1.4 Summary

4.2 Fits to individual haloes

5 COMPARISON TO OTHER MODELS

5.1 Models for the orbiting term

5.2 Individual profiles

5.3 Projected profiles and lensing signal

5.4 Dynamical models

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: DETAILS ON FITTING PROCEDURE

A1 Fitting procedure

A2 Parameter Limits

APPENDIX B: MODEL VARIANT WITH CORRECTION AT SCALE RADIUS

APPENDIX C: DERIVATIVES WITH RESPECT TO PARAMETERS

C1 Orbiting: Einasto

C2 Orbiting: New model (Model A)

C3 Orbiting: New model (Model B)

C4 Infalling: New model

APPENDIX D: DYNAMICAL PROPERTIES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

A dynamics-based density profile for dark haloes – II. Fitting function