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

Benedikt Diemer
Department of Astronomy, University of Maryland, College Park, MD 20742, USA
Email<EMAIL_ADDRESS>
(Accepted 2022 December 18. Received 2022 November 15; in original form 2022
May 6)

###### 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 halos 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,
$\rho_{\rm orb}\propto\exp\left[-2/\alpha\ (r/r_{\rm s})^{\alpha}-1/\beta\
(r/r_{\rm t})^{\beta}\right]$, with a five-parameter space of normalization,
physically distinct scale and truncation radii, and $\alpha$ and $\beta$,
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% accuracy across almost all of an expansive
parameter space in halo mass, redshift, cosmology, and accretion rate. When
fixing $\alpha=0.18$ and $\beta=3$, the formula becomes a three-parameter
model that fits individual halos 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.

###### keywords:

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

††pubyear: 2023††pagerange: A dynamics-based density profile for dark haloes –
II. Fitting function–D

##  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 et al., 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 description 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_{\rm 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 et
al. 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 et al., 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 et al.,
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 et al., 2014; More et al., 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
cut-off term with variable truncation radius, sharpness, and asymptotic slope
(a total of $6$ free parameters). This flexible function fits averaged
profiles to 5–10% 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;
Shin et al., 2021; O’Donnell et al., 2022; 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 halos 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%
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_{\rm 200m}$, which encloses an average of $200$ times the mean
density of the Universe, $\rho_{\rm m}$, including all particles (bound or
unbound). We denote the mass enclosed within this radius as $M_{\rm 200m}$. We
mostly express mass as peak height, $\nu\equiv\delta_{\rm c}/\sigma(M_{\rm
200m})$, where $\delta_{\rm c}=1.686$ is the critical density for collapse
(Gunn & Gott, 1972) and $\sigma(M_{\rm 200m})$ is the variance of the linear
power spectrum on the Lagrangian scale corresponding to $M_{\rm 200m}$. We
define the logarithmic mass accretion rate $\Gamma\equiv\Delta\ln M_{\rm
200m}/\Delta\ln a$, measured over one dynamical time (which we define as the
mass-independent halo crossing time, about $5\ {\rm 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)$,

$\rho_{\rm orb}=\rho_{\rm s}e^{S(r)}\,.$ (1)

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

$\gamma\equiv\frac{{\rm d}\ln\rho}{{\rm d}\ln r}=r\frac{{\rm d}S}{{\rm
d}r}\,.$ (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)$,

$S(r)=\int\frac{\gamma(r)}{r}{\rm d}r\,.$ (3)

Clearly, one condition for the fitting function to be useful is that
$\gamma/r$ must integrate to a reasonably simple, analytical expression. For
example, exponential cut-offs in the slope such as $\gamma\propto\exp(r/r_{\rm
t})$ are problematic because their integrals (if solvable) involve Gamma
functions or other complex mathematical expressions.

###  2.2 The Einasto profile

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; $\rho_{\rm 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_{\rm s}=0.2\ R$, $\alpha=0.15$, $r_{\rm t}=R$, $\beta=4$), and the dotted
vertical lines highlight the fiducial $r_{\rm s}$ and $r_{\rm 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 $\gamma(0)=0$ at the halo
centre, but we would need to extend the plot to much smaller radii to see this
effect.

We begin by reviewing the Einasto profile in light of the discussion above.
Its defining feature is a slope of $\gamma=-2(r/r_{\rm s})^{\alpha}$. The
parameter $\alpha$ determines how rapidly the profile steepens, and the
meaning of the scale radius is that $\gamma(r_{\rm s})=-2$. Moreover, the
Einasto profile leads to a core at small radii, $\gamma(0)=0$. The slope
integrates to $S(r)=-(2/\alpha)(r/r_{\rm s})^{\alpha}+C$ We can choose to set
$C=0$, in which case $\rho_{\rm s}$ describes the density at $r=0$, but the
profile is more commonly written with $C=2/\alpha$,

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

such that $S(r_{\rm s})=0$ and thus $\rho(r_{\rm s})=\rho_{\rm s}$. In
practice, we find that the $C\neq 0$ version of the Einasto profile behaves
better in least-squares fits because it does not allow a degeneracy between
$\rho_{\rm s}$, $r_{\rm s}$, and $\alpha$ (although the latter two can still
be degenerate at fixed $\rho_{\rm 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

$\rho_{\rm DK14}=\rho_{\rm Einasto}\left[1+\left(\frac{r}{r_{\rm
t}}\right)^{\beta}\right]^{-\frac{\gamma}{\beta}}\,,$ (5)

where the steepening term suppresses the density at $r\lower
3.01385pt\hbox{$\;\buildrel>\over{\sim}\;$}r_{\rm t}$ as a power law with
slope $-\gamma$, while $\beta$ governs the sharpness of the truncation.
Combined with the Einasto parameters, this function has six free parameters,
but DK14 suggested fixing $(\beta,\gamma)$ 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 halos). The main disadvantage of
equation (5) is that the $\beta$ and $\gamma$ parameters do not have a clear
physical meaning because the asymptotic slope of the overall orbiting profile
depends on $r_{\rm s}$, $r_{\rm t}$, $\alpha$, and $\gamma$ in a complex
fashion. Setting $\beta$ and $\gamma$ to fixed values reduces the profile to a
four-parameter fit with a meaningful truncation radius parameter $r_{\rm t}$,
but it no longer allows for a varying shape of the truncation term.

###  2.4 The new model: truncated exponential

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
$\delta_{\rm max}$ parameter can affect the overdensity at $R$ if it is
sufficiently small.

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,

$\gamma(r)=-2\left(\frac{r}{r_{\rm s}}\right)^{\alpha}-\left(\frac{r}{r_{\rm
t}}\right)^{\beta}\,.$ (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 $\gamma/r$,

$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\,.$ (7)

For simplicity, we could set $C=0$, in which case $\rho_{\rm s}$ would become
the density at the centre. A more elegant form is obtained by setting the
integration constant such that $\rho(r_{\rm s})=\rho_{\rm s}$,

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

and thus

$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]\,.$
(9)

We use this parametrisation throughout because the density at $r_{\rm 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 C.2. 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.

The model does, however, have two minor shortcomings. First, the truncation
term ‘breaks’ the meaning of the scale radius because the slope at $r_{\rm s}$
is now $\gamma(r_{\rm s})=-2-(r_{\rm s}/r_{\rm t})^{\beta}$ instead of $-2$.
The difference is small in most cases since $r_{\rm t}>r_{\rm s}$ and
generally $\beta\gg 1$, but it can manifest itself for extreme parameter
values (red lines in the right two columns of Fig. 1). We note that the
$\gamma(r_{\rm s})=-2$ condition will likely be violated regardless as soon as
an infalling profile is added. Nonetheless, in Appendix B we present a model
variant that enforces $\gamma(r_{\rm s})=-2$ at the expense of an additional
term in $S(r)$. Both models give the same fit and best-fit parameters for
virtually all profiles, but the variant can be preferable in cases where
$r_{\rm s}$ and $\alpha$ 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 $\rho$, namely the enclosed
mass $M(r)$, the projected density $\Sigma(R)$, and the lensing signal
$\Delta\Sigma$, cannot be computed analytically. Similar issues hamper even
the simpler Einasto profile, where $M(r)$ is a relatively complicated
analytical expression (Cardone et al., 2005; Retana-Montenegro et al., 2012)
and $\Sigma(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_{\rm 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,

$\rho(r)=\rho_{\rm m}(z)\left[\delta_{\rm
1}\left(\frac{r}{R}\right)^{-s}+1\right]\,,$ (10)

where $\delta_{\rm 1}$ represents the normalisation at radius $R$ and $-s$ the
slope. We will set $R=R_{\rm 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),

$\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)\,.$
(11)

This expression approaches $\rho\rightarrow\rho_{\rm m}(\delta_{\rm max}+1)$
as $r\rightarrow 0$. The exact value of $\delta_{\rm 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, $\zeta$,

$\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]\,,$ (12)

where we have defined

$Q(r)\equiv\left(\frac{\delta_{\rm 1}}{\delta_{\rm
max}}\right)^{\frac{1}{\zeta}}+\left(\frac{r}{R}\right)^{\frac{s}{\zeta}}\,.$
(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,

$\gamma(r)=-\left(1-\frac{\rho_{\rm
m}}{\rho(r)}\right)\frac{s}{Q(r)}\left(\frac{r}{R}\right)^{\frac{s}{\zeta}}\,,$
(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 $\zeta=0.5$ is a compromise that works for the majority of profiles.
The resulting function,

$\rho(r)=\rho_{\rm m}\left(\frac{\delta_{\rm 1}}{\sqrt{(\delta_{\rm
1}/\delta_{\rm max})^{2}+(r/R)^{2s}}}+1\right)\,,$ (15)

is the infalling profile that we use throughout this paper.

###  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_{\rm
s}<r_{\rm t}$. Some parameters are fixed when fitting individual halo profiles
(Table 1). We discuss the reasons behind the chosen values in detail in
Appendix A.2. In Paper III, we will also present the distributions of best-fit
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: 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 A.2 for details. Symbol | Description | Lower | Upper | Ind.
---|---|---|---|---
Orbiting term |  |
$\rho_{\rm s}/\rho_{\rm m}$ | Overdensity at scale radius | $10$ | $10^{7}$ | free
$r_{\rm s}/R_{\rm 200m}$ | Scale radius | $0.01$ | $0.45$ | free
$r_{\rm t}/R_{\rm 200m}$ | Truncation radius | $0.5$ | $3$ | $<10$
$\alpha$ | Radial evolution of slope | $0.03$ | $0.4$ | $0.18$
$\beta$ | Sharpness of truncation | $0.1$ | $10$ | $3$
Infalling term |  |
$\delta_{\rm 1}$ | Normalization at $R_{\rm 200m}$ | $1$ | $100$ | free
$s$ | Power-law slope | $0.01$ | $4$ | free
$\delta_{\rm max}$ | Central overdensity | $10$ | $2000$ | free
$\zeta$ | Smoothness of transition | $0.5$ | $0.5$ | $0.5$
Model variations |  |
$\rho_{0}/\rho_{\rm m}$ | Overdensity at $r=0$ | $10$ | $10^{20}$ | —
$\eta$ | 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 consists of $14$
simulations of $1024^{3}$ dark matter particles (see Table 1 in Paper I). The
suite covers different box sizes and resolutions, as well as two $\Lambda$CDM
and four self-similar cosmologies. The first $\Lambda$CDM cosmology is that of
the Bolshoi simulation (Klypin et al., 2011), consistent with WMAP7 (Komatsu
et al., 2011), namely a flat $\Lambda$CDM cosmology with $\Omega_{\rm
m}=0.27$, $\Omega_{\rm b}=0.0469$, $\sigma_{8}=0.82$, and $n_{\rm s}=0.95$.
The second is a Planck-like cosmology (Planck Collaboration et al., 2014,
$\Omega_{\rm m}=0.32$, $\Omega_{\rm b}=0.0491$, $h=0.67$, $\sigma_{8}=0.834$,
and $n_{\rm 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 et al., 2008) and allow us to test a wide range of cosmologies with
few simulations because $\Lambda$CDM can be seen as an interpolation between
different power spectrum slopes. For example, redshift trends in $\Lambda$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 $\Lambda$CDM simulations were generated using Camb
(Lewis et al., 2000). They were translated into initial conditions using
2LPTic (Crocce et al., 2006). All simulations were run with Gadget2 (Springel,
2005). We use the phase–space halo finder Rockstar (Behroozi et al., 2013a) to
identify haloes and subhaloes. We construct merger trees by connecting halos
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 halos
with at least $500$ particles within $R_{\rm 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 $\epsilon$ and $l(z)$ are the comoving force
resolution and inter-particle 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 et al., 2019; Mansfield
& Avestruz, 2021). We also impose a limit on the fraction of a halo’s mass
that is unbound, $M_{\rm 200m,all}/M_{\rm 200m,bnd}<1.5$, which excludes halos
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 halos near the truncation. Thus, we exclude
profile bins where $\rho<0.1\ \rho_{\rm 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 halos 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 halos, we again apply a lower radial limit of
$r_{\rm 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 halos with $N_{\rm 200m}\geq 5000$ because their larger
resolved radial range leads to less noisy distributions of the best-fit
parameters. In total, we have fitted the profiles of about 378000 halos in the
WMAP7, Planck, and self-similar cosmologies. The WMAP7 sample shrinks with
redshift from about 45000 at $z=0$ to about 1300 at $z=6$. Similarly, the
self-similar samples shrink from about 120000 halos in the $n=-1$ simulation
to about 13000 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 A.1. 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%
systematic uncertainty added to all bins. For individual halos, we use the sum
of the Poisson error due to (potentially) low number of particles and an
arbitrary 25% systematic error. Each fit is performed multiple times from a
number of initial guesses to ensure that the true minimum is found (see
Appendix A.1 for details).

##  4 Results

Figure 3: Fits to the selected median profiles of halos 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\%$ is linear;
the dotted gray lines mark 5%, 10%, and 20%. 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<\nu<3$). In the third column, this bin is
additionally split by accretion rate. The fourth column compares profiles at
fixed $\nu$ and $\Gamma$ in self-similar cosmologies. The total profiles are
fit to 5% accuracy, with the only exception of the wiggles in low-$\Gamma$
profiles that are caused by particles on their second orbit. The orbiting
profiles are fit to $\sim 10\%$ 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% and 20% accuracy, depending on the
halo sample.  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\%$
accuracy except at the truncation. The fit ignores bins where the density
falls below $\rho=0.1\ \rho_{\rm m}$; the corresponding deviations are driven
by the influence of neighbouring halos 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_{\rm 200m}$ boundary that is used
to define subhalos. However, this issue has no impact on the total profiles.

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\%$ being linear. The dotted lines
highlight differences of $5\%$, $10\%$, and $20\%$ to guide the eye.

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 $\delta_{\rm 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, $\nu$. The
corresponding mass bins range from $M_{\rm 200m}\approx 1.4\times 10^{10}\
{M_{\odot}}$ ($\nu=0.5$) to $M_{\rm 200m}>8\times 10^{14}\ {M_{\odot}}$
($\nu>3$). In the second column we investigate the redshift evolution of the
$2<\nu<3$ bin, which corresponds to $1.4\times 10^{14}<M_{\rm
200m}/{M_{\odot}}<8\times 10^{14}$ at $z=0$ and to $7\times 10^{8}<M_{\rm
200m}/{M_{\odot}}<4.2\times 10^{10}$ at $z=6$. The total profiles are fit
excellently, to within 5% 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, $\Gamma$, than on mass, redshift, or cosmology. Many
of the apparent trends with mass are actually trends in $\Gamma$, which are
convolved with the mass-dependent distribution of accretion rates (larger
halos accrete more actively when structure forms hierarchically). In the third
columns of Figs. 3 and 4, we split the $2<\nu<3$ bin by accretion rates; the
trends are the same for other peak heights (Paper I). An accretion rate of
$\Gamma\lower 3.01385pt\hbox{$\;\buildrel<\over{\sim}\;$}0.5$ corresponds to
pure pseudo-evolution due to the evolving overdensity threshold (Diemer et
al., 2013), meaning that the sample shown in black has undergone essentially
no actual changes to the density profile over the past dynamical time.
Conversely, halos with $\Gamma>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
$\Gamma$-selected profiles almost as well as their $\nu$-selected
counterparts, namely to about 5% within the statistical uncertainties in most
samples and to about 10% for the most extreme (low or high) accretion rates.
One exception are the ‘wiggles’ near the second caustic in low-$\Gamma$
profiles, where the fit can deviate by up to 20%. 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 halos by both $\nu$ and $\Gamma$ 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
$\Lambda$CDM, the total profiles are fit to 5% or better (subject to the
trends with $\Gamma$ 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 halos 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 halos, the
resulting differences are much more apparent in the mean than in the median
profiles. These issues caution us not to over-interpret fitting errors at very
low densities, which can reach $\rho=10^{-5}\rho_{\rm m}$ in some mean
profiles. We thus ignore densities below $\rho<0.1\ \rho_{\rm m}$ in fits to
the mean orbiting profiles.

The $\nu$-selected orbiting profiles are fit to 10% 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
$\Gamma$-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\lower
3.01385pt\hbox{$\;\buildrel<\over{\sim}\;$}0.5\ R_{\rm 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_{\rm t}$. However, the detailed profile shapes are diverse and depend on
$\nu$, $\Gamma$, and $n$, as evidenced by rows 7–9 of Figs. 3 and 4.

Our model fits most median infalling profiles to about 10% and some
$\Gamma$-selected samples to 20% 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% for low-mass samples. One reason is that the mean profiles are noticeably
affected by nearby neighbours, and thus by the definition of subhalos as
residing within $R_{\rm 200m}$ (Paper I). This arbitrary boundary leads to
sharp features around $R_{\rm 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% 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 halos), 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 halos

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 19). 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-fit parameters. The top row shows halos from the
WMAP7 cosmology at $z=0$ with ascending $\Gamma$, which leads to a decreasing
$r_{\rm t}/R_{\rm 200m}$. The normalization, slope, and asymptotic values of
the infalling profiles differ dramatically between individual halos. The
bottom row shows a few extreme cases. First, a halo at $z=4$ has a very high
accretion rate, $\Gamma\approx 8$, which leads to a truncation radius at the
lower limit of $r_{\rm t}\geq 0.55\ R_{\rm 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 ($\alpha$ instead of $r_{\rm t}$) can fit
low-$\Gamma$ profiles, it cannot capture the sharp truncation (dotted lines).

One requirement for our fitting function was that it should fit the total
profiles of individual halos, as well as their dynamically split components.
For each halo in our sample of 378000, we separately fit equation (9) to the
orbiting term (fixing $\alpha=0.18$ and $\beta=3$) and equation (15) to the
infalling term (Section 3.3). We then use the best-fit 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_{\rm s}$ and $\alpha$, or to reliably determine
$r_{\rm 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 $\delta_{\rm 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 $\rho_{\rm s}$ and $r_{\rm s}$ in the combined fit, meaning that
$r_{\rm t}$, $\delta_{\rm 1}$, $\delta_{\rm 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-fit 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 halos.
Our model (dashed lines) fits the profiles well within statistical error and
stochastic variations, except for the bottom-right panel which was selected to
show a poor fit. The top row shows halos with increasing accretion rates. As
discussed in Paper I, their peak height (mass) has very little influence on
the profile shapes, but the increasing $\Gamma$ clearly manifests itself in
decreasing values of $r_{\rm t}/R_{\rm 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_{\rm t}$ reaches the minimum imposed in our fits,
$0.55\ R_{\rm 200m}$. However, we have not observed examples where the true
$r_{\rm 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
$\delta_{\rm 1}\approx 45$ at $R_{\rm 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\approx 0.2\ R_{\rm
200m}$. Conversely, the third panel shows a halo with flat infalling and
orbiting profiles, which lead to $r_{\rm s}$ reaching the maximum of $0.45\
R_{\rm 200m}$. Finally, the bottom right 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.

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

Figure 6: Comparison of different fitting functions for the orbiting term. The
profiles in rows 1 and 4 are rescaled by $r^{2}$ 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 halos with
$2<\nu<3$ (left) and the mean profiles of a sub-sample additionally selected
to have $2<\Gamma<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 $\beta$ and $\gamma$
as free parameters (yellow dot-dashed), DK14 with the slopes fixed to their
suggested values for $\nu$\- and $\Gamma$-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 over-
fit by the free DK14 model (right column).  Figure 7: Comparison of the
relative fit quality of equation (9) and the NFW and Einasto models, for all
378000 halos 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 $\chi^{2}/N_{\rm dof}$ between each model and our new
model (with three free parameters, $\rho_{\rm s}$, $r_{\rm s}$, and $r_{\rm
t}$). The points with error bars mark the mean ratio and 68% 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 $\alpha$ is not sufficient to fit the truncation. When fitting
the total profiles but keeping $\alpha$ and $r_{\rm t}$ fixed based on the
separate fit (right panel), the fit quality is improved by 10% on average (and
a median of $5\%$). This result holds even when fitting only to the total term
with no knowledge of the split profiles (yellow). Fixing $\alpha$ in the
Einasto fits (light blue) degrades the orbiting fit but also reduces the
number of degrees of freedom so that the $\chi^{2}/N_{\rm dof}$ difference in
the total fits remains similar. The two-parameter NFW model compares less well
than the Einasto profile (dark blue).

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 halos 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).

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
cut-off 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\lower 3.01385pt\hbox{$\;\buildrel>\over{\sim}\;$}r_{\rm t}$. This issue
persists for profiles with other slope transitions such as Burkert (1995) or
superNFW (Lilley et al., 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% accuracy, but it fails to capture the
orbiting term because it approaches a fixed slope of $-3-2n$. Here $n\approx
11$, a value much larger than the $n\approx 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_{\rm s}$ and $r_{\rm 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_{\rm 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 et al., 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_{\rm 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 et al., 2016) or coreEinasto (Lazar et al., 2020) models, but
these modifications obviously do not change the fit at large radii. Springel &
White (1999) suggested $\rho\propto\exp[-(r-R_{\rm 200c})/r_{\rm s}]$ to model
tidal truncation, but this term does not introduce a new scale and thus cannot
capture variations in $r_{\rm t}/r_{\rm 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,
$\beta$, and the asymptotic slope of the steepening term, $\gamma$, are free
parameters (hereafter ‘DK14-6’, yellow dot-dashed lines). In the second, they
are fixed to the optimal values recommended by DK14, namely $\beta=4$ and
$\gamma=8$ for the $\nu$-selected sample in the left column and $\beta=6$ and
$\gamma=4$ for $\Gamma$-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 $\gamma$, 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
$\chi^{2}/N_{\rm dof}$ values. The conclusions depend somewhat on whether
samples are selected by $\nu$ only or also $\Gamma$, 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
describe 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 $\Delta_{\chi}\equiv\Delta\log_{10}(\chi^{2}/N_{\rm dof})$ of up
to $2$. Compared to the DK14-6 form that has one more free parameter, the
median $\nu$-selected profiles tend to be better fit by the new form, and the
mean profiles by DK14-6. Most $\Gamma$-selected samples are better fit by
DK14-6, but at the expense of an extreme range of best-fit $\beta$ and
$\gamma$. The new model better fits particularly sharp truncations in median
samples with $\Gamma\lower 3.01385pt\hbox{$\;\buildrel>\over{\sim}\;$}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 halos, 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 $\chi^{2}/N_{\rm dof}$. However, giving the DK14
model its full freedom of $6$ parameters often leads to meaningless values of
$\beta$ and $\gamma$. 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 halos 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 A.1), our best option is to compare the relative $\chi^{2}/N_{\rm
dof}$ of different fits. To this end, we repeat the procedure described in
Section 4.2 with three models: NFW, Einasto with free $\alpha$, and Einasto
with $\alpha=0.18$. In Fig. 7 we compare these fits for all 378000 halos in
our individual sample. While there are mild trends when we break up the sample
by cosmology and redshift, the overall picture remains unchanged.

The most ‘fair’ comparison is between the new model and an Einasto profile
with free $\alpha$ because both have three parameters. When fitting only the
orbiting profile, the new model has a lower $\chi^{2}/N_{\rm dof}$ by a median
factor of four (green histogram in the left 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-$\Gamma$ halos 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 panel of Fig. 7. We test two variations of the fits. First
we vary only $\rho_{\rm s}$ and $r_{\rm s}$ in the combined fit, keeping
$\alpha$, $r_{\rm t}$, and the parameters of the infalling profile to their
values from the separate fits (green). The new model now improves
$\chi^{2}/N_{\rm dof}$ by about 10% 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 $\alpha=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 $\chi^{2}/N_{\rm 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% and 34%. We have
confirmed that these results remain unchanged when directly optimizing
$\chi^{2}$ instead of the Cauchy loss function $\ln(1+\chi^{2})$ (Appendix
A.1).

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 halos and their best-fit parameters.

###  5.3 Projected profiles and lensing signal

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}M_{\odot}$ and $c_{\rm 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_{\rm t}$ are easily distinguishable by eye. The
differences are smaller but propagate to a much larger radial range in
$\Delta\Sigma$ because $\Sigma(r)$ is subtracted from the averaged surface
density.

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,
$\Sigma(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}(<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}M_{\odot}$ and $c_{\rm 200m}=5$ at
$z=0$. We set $\alpha=0.18$ and $\beta=0.3$, but we vary $r_{\rm t}$ between
values corresponding to low and high accretion rates (about $0.7<r_{\rm
t}/R_{\rm 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).

###  5.4 Dynamical models

Figure 9: Dynamical properties of our new model (dark blue), a King (1966)
profile (light blue, with a comparable cut-off radius), as well as the NFW and
Einasto forms (gray). The panels show the density $\rho$, enclosed mass $M$,
potential $\psi$, isotropic velocity dispersion $\sigma^{2}$, and distribution
function $f(\varepsilon)$ (from top to bottom). The latter falls to zero for
unbound particles with $\varepsilon<0$; conversely, $\varepsilon=\psi(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 cut-off 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.

Assuming that halos 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, $\sigma^{2}$, and the phase-space distribution function,
$f(\varepsilon)$. 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 et al., 2010; Okumura
et al., 2018; Adhikari et al., 2019; Hamabata et al., 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 $\psi(r)$ to be the (positive)
gravitational potential of the halo and $\varepsilon\equiv\psi(r)-v^{2}/2$ to
be the relative binding energy of particles. Particles with $\varepsilon<0$
would be unbound and thus lead to $f(\varepsilon<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 $\varepsilon=\psi(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
$\rho_{0}$, and the gravitational constant $G$. The results thus depend only
on the relative scales $r_{\rm s}/R$ and $r_{\rm t}/R$, as well as the
steepening parameters $\alpha$ and $\beta$. We choose a representative set of
parameters for all models, namely $r_{\rm s}/R=0.2$, $\alpha=0.18$, $r_{\rm
t}/R=1$, and $\beta=4$. We observe that $M(r)$ flattens as a consequence of
the truncation in $\rho$, and that the potential approaches $\psi\propto 1/r$
at a smaller radius than for the Einasto profile.

The density cut-off 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 et al., 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 profile 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 $\beta$ 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(\varepsilon)$ 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

$f_{\rm king}(\varepsilon)=f_{0}\left(e^{(\varepsilon-\varepsilon_{\rm
t})/\sigma^{2}}-1\right)\,,$ (16)

where $f_{0}$ is a normalization, $\sigma$ the velocity dispersion of the
isothermal sphere, and $\varepsilon_{\rm t}$ the energy where the object is
truncated ($f\rightarrow 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_{\rm
t}\equiv\varepsilon_{\rm t}/[\sigma^{2}\psi(0)]$. The light-blue lines in Fig.
9 show a King profile with $Z_{\rm t}=0.2$, a cut-off energy chosen to roughly
match the real-space truncation radius of $r_{\rm t}=0.2\ R$. The position and
sharpness of the truncation can be adjusted using $Z_{\rm 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\ll 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 $\gamma=-4$
at large radii, regardless of the $\phi_{0}$ parameter that controls the
distribution of particle energies (Williams & Hjorth, 2010; Williams et al.,
2010). Similarly, Pontzen & Governato (2013) derive a distribution function by
maximising 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, $\rho_{\rm
orb}\propto\exp\left[-2/\alpha\ (r/r_{\rm s})^{\alpha}-1/\beta\ (r/r_{\rm
t})^{\beta}\right]$. 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 $\Delta\Sigma$ integrals. Our
main conclusions are as follows.

  1. 1.

The new model fits mean and median total profiles to roughly $\pm 5\%$ across
a vast range of radius, halo mass, redshift, and cosmology. When halos are
also selected by accretion rate, the fit quality can degrade to $\pm 20\%$
near the transition between the orbiting and infalling profiles.

  2. 2.

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.

  3. 3.

By fixing $\alpha=0.18$ and $\beta=3$, the orbiting model becomes a
3-parameter fit that captures the profiles of individual halos more accurately
than Einasto profiles (on average).

  4. 4.

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.

  5. 5.

Different truncation radii lead to clear differences in the projected surface
density profiles and in the lensing signal $\Delta\Sigma$.

  6. 6.

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-fit 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 subhalos
(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 et al. (2014) Adhikari S., Dalal N., Chamberlain R. T., 2014, JCAP, 11, 19
  * Adhikari et al. (2019) Adhikari S., Dalal N., More S., Wetzel A., 2019, ApJ, 878, 9
  * Adhikari et al. (2021) Adhikari S., et al., 2021, ApJ, 923, 37
  * Amorisco (2021) Amorisco N. C., 2021, arXiv e-prints, p. arXiv:2111.01148
  * An & Zhao (2013) An J., Zhao H., 2013, MNRAS, 428, 2805
  * Aung et al. (2021) Aung H., Nagai D., Rozo E., García R., 2021, MNRAS, 502, 1041
  * Aung et al. (2022) Aung H., Nagai D., Rozo E., Wolfe B., Adhikari S., 2022, arXiv e-prints, p. arXiv:2204.13131
  * Baes & Camps (2021) Baes M., Camps P., 2021, MNRAS, 503, 2955
  * Baes & Ciotti (2019) Baes M., Ciotti L., 2019, A&A, 626, A110
  * Baltz et al. (2009) Baltz E. A., Marshall P., Oguri M., 2009, J. Cosmology Astropart. Phys., 1, 15
  * Baxter et al. (2017) Baxter E., et al., 2017, ApJ, 841, 18
  * Becker & Kravtsov (2011) Becker M. R., Kravtsov A. V., 2011, ApJ, 740, 25
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013b, ApJ, 763, 18
  * Beraldo e Silva et al. (2014) Beraldo e Silva L., Lima M., Sodré L., Perez J., 2014, Phys. Rev. D, 90, 123004
  * Bertschinger (1985) Bertschinger E., 1985, ApJS, 58, 39
  * Betancort-Rijo et al. (2006) Betancort-Rijo J. E., Sanchez-Conde M. A., Prada F., Patiri S. G., 2006, ApJ, 649, 579
  * Bianconi et al. (2021) Bianconi M., Buscicchio R., Smith G. P., McGee S. L., Haines C. P., Finoguenov A., Babul A., 2021, ApJ, 911, 136
  * Binney (1982) Binney J., 1982, MNRAS, 200, 951
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Bond et al. (1996) Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603
  * Bose & Loeb (2021) Bose S., Loeb A., 2021, ApJ, 912, 114
  * Brown et al. (2020) Brown S. T., McCarthy I. G., Diemer B., Font A. S., Stafford S. G., Pfiefer S., 2020, MNRAS, 495, 4994
  * Burkert (1995) Burkert A., 1995, ApJ, 447, L25
  * Callingham et al. (2020) Callingham T. M., Cautun M., Deason A. J., Frenk C. S., Grand R. J. J., Marinacci F., Pakmor R., 2020, MNRAS, 495, 12
  * Cardone et al. (2005) Cardone V. F., Piedipalumbo E., Tortora C., 2005, MNRAS, 358, 1325
  * Chang et al. (2018) Chang C., et al., 2018, ApJ, 864, 83
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * Courteau et al. (2014) Courteau S., et al., 2014, Reviews of Modern Physics, 86, 47
  * Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369
  * Dacunha et al. (2022) Dacunha T., Belyakov M., Adhikari S., Shin T.-h., Goldstein S., Jain B., 2022, MNRAS, 512, 4378
  * Dehnen (1993) Dehnen W., 1993, MNRAS, 265, 250
  * Dekel et al. (2017) Dekel A., Ishai G., Dutton A. A., Maccio A. V., 2017, MNRAS, 468, 1005
  * Dhar (2021) Dhar B. K., 2021, MNRAS, 504, 4583
  * Dhar & Williams (2010) Dhar B. K., Williams L. L. R., 2010, MNRAS, 405, 340
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415
  * Diemand & Kuhlen (2008) Diemand J., Kuhlen M., 2008, ApJ, 680, L25
  * Diemer (2017) Diemer B., 2017, ApJS, 231, 5
  * Diemer (2018) Diemer B., 2018, The Astrophysical Journal Supplement Series, 239, 35
  * Diemer (2020) Diemer B., 2020, ApJS, 251, 17
  * Diemer (2022) Diemer B., 2022, MNRAS, 513, 573
  * Diemer & Kravtsov (2014) Diemer B., Kravtsov A. V., 2014, ApJ, 789, 1
  * Diemer & Kravtsov (2015) Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108
  * Diemer et al. (2013) Diemer B., More S., Kravtsov A. V., 2013, ApJ, 766, 25
  * Dooley et al. (2014) Dooley G. A., Griffen B. F., Zukin P., Ji A. P., Vogelsberger M., Hernquist L. E., Frebel A., 2014, ApJ, 786, 50
  * Doroshkevich & Shandarin (1978) Doroshkevich A. G., Shandarin S. F., 1978, Soviet Ast., 22, 653
  * Drakos et al. (2017) Drakos N. E., Taylor J. E., Benson A. J., 2017, MNRAS, 468, 2345
  * Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
  * Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
  * Eckert et al. (2022) Eckert D., Ettori S., Robertson A., Massey R., Pointecouteau E., Harvey D., McCarthy I. G., 2022, A&A, 666, A41
  * Eddington (1916) Eddington A. S., 1916, MNRAS, 76, 572
  * Efstathiou et al. (1988) Efstathiou G., Frenk C. S., White S. D. M., Davis M., 1988, MNRAS, 235, 715
  * Einasto (1965) Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Einasto (1969) Einasto J., 1969, Astrophysics, 5, 67
  * Erwin (2015) Erwin P., 2015, ApJ, 799, 226
  * Fielder et al. (2020) Fielder C. E., Mao Y.-Y., Zentner A. R., Newman J. A., Wu H.-Y., Wechsler R. H., 2020, MNRAS, 499, 2426
  * Fong & Han (2021) Fong M., Han J., 2021, MNRAS, 503, 4250
  * Freundlich et al. (2020) Freundlich J., et al., 2020, MNRAS, 499, 2912
  * Fukushige & Makino (2001) Fukushige T., Makino J., 2001, ApJ, 557, 533
  * Gao et al. (2008) 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
  * García et al. (2021) García R., Rozo E., Becker M. R., More S., 2021, MNRAS, 505, 1195
  * García et al. (2022) García R., Salazar E., Rozo E., Adhikari S., Aung H., Diemer B., Nagai D., Wolfe B., 2022, arXiv e-prints, p. arXiv:2207.11827
  * Genina et al. (2018) Genina A., et al., 2018, MNRAS, 474, 1398
  * Ghigna et al. (2000) Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 2000, ApJ, 544, 616
  * Graham et al. (2006) Graham A. W., Merritt D., Moore B., Diemand J., Terzić B., 2006, AJ, 132, 2701
  * Gunn & Gott (1972) Gunn J. E., Gott III J. R., 1972, ApJ, 176, 1
  * Hamabata et al. (2019) Hamabata A., Oguri M., Nishimichi T., 2019, MNRAS, 489, 1344
  * Hansen & Moore (2006) Hansen S. H., Moore B., 2006, New Astron., 11, 333
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hjorth & Williams (2010) Hjorth J., Williams L. L. R., 2010, ApJ, 722, 851
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Jaffe (1983) Jaffe W., 1983, MNRAS, 202, 995
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Knollmann et al. (2008) Knollmann S. R., Power C., Knebe A., 2008, MNRAS, 385, 545
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Lazar et al. (2020) Lazar A., et al., 2020, MNRAS, 497, 2393
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Lilley et al. (2018) Lilley E. J., Evans N. W., Sanders J. L., 2018, MNRAS, 476, 2086
  * Łokas & Mamon (2001) Łokas E. L., Mamon G. A., 2001, MNRAS, 321, 155
  * Lucie-Smith et al. (2022a) Lucie-Smith L., Peiris H. V., Pontzen A., Nord B., Thiyagalingam J., Piras D., 2022a, Phys. Rev. D, 105, 103533
  * Lucie-Smith et al. (2022b) Lucie-Smith L., Adhikari S., Wechsler R. H., 2022b, MNRAS, 515, 2164
  * Ludlow & Angulo (2017) Ludlow A. D., Angulo R. E., 2017, MNRAS, 465, L84
  * Ludlow et al. (2011) Ludlow A. D., Navarro J. F., White S. D. M., Boylan-Kolchin M., Springel V., Jenkins A., Frenk C. S., 2011, MNRAS, 415, 3895
  * Ludlow et al. (2013) Ludlow A. D., et al., 2013, MNRAS, 432, 1103
  * Ludlow et al. (2016) 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
  * Ludlow et al. (2019) Ludlow A. D., Schaye J., Bower R., 2019, MNRAS, 488, 3663
  * Mamon & Łokas (2005) Mamon G. A., Łokas E. L., 2005, MNRAS, 362, 95
  * Mamon et al. (2010) Mamon G. A., Biviano A., Murante G., 2010, A&A, 520, A30
  * Mansfield & Avestruz (2021) Mansfield P., Avestruz C., 2021, MNRAS, 500, 3309
  * Merritt et al. (2006) Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ, 132, 2685
  * Michie (1963) Michie R. W., 1963, MNRAS, 125, 127
  * Moore et al. (1999) Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999, MNRAS, 310, 1147
  * More et al. (2015) More S., Diemer B., Kravtsov A. V., 2015, ApJ, 810, 36
  * More et al. (2016) More S., et al., 2016, ApJ, 825, 39
  * Murata et al. (2020) Murata R., Sunayama T., Oguri M., More S., Nishizawa A. J., Nishimichi T., Osato K., 2020, PASJ, 72, 64
  * Navarro et al. (1995) Navarro J. F., Frenk C. S., White S. D. M., 1995, MNRAS, 275, 720
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Navarro et al. (2004) Navarro J. F., et al., 2004, MNRAS, 349, 1039
  * Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
  * O’Donnell et al. (2022) O’Donnell C., Behroozi P., More S., 2022, MNRAS, 509, 3285
  * O’Neil et al. (2021) O’Neil S., Barnes D. J., Vogelsberger M., Diemer B., 2021, MNRAS, 504, 4649
  * O’Neil et al. (2022) O’Neil S., Borrow J., Vogelsberger M., Diemer B., 2022, MNRAS, 513, 835
  * Oguri & Hamana (2011) Oguri M., Hamana T., 2011, MNRAS, 414, 1851
  * Okumura et al. (2018) Okumura T., Nishimichi T., Umetsu K., Osato K., 2018, Phys. Rev. D, 98, 023523
  * Peng et al. (2010) Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2010, AJ, 139, 2097
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Pontzen & Governato (2013) Pontzen A., Governato F., 2013, MNRAS, 430, 121
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Prada et al. (2006) Prada F., Klypin A. A., Simonneau E., Betancort-Rijo J., Patiri S., Gottlöber S., Sanchez-Conde M. A., 2006, ApJ, 645, 1001
  * Read et al. (2016) Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS, 459, 2573
  * Retana-Montenegro et al. (2012) Retana-Montenegro E., van Hese E., Gentile G., Baes M., Frutos-Alfaro F., 2012, A&A, 540, A70
  * Ricotti et al. (2007) Ricotti M., Pontzen A., Viel M., 2007, ApJ, 663, L53
  * Schneider et al. (2019) 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
  * Shapiro et al. (1999) Shapiro P. R., Iliev I. T., Raga A. C., 1999, MNRAS, 307, 203
  * Shin & Diemer (2022) Shin T.-h., Diemer B., 2022, arXiv e-prints, p. arXiv:2210.14262
  * Shin et al. (2019) Shin T., et al., 2019, MNRAS, 487, 2900
  * Shin et al. (2021) Shin T., et al., 2021, MNRAS, 507, 5758
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & White (1999) Springel V., White S. D. M., 1999, MNRAS, 307, 162
  * Stadel et al. (2009) Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21
  * Sugiura et al. (2020) Sugiura H., Nishimichi T., Rasera Y., Taruya A., 2020, MNRAS, 493, 2765
  * Tasitsiomi et al. (2004) Tasitsiomi A., Kravtsov A. V., Gottlöber S., Klypin A. A., 2004, ApJ, 607, 125
  * Tavio et al. (2008) Tavio H., Cuesta A. J., Prada F., Klypin A. A., Sanchez-Conde M. A., 2008, arXiv:0807.3027,
  * Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
  * Tomooka et al. (2020) Tomooka P., Rozo E., Wagoner E. L., Aung H., Nagai D., Safonova S., 2020, MNRAS, 499, 1291
  * Udrescu et al. (2019) Udrescu S. M., Dutton A. A., Macciò A. V., Buck T., 2019, MNRAS, 482, 5259
  * Umetsu (2020) Umetsu K., 2020, A&ARv, 28, 7
  * Umetsu & Diemer (2017) Umetsu K., Diemer B., 2017, ApJ, 836, 231
  * Velliscig et al. (2014) 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
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wagner (2020) Wagner J., 2020, General Relativity and Gravitation, 52, 61
  * Wang et al. (2020) Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2020, Nature, 585, 39
  * Wechsler & Tinker (2018) Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Williams & Hjorth (2010) Williams L. L. R., Hjorth J., 2010, ApJ, 722, 856
  * Williams et al. (2010) Williams L. L. R., Hjorth J., Wojtak R., 2010, ApJ, 725, 282
  * Xhakaj et al. (2020) Xhakaj E., Diemer B., Leauthaud A., Wasserman A., Huang S., Luo Y., Adhikari S., Singh S., 2020, MNRAS, 499, 3534
  * Xhakaj et al. (2022) Xhakaj E., Leauthaud A., Lange J., Hearin A., Diemer B., Dalal N., 2022, MNRAS, 514, 2876
  * Zhao (1996) Zhao H., 1996, MNRAS, 278, 488
  * Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Börner G., 2009, ApJ, 707, 354
  * de Blok (2010) de Blok W. J. G., 2010, Advances in Astronomy, 2010, 789293

##  Appendix A Details on fitting procedure

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

###  A.1 Fitting procedure

A first important choice is the definition of the residual. The conventional
approach of minimizing the square of $(\rho_{\rm fit}-\rho_{\rm data})/\sigma$
works poorly in the outer profiles because the summed $\chi^{2}$ is dominated
by the much higher densities at small radii. Some authors have chosen to
minimize $r^{2}(\rho_{\rm fit}-\rho_{\rm data})/\sigma$ 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

$\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})}\,,$ (17)

where $\rho_{\rm i,fit}$ is the fit density in bin $i$, $\rho_{\rm i,data}$ is
the density of the simulated profile (individual or averaged), and
$\sigma_{\rm i}$ is the uncertainty in $\rho_{\rm 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 $\sigma_{\rm i}$ of each
profile bin. The bootstrap estimate, $\sigma_{\rm 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% in
quadrature,

$\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}}\,.$
(18)

The systematic error does have some effect on on the best-fit 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% accuracy, this value is a reasonable approximation to the systematic
uncertainty.

We minimize $\chi^{2}$ 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 $\chi^{2}$ 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^{-3}$. 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 $\chi^{2}$,
we set the new results as the initial guess. If the total gradient at the
solution (‘optimality’) is larger than $10^{-3}$, 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
$\chi^{2}$ of all attempts and re-run the fit with a higher accuracy of
$10^{-5}$ 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 $\chi^{2}/N_{\rm 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, $\chi^{2}/N_{\rm 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,

$\sigma_{\rm i,ind}=\rho_{\rm i,data}\sqrt{\frac{1}{N_{\rm i,ptl}}+\sigma_{\rm
sys,ind}^{2}}\,,$ (19)

where $N_{\rm i,ptl}$ is the number of particles in the given bin and
$\sigma_{\rm sys,ind}$ is a fixed systematic error. Since most bins contain
$N_{\rm i,ptl}\gg 1$ particles, the normalization of $\chi^{2}$ is largely
determined by the somewhat arbitrary value for $\sigma_{\rm 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 be dominated by a few bins with high $N_{\rm ptl}$. We set
$\sigma_{\rm sys,ind}=0.25$, meaning that the systematic error dominates for
bins with $N_{\rm i,ptl}>16$. With this value, the median $\chi^{2}/N_{\rm
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-fit parameters.

###  A.2 Parameter Limits

Figure 10: Same as Fig. 1, but showing the modified Model B with $\eta=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 $\gamma(r_{\rm s})=-2$
regardless of the values of $r_{\rm t}$ and $\beta$, Model A can exhibit
slight deviations from this value if $r_{\rm t}$ and/or $\beta$ are small.

The normalization, $\rho_{\rm s}$, is tightly related to the formation
redshift of halos (Navarro et al., 1997; Ludlow et al., 2013) and thus cannot
vary arbitrarily. Nonetheless, we allow a generous range because $\rho_{\rm
s}$ is well-constrained in all fits. If we instead use the density at the
centre, $\rho_{0}$, 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.,
$\gamma\approx-2$ over a wide range of radii. In such cases, the fit tends to
shift $r_{\rm s}$ towards large values that can even exceed $r_{\rm t}$, which
is technically a valid solution because of the symmetry between the $r_{\rm
s}$ and $r_{\rm 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_{\rm t}>r_{\rm s}$. In particular, we set $0.01<r_{\rm
s}/R_{\rm 200m}<0.45$ and $0.5<r_{\rm t}/R_{\rm 200m}<3$. The range of $r_{\rm
s}$ limits the concentration to $2.2<c_{\rm 200m}<100$, but halos outside of
this range would correspond to extreme outliers in the population (e.g.,
Diemer & Kravtsov, 2015); we never observe a preference for $r_{\rm s}\lower
3.01385pt\hbox{$\;\buildrel>\over{\sim}\;$}0.4\ R_{\rm 200m}$ in well-
constrained fits. The lower limit corresponds roughly to the lowest radii
where our averaged profiles are resolved, meaning that $r_{\rm s}$ cannot be
measured reliably near $0.01\ R_{\rm 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_{\rm 200m}$ in any
profiles, and we generally find little evidence for $r_{\rm t}>2\ R_{\rm
200m}$. If $r_{\rm t}$ pushes against the upper limit, the profile does not
experience a clear truncation, meaning that the values of $r_{\rm t}$ and
$\beta$ 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<\alpha<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 halos with high accretion rates (Paper I). The Einasto profile
can suffer from a well-known degeneracy between $r_{\rm s}$ and $\alpha$ if
profiles do not extend to small radii or if their slope evolves slowly with
radius (e.g., Ricotti et al., 2007; Udrescu et al., 2019). While virtually all
averaged profiles do constrain $\alpha$ to some extent, most individual halo
profiles do not. In particular, the average fit quality remains roughly
constant within the range $0.1<\alpha<0.2$, highlighting that individual
profiles simply do not contain sufficient information to determine $\alpha$.
We could set it based on an $\nu$-$\alpha$ 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 $\alpha\approx 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 halos with very high $\Gamma$ do indeed show a mild
preference for low $\alpha$, the differences in fit quality are modest at
best. Thus, we simply fix $\alpha=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., $\alpha=0.16$, Wang et al., 2020).

We have no intuitive prior for the transition sharpness $\beta$, but the
profiles shown in Paper I range from gradual to very sharp transitions. Thus,
we allow $\beta$ to fluctuate between $0.1$ and $10$. While some fits
technically prefer even larger values, those correspond to an essentially
instantaneous cut-off, where $\beta$ is not well defined. We might worry that
the $r_{\rm t}$-$\beta$ term suffers from a degeneracy similar to that of the
$r_{\rm s}$-$\alpha$ term, but we find that $r_{\rm t}$ and $\beta$ are not
intrinsically degenerate and both well-constrained for most averaged profiles
(Paper III). In individual halos, however, the orbiting profiles commonly
become unresolved at the radii that most constrain $\beta$. We thus fix
$\beta=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_{\rm t}/R_{\rm 200m}\leq 10$.

In the infalling profile, we limit the normalization to $1<\delta_{\rm
1}<100$. Added to the mean density of the Universe, the lower limit
corresponds to $\rho_{\rm inf}(R_{\rm 200m})=2\rho_{\rm 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 $\delta_{\rm max}$ reaches values between $10$ and $1000$ in the
profiles where it is constrained, but we allow $1<\delta_{\rm max}<10^{4}$ 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 $\delta_{\rm max}$ because the infalling profile is subdominant to the
orbiting term at $r<R_{\rm 200m}$, meaning that $\delta_{\rm max}$ is
unconstrained. In a situation where only the total profile is known,
$\delta_{\rm 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 $\gamma(r_{\rm 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 $\gamma(r_{\rm
s})=-2$, but we need to be careful. For example, if we add $(r_{\rm s}/r_{\rm
t})^{\beta}$ into the slope of Model A, $\gamma(r)$ is still integrable and
adds the term $(r_{\rm s}/r_{\rm t})^{\beta}\ln(r)$ to $S(r)$ from equation
(7). While this new model would now satisfy $\gamma(r_{\rm s})=-2$, $\gamma$
has been renormalized at all radii with a fixed term that survives all the way
to the centre of the halo, $\gamma(0)=(r_{\rm s}/r_{\rm t})^{\beta}$. 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_{\rm s}/r_{\rm
t})^{\beta}$ at $r=r_{\rm s}$,

$\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}\,.$ (20)

We integrate $\gamma/r$ to find

$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\,.$ (21)

Once again, we can set $C=0$ or ensure $S(r_{\rm s})=0$ by setting

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

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

$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]+\frac{1}{\eta}\left(\frac{r_{\rm s}}{r_{\rm
t}}\right)^{\beta}\left[\left(\frac{r}{r_{\rm s}}\right)^{\eta}-1\right]\,.$
(23)

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, $\eta$ is a nuisance
parameter that determines how fast this decay happens: the larger $\eta$, the
faster the slope approaches zero at small radii. On the other hand, $\eta$
also introduces an undesirable change in the slope at larger radii (where
$r/r_{\rm s}\lower 3.01385pt\hbox{$\;\buildrel>\over{\sim}\;$}r_{\rm s}/r_{\rm
t}$). However, as long as $\beta\gg 1$ and $\eta<1$, this correction is barely
noticeable. We fix $\eta=0.1$, but the Model B fits are insensitive in the
range $0.01<\eta<0.2$. In some fits to mean orbiting-only profiles, the best-
fit value of $\alpha$ depends slightly on $\eta$, but this dependence
disappears when fitting the entire profile or median orbiting profiles.

Fig. 10 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_{\rm t}$
and/or small $\beta$. The example of $r_{\rm t}/R_{\rm 200m}=0.3$ is outside
of the allowed parameter range for our fits ($r_{\rm t}\geq 0.5R_{\rm 200m}$),
but we do allow very small values of $\beta$. Model A fixes the slope at
$r_{\rm t}$ to be $\gamma(r_{\rm t})=-2(r_{\rm t}/r_{\rm s})^{\alpha}-1$,
which lies between $-3$ and $-4$ for most parameter values. Conversely, Model
B fixes the slope at $r_{\rm s}$ to be $\gamma=-2$ regardless of the other
parameters, but the slope at $r_{\rm t}$ becomes $\gamma(r_{\rm t})=-2(r_{\rm
t}/r_{\rm s})^{\alpha}-1+(r_{\rm t}/r_{\rm s})^{\eta-\beta}$. The
$\beta$-dependent correction term is negligible except for the smallest
$\beta$.

The advantage of Model B is that $r_{\rm s}$ takes on its usual meaning as the
radius where $\gamma=-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-fit
parameters occur for $r_{\rm s}$ in low-$\Gamma$ averaged samples and in some
individual halo profiles. The values of $r_{\rm t}$ and $\beta$ 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-fit 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-fit 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
$\theta$. For orbiting profiles that follow the form $\rho=\rho_{\rm
s}e^{S(r)}$, those derivatives become

$\frac{{\rm d}\ln\rho}{{\rm d}\ln\theta}=\theta\frac{{\rm d}S}{{\rm
d}\theta}\,.$ (24)

The derivative with respect to the normalisation is ${\rm d}\ln\rho/{\rm
d}\ln\rho_{\rm 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.

###  C.1 Orbiting: Einasto

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

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm s}}$
$\displaystyle=2\left(\frac{r}{r_{\rm s}}\right)^{\alpha}$
$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\alpha}$
$\displaystyle=\frac{2}{\alpha}\left(\frac{r}{r_{\rm
s}}\right)^{\alpha}\left[1-\alpha\ln\left(\frac{r}{r_{\rm
s}}\right)\right]\,.$ (25)

If we set $C=2/\alpha$ as in equation (4), the derivative with respect to
$r_{\rm s}$ remains the same but

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\alpha}$
$\displaystyle=\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]\,.$ (26)

###  C.2 Orbiting: New model (Model A)

If $C=0$, the derivatives with respect to $r_{\rm s}$ and $\alpha$ are the
same as for the Einasto profile with $C=0$ (equation C.1), and, analogously,

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm t}}$
$\displaystyle=\left(\frac{r}{r_{\rm t}}\right)^{\beta}$
$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\beta}$
$\displaystyle=\frac{1}{\beta}\left(\frac{r}{r_{\rm
t}}\right)^{\beta}\left[1-\beta\ln\left(\frac{r}{r_{\rm t}}\right)\right]\,.$
(27)

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

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm s}}$
$\displaystyle=2\left(\frac{r}{r_{\rm s}}\right)^{\alpha}+\left(\frac{r_{\rm
s}}{r_{\rm t}}\right)^{\beta}$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm
d}\ln\alpha}$ $\displaystyle=\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]$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln
r_{\rm t}}$ $\displaystyle=\left(\frac{r}{r_{\rm
t}}\right)^{\beta}-\left(\frac{r_{\rm s}}{r_{\rm t}}\right)^{\beta}$
$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\beta}$
$\displaystyle=\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]\,.$ (28)

###  C.3 Orbiting: New model (Model B)

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

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm s}}$
$\displaystyle=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]$ $\displaystyle\frac{{\rm
d}\ln\rho}{{\rm d}\ln\alpha}$
$\displaystyle=\frac{2}{\alpha}\left(\frac{r}{r_{\rm
s}}\right)^{\alpha}\left[1-\alpha\ln\left(\frac{r}{r_{\rm s}}\right)\right]$
$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm t}}$
$\displaystyle=\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}$
$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\beta}$
$\displaystyle=\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)$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\eta}$
$\displaystyle=\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]\,.$ (29)

When setting $C$ as in equation (22), we get

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln r_{\rm s}}$
$\displaystyle=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]$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm
d}\ln\alpha}$ $\displaystyle=\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]$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln
r_{\rm t}}$ $\displaystyle=\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]$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm
d}\ln\beta}$ $\displaystyle=\left(\frac{r}{r_{\rm
t}}\right)^{\beta}\left[\frac{1}{\beta}-\ln\left(\frac{r}{r_{\rm
t}}\right)\right]$ $\displaystyle-\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]$ $\displaystyle\frac{{\rm
d}\ln\rho}{{\rm d}\ln\eta}$ $\displaystyle=\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]\,.$ (30)

###  C.4 Infalling: New model

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

$\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\delta_{\rm 1}}$
$\displaystyle=f_{\rm m}\left[1-\frac{1}{Q(r)}\left(\frac{\delta_{\rm
1}}{\delta_{\rm max}}\right)^{1/\zeta}\right]$ $\displaystyle\frac{{\rm
d}\ln\rho}{{\rm d}\ln s}$ $\displaystyle=-f_{\rm
m}\frac{s}{Q(r)}\left(\frac{r}{R}\right)^{s/\zeta}\ln\left(\frac{r}{R_{\rm
200m}}\right)$ $\displaystyle\frac{{\rm d}\ln\rho}{{\rm d}\ln\delta_{\rm
max}}$ $\displaystyle=f_{\rm m}\frac{1}{Q(r)}\left(\frac{\delta_{\rm
1}}{\delta_{\rm max}}\right)^{1/\zeta}$ $\displaystyle\frac{{\rm
d}\ln\rho}{{\rm d}\ln\zeta}$ $\displaystyle=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]\,.$
(31)

All derivatives share the same $f_{\rm m}\equiv 1-\rho_{\rm m}/\rho(r)$
prefactor, which ensures that they vanish at large radii where $\rho$
asymptotically approaches the mean density. This factor is unity when taking
the logarithmic derivative of $\rho-\rho_{\rm 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 $\zeta=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

$\psi(r)=4\pi
G\left[\frac{M(r)}{r}+\int_{r}^{\infty}\rho(r^{\prime})r^{\prime}{\rm
d}r^{\prime}\right]\,.$ (32)

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

$\sigma^{2}(r)=\frac{1}{\rho(r)}\int_{r}^{\infty}\frac{\rho(r^{\prime})M(r^{\prime})}{r^{\prime
2}}{\rm d}r^{\prime}\,.$ (33)

Analytical expressions for $\sigma$ 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 ${\rm d}\rho/{\rm d}\psi$, 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),

$f(\varepsilon)=\frac{1}{\sqrt{8}\pi^{2}}\int_{r_{\varepsilon}}^{\infty}\frac{g(r)}{\sqrt{\varepsilon-\psi(r)}}{\rm
d}r\,,$ (34)

where $r_{\varepsilon}$ is the radius where
$\psi(r_{\varepsilon})=\varepsilon$ and

$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]\,.$ (35)

The problem of finding $f(\varepsilon)$ has now been reduced to computing
${\rm d}\rho/{\rm d}r$ and ${\rm d}^{2}\rho/{\rm d}r^{2}$. We can generally
write

$\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)\,.$
(36)

For the Einasto profile, we compute

$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}\,.$ (37)

For Model A, the expressions become somewhat more complicated,

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

and

$\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)\,.$
(39)

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,

$\rho(r)=4\pi\int_{0}^{\psi(r)}f(\varepsilon)\sqrt{2[\psi(r)-\varepsilon]}\
{\rm d}\varepsilon\,.$ (40)