From particles to orbits: precise dark matter density profiles using dynamical information

We introduce a new method to calculate dark matter halo density profiles from simulations. Each particle is 'smeared' over its orbit to obtain a dynamical profile that is averaged over a dynamical time, in contrast to the traditional approach of binning particles based on their instantaneous positions. The dynamical and binned profiles are in good agreement, with the dynamical approach showing a significant reduction in Poisson noise in the innermost regions. We find that the inner cusps of the new dynamical profiles continue inward all the way to the softening radius, reproducing the central density profile of higher resolution simulations within the 95$\%$ confidence intervals, for haloes in virial equilibrium. Folding in dynamical information thus provides a new approach to improve the precision of dark matter density profiles at small radii, for minimal computational cost. Our technique makes two key assumptions: that the halo is in equilibrium (phase mixed), and that the potential is spherically symmetric. We discuss why the method is successful despite strong violations of spherical symmetry in the centres of haloes, and explore how substructures disturb equilibrium at large radii.


INTRODUCTION
The observationally inferred density distribution of dark matter in haloes around galaxies offers a crucial hint as to the nature of the elusive substance.However, the observations must be carefully compared with theoretical predictions based largely on numerical simulations (for reviews see e.g.Frenk & White 2012;Vogelsberger et al. 2020;Angulo & Hahn 2022).Dark-matter-only (DMO) simulations have shown that the spherically-averaged density profiles of haloes in the Cold Dark Matter (CDM) paradigm follow approximately the Navarro-Frenk-White (NFW) profile (Dubinski & Carlberg 1991;Navarro et al. 1996bNavarro et al. , 1997;;Dutton & Macciò 2014) described by a divergent cusp ( ∼  −1 ) at small radii, and by a steeper power law ( ∼  −3 ) in the outer regions.The NFW profile has two free parameters which may be fitted to the density structure of simulated haloes for most of the radial extent, but the fit becomes poor in the innermost parts and in the outskirts of the haloes (e.g.Navarro et al. 2004;Diemer & Kravtsov 2014;Fielder et al. 2020;Wang et al. 2020;Lucie-Smith et al. 2022).
Over time, a variety of fitting functions have been proposed to better represent the profile's inner slope, such as Einasto models (Einasto 1965;Chemin et al. 2011) or other forms of double-power law (e.g.Hernquist 1990;Burkert 1995;Zhao 1996;Salucci et al. 2007;Hague & Wilkinson 2013;Oldham & Auger 2016;Hayashi et al. 2020).However, the central regions of the profiles remain ★ E-mail: claudia.muni.21@ucl.ac.uk notoriously difficult to probe due the finite number of particles and consequent need to 'soften' the potential (e.g.Power et al. 2003;Diemand et al. 2004;Dehnen & Read 2011), causing the cusp to be numerically flattened (e.g.Navarro et al. 1996b;Ghigna et al. 2000;Fukushige & Makino 2001;Wang et al. 2020).Constraining the central asymptotic behaviour of the profile therefore remains largely dependent on the number of particles concentrated at small radii.
While the focus in the present work will be on DMO simulations, we note that when baryons are added into simulations, effects such as supernova feedback and enhanced dynamical friction can cause the central cusp to turn into a flattened density 'core' (e.g.Navarro et al. 1996a;Read & Gilmore 2005;Pontzen & Governato 2012;Read et al. 2016;El-Zant et al. 2001;Nipoti & Binney 2014;Popolo & Pace 2016;Orkney et al. 2022).Ultimately, understanding the predicted distribution of dark matter does require such baryonic simulations, especially since there are strong indications of flattened central cores in observations; see e.g.Flores & Primack (1994); de Blok et al. (2001); Marchesini et al. (2002); Battaglia et al. (2008); Walker & Peñarrubia (2011); Oh et al. (2015); Read et al. (2017Read et al. ( , 2019)); Zoutendĳk et al. (2021); De Leo et al. (2023), or for countering views see Pineda et al. (2016); Genina et al. (2017); Oman et al. (2018).The focus in the present work is nonetheless on understanding how DMO predictions can be improved and better understood; we will consider baryonic effects in a future paper.
In the outskirts of haloes, density profiles scatter significantly due to the presence of surrounding substructures and the out-ofequilibrium dynamics of accreting material.For instance, the caustics generated by the infalling particles on their first apocentre passage sets the scale for the splashback radius, which creates an observable signature in the outer regions of halo profiles (Diemer & Kravtsov 2014;Adhikari et al. 2014;More et al. 2015;Shin et al. 2019).Recently, Lucie-Smith et al. (2022) showed that a good fit to the diversity of halo profiles out to two virial radii can be obtained using only three free parameters (i.e., one additional parameter is sufficient to capture the diversity of these outer regions).This relatively simple behaviour may be linked to the typical orbits on which material accretes into a halo, further motivating a study of how the instantaneous profile relates to a dynamically-generated equilibrium profile (e.g.Diemer 2022a,b;Shin & Diemer 2023).
In this work, we present and study a method to calculate dark matter density profiles from simulated halos using dynamical information.This possibility has been discussed before -notably in appendices to Read & Gilmore (2005) and Pontzen & Governato (2013), and in Callingham et al. (2020) -but its possible application to reducing the noise in numerical density estimates has not been explored in detail.Specifically, the technique 'smears' particles in a snapshot along their orbits, spreading the mass of each across multiple density bins.Such a dynamical approach shares some similarities with certain classical mass modelling techniques (Schwarzschild 1979;Syer & Tremaine 1996) but, unlike these, it does not attempt to match observational constraints to underlying orbits and potentials; rather it constructs these from a simulation snapshot.The result is a profile which is averaged over a dynamical time, and which consequently has reduced Poisson noise compared to traditional binned estimates at the same resolution.This, in turn, makes it possible to probe further into the behaviour of the inner regions, at radii where there are very few particles present.
Calculating a density profile through this averaging process inherently assumes an equilibrium, phase-mixed distribution function.This assumption is expected to be significantly broken in the outer parts of a halo approaching the virial radius or beyond.Furthermore, for a practical calculation, we will also assume spherical symmetry (although this assumption could in principle be relaxed).The gravitational potentials of real and simulated haloes are far from being perfectly spherical.Their shapes tend to be closer to triaxial, especially towards the centre (e.g.Frenk et al. 1988;Jing & Suto 2002;Allgood et al. 2006;Orkney et al. 2023); however it has previously been argued using Hamiltonian perturbation theory that approximating the true triaxial potential by a spherically-averaged version should make little difference to dynamical density estimates if the system is in equilibrium (Pontzen et al. 2015).We will return to this point in our discussion.Our results focus on the innermost and the outermost regions of haloes to investigate the limits of dynamical halo profiles subject to these coupled assumptions of equilibrium and spherical symmetry.
The rest of the paper is structured as follows.In Section §2, we explain the procedure used to generate the dynamical density profiles.In Section §3, we describe the simulation suites and the selection of snapshots analysed in this work.In Section §4, we present the main results for the dynamical profiles, focusing on the inner and outer regions, and comparing our dynamical technique to traditional binned methods.In Section §5, we discuss the implication of our results and outline possible further work.

METHODS
We now describe the methods used to construct dynamical profiles.Section §2.1 considers the construction of a spherically-averaged gravitational potential starting from a simulation snapshot; the calculation of particle orbits within that potential; and finally the computation of the dynamical density profile.In Section §2.2, we introduce a refinement to the method which improves the accuracy of the orbit integration around apocentre and pericentre.Then, in Section §2.3, we describe an iterative process via which a self-consistent densitypotential pair may be generated.

Creating the dynamical density profiles
We start by assuming that we have a snapshot containing only dark matter particles, centred on the target halo.The spherically-averaged gravitational potential given by all the particles in the snapshot is then calculated in bins of width Δ according to the discretized integral where  is an index over the bins,  is the bin number for which the potential is being calculated, and   is the radius in the centre of the th bin, taking the value   = (  − 1/2)Δ.In addition,  (<   ) is the mass enclosed within radius   , and  is the gravitational constant.Although the potential for each bin  is evaluated from quantities at the centre of the bin, the values are assigned to the right edge of the corresponding bins, since Φ(  ) represents the average of the potential over the entire bin .The zero point of the potential is set at  = 0 (the left edge of the first bin).
Equation ( 1) is the simplest of several possible choices to perform numerical integration.We tested that adopting a more sophisticated method does not significantly affect the final results.Therefore, we adopted the simple approach for transparency.
The total number of bins over which Φ is calculated is determined by the radius of a 'bounding sphere' centred around the halo.In addition to choosing the radius at which to truncate the potential, we must also decide how to treat particles whose orbits cross this boundary.In keeping with the core assumption of equilibrium, we make the boundary reflecting, i.e. particles bounce elastically off it.
One may equivalently imagine the potential as having an infinite potential step at the truncation radius.While this is unphysical for any individual particle considered in isolation, across the population it is equivalent to the much more reasonable assumption that the outwards flux through the sphere is balanced by a matching inwards flux.This assumption can be tested by changing the truncation radius; the halo virial radius is a natural first choice, and we will explore the effects of other choices on the final density profile in Section §4.2.2.
Assuming equilibrium, the probability density   () of finding particle  at radius  is proportional to the time spent by the particle in the infinitesimal interval around that radius: where   () describes the radius as a function of time for the particle (on its spherical idealised orbit),   is the period of the orbit,   is its specific energy, and   is its specific angular momentum.Rather than calculate   directly we first calculate an unnormalised version of the probability,  , ≡ (  /2)   (  ).Here  indexes the particles, and  indexes the spatial bins.By writing the specific energy of a particle as the sum of the potential energy, the kinetic energy due to the angular momentum, and the kinetic energy due to the radial motion, we can solve for  and obtain (3)  3) evaluated for a typical particle (light-blue bins), with a bin size Δ =  /2, compared with the analytic integrand (black line).The integrand is well behaved for most of the radial range of the orbit, and therefore well approximated by the binned density.However, it has two integrable divergences at pericentre and apocentre (here located at ∼ 2.2 kpc and ∼ 8.7 kpc, respectively).Even if the particle never reaches the centre of one of these extremal bins, it may still spend significant time within the bin.Capturing this effect correctly in the binned probability requires the special treatment explained in the text.The dark-blue shaded areas represent the analytical corrections added at the pericentre and apocentre for this orbit.
Note that this expression is only valid between pericentre and apocentre; outside this radial range, it becomes imaginary.However the true probability of finding the particle outside the extrema of its orbit is zero by definition, and therefore one may make Eq.( 3) true for all radial bins by taking its real part.We produce a normalized probability for each bin  and particle  according to where ℜ denotes the real part.If a particle  is on an almost perfectly circular trajectory, it may remain within a single radial bin  for its entire orbit; in this case, the equation above fails and instead a unit probability is assigned to the bin enclosing the original position of the particle in the snapshot,  , = 1.The density at the centre of bin  can then be estimated from the set of  , as where   is the mass of each particle , and there are  particles in total.
The statistical errors in the dynamical density profile are estimated using bootstrapping.For each of 100 bootstrap samples, we create a mock set of particles by sampling (with replacement) from the actual set of particles in the halo; we then perform the full dynamical density estimate on the mock set of particles.We determined that 100 bootstrap samples was sufficient to achieve convergence on the 95% confidence interval; in Section §4, our results are shown with these uncertainties as a shaded band.

Improving accuracy at apocentre and pericentre
The function in Eq. (3) has two integrable divergences located at the pericentre and apocentre of each orbit (Figure 1).Unless the bins are infinitesimally small, the probability of finding the particle in the bin  , containing such a divergence might be misestimated.To correct for this, in these two bins we use an approximation scheme based on a local Taylor expansion of the potential.We define the effective potential as Φ eff = Φ +  2 /(2 2 ), and expand Φ eff ( 0 + ) around  0 , where  0 is the divergence point (pericentre or apocentre) for every orbit, i.e., a root of Eq. (3).We now consider the case of a pericentre where the divergence  0 is inside the th bin (i.e., ( − 1)Δ <  0 < Δ), as an example.The mean value of ℜ  −1 across the entire bin may be calculated as q ≡ 1 Δ Here we have also used the fact that Φ eff ( 0 ) =   , by definition.We can furthermore approximate dΦ eff /d |  0 ≈ dΦ eff /d |   to avoid having to calculate the exact location of the divergences; this will give us a correction that is accurate to first order.The integration is then analytically tractable, giving This analytical estimate of the mean value is then used to represent the value of the probability density function within the pericentre bin   .The apocentre bin is treated in the same way, and both corrections are included before producing the normalized probability according to Eq. (4).There are two cases in which these corrections cannot be evaluated.One of them is when an orbit is unresolved (i.e. its probability function only spans one bin), since in that case pericentre and apocentre are coincident.As previously stated, when this occurs, the particle is given unit probability to be found within the single bin, and corrections are not required.The apocentre corrections are also ignored when the particle's apocentre falls outside of the radius of the 'reflecting wall' which serves as the boundary for the halo.Since the particles can be thought of as being reflected back once they hit the boundary, their radial paths are truncated at the location of the wall, and no apocentre corrections are required.

Iterating the potential
The dynamical density profile given by Eq. ( 5) implies also a mass profile  (< ) and therefore a potential Φ() through Eq. (1).However, the potential used in producing the density estimate was initialized directly using the particle radii from the original snapshot.The overall procedure, therefore, results in an inconsistent potentialdensity pair.The difference between the mass distribution is especially evident in the inner regions because our potential is calculated without softening, and the pericentres of orbits can therefore reach radii closer to the centre of the haloes.To resolve this discrepancy, we iterate until a self-consistent density-potential configuration pair is reached.Over the course of the iterations, the gravitational potential from the simulation is gradually transformed into the potential inferred from the dynamical density profile.This technique also removes any discontinuities in the derivatives of the potential at small radii due to the finite particle number.
The iteration process involves a series of steps: (i) A dynamical density profile is first obtained as described in Sections §2.1 and §2.2.
(ii) The mass distribution implied by the dynamical profile is calculated according to The mass at the centre of the bin (  ) is then obtained by averaging the mass at adjacent edges.
(iii) The new mass distribution is inserted into Eq.( 1) to evaluate a new gravitational potential.
(iv) The angular momenta of the particles is assumed to be unchanged, and the energies are updated by keeping the radial action constant at first order (see below).
(v) The cycle is repeated, starting from point (ii) and using the updated dynamical profile, until convergence in the dynamical profile is reached.
Evolving the gravitational potential into the new configuration will affect the phase space distribution of the particles.Hence we require the energies of the particles to change accordingly.In step (iv) the updated energies are calculated by keeping the radial action of each particle constant to first order, for each particle , where  new, , Φ new and  old, , Φ old are the specific energy and the potential after and before the iteration respectively, and ΔΦ = Φ new − Φ old .We keep   constant since we can interpret each iteration as a small change to the potential of the halo, akin to an adiabatic relaxation.This process does not correspond to a literal physical evolution of the halo in time, but an adiabatic transformation is nonetheless the most conservative way to map orbits from the potential at each iteration to the next.In other words, we assume that the action distribution of the particles in the simulation is sampled from an underlying 'true' distribution (as would be attained by a simulation of infinite resolution).We then recover the potential implied by the dynamical profile given this action distribution.If we assume that the change to the potential between iterations is sufficiently small, we only need update the actions at first order in the potential change, i.e.Eq.( 9).The definition of the radial action is With this in hand, we solve Eq. ( 9) to first order in the quantities ΔΦ and Δ  =  new, −  old, .By Taylor expanding, we find (11) i.e. the change in energy is equal to the average of the change in potential, weighted by the probability of finding the particle at a given radius.(At first order, the changes to the values of apocentre and pericentre of the orbit do not contribute to Δ, and can therefore be neglected.) The first iteration produces a significant change in the inner density distribution but after approximately 3 iterations, convergence in the dynamical profile is reached (i.e. the changes in the density profiles become significantly smaller than the bootstrap-determined uncertainties).We will discuss this further in Section §4.1.2below.

THE SIMULATION SNAPSHOTS
We analyse a selection of seven snapshots drawn from cosmological zoom simulations of dark matter haloes spanning a wide range of masses, from ∼10 9 M ⊙ to ∼10 12 M ⊙ (see Table 1).
The five smallest haloes are part of the Engineering Dwarfs at Galaxy Formation's Edge (EDGE) project (Agertz et al. 2019;Rey et al. 2019Rey et al. , 2020;;Orkney et al. 2021); the two largest haloes were taken from the vintergatan-gm project, which in turn uses the initial conditions described by Rey & Starkenburg (2021) As previously stated, we consider the dark-matter-only simulations from these suites, i.e. they do not contain any baryonic components; hence steep cusps are expected in the central regions of the density profiles.
The selected haloes were re-simulated at two different resolutions; the particle mass ratio between the lower and in the higher resolution runs is 64 (for EDGE) and 8 (for vintergatan-gm).Both suites of simulations are generated using the adaptive mesh refinement (AMR) code ramses (Teyssier 2002).The mesh is refined whenever a grid cell contains more than 8 particles; consequently, the softening lengths are adaptive and we provide a softening scale estimate  equal to the size of the smallest grid cell used for gravity calculations.We call low resolution the simulations with softening scale of 0.095 kpc (0.142 kpc for the vintergatan-gm haloes), and high resolution the ones with softening of 0.012 kpc (0.035 kpc for the vintergatangm haloes).Ultra-high resolution runs with softening scale ∼ 0.006 kpc are also available for some EDGE simulations.All the snapshots analysed in the current work are taken at the present day ( = 0).
Simulation snapshots are loaded using pynbody (Pontzen et al. 2013).Before processing, each halo is centred using the shrinkingsphere method of Power et al. (2003); the central 1 kpc is used to calculate a centre of mass velocity, which is then subtracted from all particles.We also calculate a virial radius,  vir , defined to be the radius at which the enclosed mean density is equal to 178 times the cosmic mean.
All particles interior to the reflecting wall at the time of the snapshot are included in the calculations.Some of the selected haloes contain large substructures, especially in their outskirts; these are deliberately retained in our analysis in order to test the limits of the assumption of equilibrium.The reflecting boundary described in Section §2.1 was placed at 120 kpc for the haloes with mass ≲ 5 × 10 9 M ⊙ .This is between 2 and 3 times the size of their virial radii, a choice which allows us to explore how the dynamical information affects the density distribution in their outer regions.The boundary for the two largest haloes was placed at 350 kpc, which is approximately the location of their virial radii, and was not extended to larger radii in this work because the 'zoomed' region of these haloes is only twice the virial radius, beyond which low resolution particles are present.For efficiency, the dynamical profiles of the two largest haloes are generated using only a randomly selected fraction (a third) of the particles.
While it is not possible to recreate precisely the in-simulation softening  with a spherical approximation, it is clear that the bin width Δ must be comparable to  in order that the potential is meaningful.We found that our results were insensitive to the precise bin width chosen, provided that it is of this order, and therefore chose to fix Δ = /2.This choice of bin width is sufficiently small to allow investigation of the dynamically-inferred density profile close to the halo centre.We note that for  ≲ 3 ≡  conv the effect of spurious relaxation in simulation becomes important and a profile constructed through direct particle binning is poorly resolved.Detailed studies of convergence (e.g.Power et al. 2003;Gao et al. 2012;Ludlow et al. 2019) show that the value of  conv must be determined empirically for each simulation setup, and any relation to the softening length  is approximate; the scale is mainly dictated by the number of particles present in the innermost regions.Our comparisons of binned profiles between high and low resolution simulations below confirm that  conv ∼ 3 gives a sufficiently good approximation to the innermost reliable radius of the low resolution binned profiles 1 .

RESULTS
In this Section, we present and discuss the dynamical density profiles of our dark matter haloes.In each case, we calculate dynamical profiles from the low resolution snapshots and compare them with binned profiles from both low resolution and high resolution snapshots.The profiles are shown in Figures 2, 3 and 4 (for lowest-mass dwarf, intermediate-mass dwarf and Milky-Way-mass haloes respectively), alongside images of the haloes' dark matter density projected down the  axis.We compare our dynamical profiles (blue lines) to the traditional binned estimates from both the high and the low resolution snapshots (black and pink points respectively), which are 1 Using Eq.( 14) in Ludlow et al. (2019) to calculate the convergence radius, with the constant  03 calibrated using a subset of our haloes at different resolutions, we obtained estimates between 2.6 and 4.5 across all snapshots.
Our  conv = 3 estimate, therefore, falls within this range.
plotted down to their estimated softening length (see Table 1).Inset panels show the inner density profile in greater detail.
Overall, the dynamical profiles (blue lines), obtained from the low resolution simulations, agree well with the low resolution binned profiles (pink points) for the majority of the radial extent of the haloes.The 95% bootstrap-determined uncertainties on the dynamical profiles are shown as shaded blue bands, and are significantly smaller than the 95% Poisson noise on direct binned estimates at the same resolution (pink error-bars).This follows from the fact that the particles in the original snapshot are now spread across multiple density bins, hence providing better statistics.
By dividing the total volume occupied by each halo into thin shells, we can also calculate the average radial velocities of the particles contained within the shells.These are shown for the low resolution simulations in the panels below the density profiles in Figures 2 -4.These values will help us discuss below how well the assumption of equilibrium holds for each halo.
We will first discuss the behaviour of the dynamical profiles in the inner regions (around or even interior to the traditional convergence radius; Section §4.1), then in the outer regions (around and beyond the virial radius; Section §4.2).

Inner regions
The direct comparison of dynamical profiles (blue lines) with binned profiles from higher resolution simulations (black points) is of considerable interest: it addresses the question of whether our technique can partially correct for finite particle number in the innermost regions of the halo.
At radii below the approximate convergence radius of the low resolution binned profiles ( conv = 3, indicated by the pink arrows in Figures 2, 3, and 4), our dynamical density cusps are steeper than the traditional binned profiles at the same resolution.This is particularly clear in the case of the Milky-Way-mass halos (Figure 4).Comparing our results to the binned distribution of the high resolution simulations (black points), we see that the dynamical method is, in nearly all cases, able to predict the 'cuspier' behaviour of higher resolution simulations below  conv .This is especially evident in the larger haloes due to the smaller Poisson noise in the central regions.
Halo600 is an exception in which the dynamically predicted density is substantially lower than that in the high resolution simulation; Section §4.1.1 considers that case in some detail, and more broadly discusses caveats about making comparisons between low and high resolution simulations.Nonetheless, in the other cases studied, the dynamically predicted cusp extends below  conv of the low resolution simulations, where very few particles are present at the time of the snapshot 2 .As well as being less biased than the binned profiles, our dynamical profiles also have lower numerical noise.On 2 For the EDGE haloes there are only 150-250 particles below  conv , and 350-450 for the vintergatan-gm haloes; this is 0.38% and 0.007% of the average across all halos, the uncertainties at small radii (between  and  conv ) are reduced by a factor of 12 compared to traditional binned estimates.Thus, our technique uses information about the entire phase-space of the particles to produce more precise central density profiles which partially correct for the effects of softening and which are less subject to Poisson noise.
Poisson noise could also be mitigated by stacking binned profiles from adjacent snapshots (similarly to the procedure outlined in Vasiliev 2014).Figure 5 shows an example of the binned profile obtained by stacking 6 adjacent snapshots of Halo1459.This is compared to our dynamical density profile (blue line) and to the binned profile obtained from a single snapshot at  = 0 (pink points).number of particles enclosed by the virial radius of the two simulation suites, respectively.Similarly to the other cases, the dynamical density profiles from the low resolution snapshots agree well with both binned profiles.Halo600 is an outlier since it recently had a merger close to the halo's centre which disrupted the equilibrium in the inner regions; as a result the plot of v /v vir shows significant deviations from zero at small radii.Halo624 has a large substructure within its virial radius which will reach the centre of the main halo and merge with it in the next ∼ 500 Myrs.(The structure is found slightly closer to the centre in the high resolution simulation.)The significant disruption caused by this substructure to the halo's equilibrium is also evident in the average radial velocity panel, but our dynamical method nonetheless recovers a sensible 'smoothed' density profile.Similarly to the other haloes, the dynamical density profiles from the low resolution snapshots agree well with both the low and high resolution binned profiles.For efficiency, the dynamical profiles for these haloes were generated using only a randomly selected fraction (a third) of all the particles within the halo and therefore even smaller errors on the dynamical density profile are achievable in principle.In these examples, all substructures are small (less than 1% of the halo mass) and do not have a visible effect on the density profiles.
Stacking the profiles results in considerable reduction in shot noise, similar to the effect observed in the dynamical profile.However, the method fails to reproduce the steeper central gradient observed in the dynamical profile below the convergence radius, which implies a significant disagreement with the binned profile from the high resolution runs.This is due to the fact that the stacked profile retains the effects of gravitational softening and of relaxation caused by encounters between the particles in the low resolution simulation as it evolves over time.By contrast, in the dynamical profile no softening is used and the orbits are integrated independently of each allowing the iteration process (Section §2.3) to correctly recover a steep central cusp.While in principle the stacked profile could also be iterated by combining it with our dynamical method, this would entail significant complexity due to the starting potentials in each snapshot differing from each other, as well as from the final combined potential.We therefore leave any investigation of such a combined stacked-dynamical profile to future work.At radii just larger than  conv , we notice a small but statistically significant density excess in both the binned and dynamical low resolution profiles when compared with the high resolution binned profiles.This excess only covers a few density bins and is more evident for some haloes (e.g.Halo605 and 624) than others; see the inset panels zoomed in on this radius in Figure 3. Since this feature is also present when using binned methods, it must be unrelated to the inclusion of dynamical information into the calculations.We therefore leave investigation to a future study.

The challenge of direct comparisons between differing resolutions
Overall, the improvement offered by dynamical profiles over binned profiles is significant: the uncertainties at small radii are significantly mitigated compared to binned estimates, making it a substantially more precise technique.Qualitatively, it is clear that the dynamical profiles reproduce steeper profiles which appear to be in agreement with higher resolution simulations within the 95% error bounds.However, quantifying how accurate the dynamical estimates are compared to the true density distributions (i.e. the density pro- files that would be obtained from simulations of infinite resolution) is difficult for two reasons.The first is the problem of formulating a suitable comparison summary statistic; the second is the impact of small differences in halo formation and merger history on the final profile.We will describe each of these in turn.
The most natural way to measure the accuracy of a low resolution density profile would be to construct a chi-squared test to decide whether the binned or dynamical profiles more accurately predict the high resolution result.However, the size of the statistical errors on the dynamical profile are substantially smaller than those on the binned profile, putting the dynamical profiles at an automatic disadvantage in such a test.Even if one were to artificially inflate the dynamical profile error estimates, the results would remain very sensitive to the precise radial range over which the statistic is calculated.The dynamical profiles clearly predict more accurate densities interior to  conv , but outside this radius the situation is more nuanced.In particular, at large radii, the dynamical profiles' tendency to wash out substructure would lead to a heavy  2 penalty (as will be discussed in Section §4.2 below).There is therefore no straight-forward quantitative measurement of the improvement offered by dynamical density profiles, despite the clear qualitative advantages in the cusp region.
The second challenge relates to recent events in the formation and merger history, and is most clearly seen in the case of Halo600 (shown at the top of Figure 3).As with the other examples, the gradient of the dynamical profile interior to  conv is steeper than the low resolution binned profile; however, unlike the other cases, the steepening in Halo600 is insufficient to reach agreement with the high resolution binned profile.The reason can be traced to the halo's recent history in the respective simulations.The low resolution version of Halo600 underwent a minor merger at  = 0.03 (∼ 70 Myrs before present day).This merger only occurred in the low resolution version of the simulation.Although the mass of the merger is relatively small (∼10 8 M ⊙ , around 2% of the total host mass), its centre of mass before disruption is located within 1 kpc of the centre of mass of the main halo.By tracking the particles that formed the subhalo to  = 0, we find that they have traversed the halo from one side to the other, and remain in disequilibrium.The out-of-equilibrium behaviour is also Dynamical density (after iterations) Dynamical density (before iterations) Binned profile (high res) Figure 6.Dynamical density profile multiplied by  2 before (yellow) and after (blue) the dynamical iteration process compared to the high resolution binned profile (black points), shown here for the example of Halo1459.The pink arrow marks the convergence radius of the low resolution simulation binned profile (which, for clarity, is not itself shown).The effect of the iterations is especially evident at small radii, where they act to make the central regions moderately denser, in better agreement with the high resolution profile.
visible as large fluctuations in the binned radial velocities as seen in the lower panel of the Halo600 plot in Figure 3.Despite this, note that the dynamical density profile still performs somewhat better than the binned profile.
From the above analysis, we deduce that even a relatively small merger might affect the equilibrium of a halo.A statistical study on a larger sample of haloes is necessary to constrain the exact relationship between merger-to-main halo mass ratio and the effect that the merger events have on the dynamical profile.Other features will also play a role, such as the object infall velocity or the angle of collision.The investigation of these effects is beyond the scope of this work.

Effect of potential iterations
Having established that dynamical profiles offer an accuracy improvement over binned profiles near the centres of halos, albeit one that is hard to quantify, we now consider the effect of the iterative part of our algorithm (Section §2.3) in achieving this.
Figure 6 shows the effect that the iteration process outlined has on the dynamical profile.After the iterations, the profile's central gradient becomes moderately steeper.This can be understood by considering that the particles previously located at larger radii are now allowed to extend further inwards compared to their original positions in the snapshot, hence increasing the density in the inner regions.Note that the increase in central density may appear to violate mass conservation, since the total mass of the halo should be unaltered.However we verified that the mass enclosed converges to the same value at the virial radius; the volume of the sphere inside  conv is just 3 × 10 −5 % of the total volume inside the virial radius, and therefore a very small reduction in density across a large range of radii is able to provide the mass for an increased density cusp.
Overall, we therefore conclude that the iterative component of the algorithm is important not just for self-consistency (as argued in Section §2.3) but also to achieve the increased densities interior to the binned profile's convergence radius.Given that we kept actions fixed (to first order) during the iterations, one can envisage them as adiabatically transforming away some numerical effects of softening. .Dynamical density profile multiplied by  2 (blue line) obtained from the high resolution simulation of Halo1459, compared to the binned density profiles of the high (black points) and ultra-high (green points) resolution snapshots.The binned profile obtained from the low resolution snapshot is shown for reference (pink points).The black arrow indicates the approximate convergence radius of the high resolution binned profile (3 ).The dynamical density profile from the high resolution simulation predicts the ultra-high resolution simulation well, underscoring how the method can be applied at any resolution to extract additional information.

Comparison at ultra-high resolution
So far, we have applied our dynamical method to the low resolution snapshots and compared our results against the binned profiles obtained from the high resolution versions of the simulations.In order to understand whether this improvement is independent of resolution, we now test the dynamical approach on the high resolution simulations and compare the results to ultra-high resolution snapshots.
Figure 7 shows the dynamical density profile calculated from the high resolution simulation of Halo1459 compared to the binned distribution from an ultra-high resolution simulation with  ≃ 6 pc (half the softening length of the high resolution snapshots previously analysed).We take Halo1459 as an example, but similar results are observed for the other haloes.
All the conclusions drawn in the case of the low resolution dynamical profile are still valid when the code is applied to the high resolution snapshot: the dynamical density shows smaller uncertainties, a steeper cusp that extends further inwards and approximately follows the higher resolution binned profile, and a small density excess at ∼ conv in the lower resolution profile.Overall, this confirms that the improvements obtained by adding dynamical information to the profiles continue even for increasingly precise simulations, making them resolution-independent.
In Figure 8 we show the dynamical profile obtained from the high resolution simulation of Halo600.When the dynamical code was previously applied to the low resolution simulation (top of Figure 3), we saw that the steepening in the cusp was insufficient to reach agreement with the high resolution binned profile.This is not the case when the dynamical profile is calculated from the high resolution snapshot: the cusp of the dynamical profile is entirely consistent with the ultra-high resolution binned profile.This provides further evidence that the disagreement between the dynamical and binned profiles at small radii in the low resolution case is a result of the disequilibrium caused by the merger event, which did not occur in the high resolution version.7 but for Halo600.The dynamical profile from the high resolution simulation of this halo shows a steep cusp consistent with the ultra-high resolution binned profile.The high resolution simulation, unlike the low resolution version, did not recently undergo a merger close to the halo's centre.This provides further evidence that the disagreement between the dynamical and binned profiles seen at small radii in the low resolution case is due to disequilibrium caused by the merger event.

Outer regions
Having shown that the dynamical profile technique performs well in suppressing numerical noise at small radii (comparable to the convergence radius), we next consider its predictions at large radii (comparable to the virial radius  vir ).At such large radii, finite particle number is unlikely to be a limiting factor in drawing physical conclusions and therefore the motivation for studying the dynamical profile is different.Specifically, we are interested in understanding the degree to which halos may be considered equilibrium structures; departure from such equilibrium invalidates our assumptions and therefore should lead to an inaccurate profile.The virial radius roughly defines the point past which most particles are no longer gravitationally bound to the halo, such that infalling particles from the halo's environment begin to dominate.
We are able to study the dynamical profiles beyond  vir for dwarfscale haloes, since the zoom region extends several times further out.Beyond the virial radius we find, as expected, that the dynamical profiles are typically inaccurate; see Halo1445 and 1459 in Figure 2 for particularly clear cases.
This provides one clear signature of out-of-equilibrium dynamics.However, another way to measure departures from equilibrium is via the binned average radial velocities of the particles ( v ), which should be consistent with zero in equilibrium.Measured values of v are shown in the panels below the density profiles in Figures 2,  3, and 4. As expected, these values deviate strongly from zero outside the virial radius, confirming our interpretation above.However, more surprisingly, the mean velocity values deviate from zero even interior to the virial radius, in regions where the binned and dynamical profiles fully agree (e.g. in Halo600, 605, 1459 over the radial range 1 <  < 40 kpc).The root-mean square deviation of the radial velocities of all haloes (excluding Halo624) in the region  <  vir is of order ∼ 5% of the virial velocity.These deviations are statistically important, and yet do not appear to have a significant effect on the overall density structure which is in good agreement with the binned estimates.This suggests that the dynamical profiles are robust to even significant violations of their equilibrium assumption.

The role of substructures
Although dynamical profiles remain robust despite the existence of smooth inflows detectable well interior to the virial radius, a more difficult challenge is posed by substructures.Most haloes have spikes in the binned density distribution at certain radii: for Halo600, 1445, and 1459 (Figure 2, and top of Figure 3) these can be seen beyond the virial radius at  ∼ 90 − 100 kpc, while for Halo624 (bottom of Figure 3) we see them much closer to the centre at  ∼ 10 − 20 kpc.We refer to the locations of these features as  spike .We verified that these local density spikes are indeed caused by substructures (see brown circles in the haloes density images in Figures 2, 3, 4), which each contain between 3% and 9% of the mass of the main halo.All the other substructures present within the reflecting boundary have masses below 0.5% of the main halo's mass.
The dynamical density profile does not reproduce spikes associated with substructure; by design, it smears them out along their orbit without taking into account the self-binding of the substructure.This leads to systematic differences between the binned and dynamical profiles, since the spike is smoothed out while conserving the total mass.This effect is especially evident outside the virial radius in Halo1445 and 1459 (Figure 2).In these cases, substructures (indicated by brown arrows at the appropriate radii on the density plots) coincide with significant disagreements between binned and dynamical halo profiles.
Halo624 contains a large substructure of mass ∼1.4 × 10 9 M ⊙ within its virial radius (at ∼20-25 kpc).This is clearly visible in the density image at the bottom of Figure 3.The substructure will reach the centre of the main halo and merge in the next ∼ 500 Myrs (based on its estimated infall velocity at  = 0), and the disruption to the halo's equilibrium caused by the presence of substructure is also evident in the large deviations from zero in the average radial velocity panel.Despite this, the dynamical profile still faithfully represents the density distribution at radii between the centre of the halo and the location of the substructure.This shows that the effects of the dark matter spike are localised to the area around the substructure, and our method can represent the correct density distribution in other regions of the halo.
Halo605 provides an example with no large substructures present within the entire volume analysed.Despite fluctuations of the binned mean velocity, the dynamical profile agrees with the binned profile up to radii of 100 kpc which is around 2 vir .Taken with the discussion above, this counterexample strongly suggests that substructures, rather than smooth radial flows, are the dominant factor in determining whether binned and dynamical profiles differ significantly, and that the effect of substructures on the profile is always localised.

Effect of the reflecting boundary
As described in Section §2.1, the dynamical density profile requires an outer boundary condition.We have assumed a perfectly reflecting wall, which is equivalent to assuming that the particles flowing inwards across the boundary are exactly balanced by the flux outwards, in keeping with our broader assumption of dynamical equilibrium.However, there remains the freedom to move the reflecting wall to an arbitrary location.We carried out a number of experiments to determine the effect of this choice.If, for example, a boundary is placed inside the virial radius we found that the dynamical density profile is insensitive to the particular choice of location.However, in order to probe the outer parts of the halo the results above were all presented with the boundary outside the virial radius.In this case, there is more sensitivity to the particular choice of location.
An example is shown in Figure 9 for Halo605.As usual, the binned profile is shown by pink points with error bars while dynamical profiles are represented as lines.Here, however, we show two alternative dynamical profiles: one with the reflecting boundary moved inwards to 100 kpc (≃ 2 times the virial radius, as previously adopted, and illustrated here with a blue line) and one with the reflecting boundary moved outwards to 200 kpc (≃ 4 times the virial radius, illustrated with a grey line).This shift causes the dynamical profile to deviate from the binned density in the range  vir <  < 2 vir , where there was previously agreement.
The change is caused by particles that, at the time of the snapshot, are exterior to 2 vir but infalling, such that they spread to lower radii when the equilibrium assumption is imposed.The binned profile shows a 'kink' at ≃ 100 kpc which means that, in this particular case, there is a relatively large mass in such infalling particles.When the reflecting wall is located at 2 vir , these particles are safely isolated outside of the boundary, and therefore cannot affect the density profile.
In a sense, moving the reflecting wall to increasingly large radii provides a prediction of the future profile, since it extrapolates to a time when far-out particles have been able to fall into the inner regions.However, we did not study to what extent this can actually be used to make meaningful predictions and we caution that the actual process via which infalling particles relax into virial equilibrium is unlikely to be fully captured; in effect, our algorithm assumes conservation of their adiabatic invariants which is unlikely to be correct in detail.
For practical purposes, the most conservative choice of reflecting wall boundary is at the virial radius, but our results show that it is entirely possible to obtain accurate profiles out to twice the virial radius.Beyond this, dynamical profiles with extended radial range may be of interest for understanding the accretion processes of halos and 'splashback' features (Diemer & Kravtsov 2014;Adhikari et al. 2014;More et al. 2015;Shin & Diemer 2023;Lucie-Smith et al. 2022), something we will investigate in the future.

CONCLUSIONS AND DISCUSSION
We presented a new method to estimate spherically-averaged densities in cosmological dark matter haloes.Instead of binning the particle in a snapshot by radius, which is the most obvious and prevalent approach, we use the velocity information in the snapshot to 'smear' each particle along a trajectory, substantially reducing Poisson noise.Such a method has been proposed before (Read & Gilmore 2005;Pontzen & Governato 2013), but our work is the first systematic investigation of the approach.Additionally, we derive new corrections to take into account the integrable singularities at apocentre and pericentre, and introduce an iterative process to obtain a self-consistent potential-density pair.After iteration, we obtain central density estimates which (except in one case, Halo600, where a recent merger has occurred) follow the trend set by higher-resolution simulations.The agreement persists interior to the binned profile convergence radius, and all the way down to the simulation softening length.This highlights how our technique can squeeze extra information about the central regions of halos from existing simulations.
In the outer regions, the dynamical profiles continue to agree with the binned profiles even out to several times the virial radius, provided that no substructures are present.If substructures are present, the assumption of equilibrium is locally broken and the profiles in the vicinity of the substructure are 'smoothed' relative to the binned profiles.Nonetheless, the overall profiles remain accurate.Eventually, at  3) when the reflecting boundary is placed at 100 kpc (blue line) and then moved to 200 kpc (grey line) compared to the low resolution binned profile (pink points).The dynamical profile agrees well with the binned one when the boundary is placed anywhere up to 100 kpc, around twice the virial radius, but differs once contributions from particles out to 200 kpc are included in the calculations.These discrepancies propagate inwards to smaller radii, even below the virial radius (55 kpc, indicated by the vertical dashed line).This behaviour reflects our algorithm's extrapolation of how particles and substructures in the outskirts, while currently unbound, will ultimately fall into the halo at later times in the simulation, altering the density distribution.
approximately  ∼ 4 vir , effects from the haloes' environments start to dominate, bringing the haloes too far out of equilibrium for the dynamical profiles to give meaningful density estimates.Including particles from these distant halo outskirts can produce changes to the dynamical profiles, sometimes even at radii below the virial radius.This is not a surprising result since the particles at large radii will eventually fall into the halo at future times in the simulation, and the dynamical approach is extrapolating the orbits of these particles accordingly.However, whether the resulting profile can be considered a 'prediction' of the growth of the dark matter distribution at later times remains to be investigated.
These effects in the outer parts of the halo relate to the departure from perfect equilibrium (or phase-mixing), which is one of two key assumptions underlying the method.The second assumption is that the potential is spherically symmetric; this assumption is, in fact, broken by all our simulated halos, since they have triaxial equipotential surfaces.The fact that the dynamical profiles are accurate despite this broken assumption warrants further discussion.
? estimated the shapes of the five least massive dark matter haloes studied in this work by calculating the intermediate-to-major and minor-to-major axial ratios (/ and /) up to approximately 20 kpc in radius.The exact shape of each halo is not constant with radius: the / ratio for all the haloes varies within the interval 0.4-0.8(ratios of exactly 1 indicate perfect sphericity).The DMO haloes are generally the least spherical near their centre, becoming increasingly more spheroidal at radii beyond the cusp (≳ 1 kpc).Nevertheless, the dynamical density profiles are able to correctly represent the density distributions for the entire radial extent of the haloes.
The nature of the particles' orbits in an aspherical system is very different from the orbits that would be observed in a sphericallyaveraged version of the same potential.In the spherical case the angular momentum of individual particles is always constant; this is not the case in aspherical systems where only the total angular momentum of the entire system is conserved.This allows specific types of orbits (which would not be allowed in a spherical potential) to exist, such as box orbits which plunge through the centre of the halo.Therefore, the fact that we are able to infer reliable results about the haloes' properties using only an artificial version of the dynamics which does not correspond to the real trajectories of the particles is not a straightforward outcome.
However, such an outcome was previously predicted by relying on having a distribution function of particles in equilibrium (Pontzen et al. 2015).For every particle that is on an orbit losing angular momentum, there must be another particle on an orbit gaining angular momentum.To put it another way, the net flux of particles through the spherical action space must be everywhere zero, and so in a statistical sense, averaged across all particles, the spherical orbits remain a good approximation.For a more technical discussion, see Pontzen et al. (2015).The present work provides additional evidence that this mapping from a real triaxial system onto an effective spherical system is able to give accurate insights into dark matter halo structure.That said, the dynamical density method could be readily extended beyond the assumption of spherical symmetry, similarly to other mass modelling techniques (Schwarzschild 1979;Syer & Tremaine 1996).
Overall, our dynamical method for the evaluation of dark matter density profiles is a powerful tool which can represent the correct mass distribution even when its fundamental assumptions are partially broken, making it largely applicable to a wide range of systems.
However, dark matter halos in the real universe have potentially been altered by baryonic effects, something which we have not investigated at all in the present paper.In forthcoming work, we will apply our dynamical density code to hydrodynamical simulations.Adding baryons to the simulations will likely alter the shape of the profile's inner regions, transforming the cusp into a flatter core.At a technical level, the gravitational potential can no longer be made fully self-consistent with the dark matter density distribution, and the potential will need to be evaluated directly from the snapshot for the baryonic component.The iterative procedure that we have outlined will therefore need to be refined before we can use it in such cases.

Figure 1 .
Figure 1.The binned probability density implied by Eq. (3) evaluated for a typical particle (light-blue bins), with a bin size Δ =  /2, compared with the analytic integrand (black line).The integrand is well behaved for most of the radial range of the orbit, and therefore well approximated by the binned density.However, it has two integrable divergences at pericentre and apocentre (here located at ∼ 2.2 kpc and ∼ 8.7 kpc, respectively).Even if the particle never reaches the centre of one of these extremal bins, it may still spend significant time within the bin.Capturing this effect correctly in the binned probability requires the special treatment explained in the text.The dark-blue shaded areas represent the analytical corrections added at the pericentre and apocentre for this orbit.

Figure 2 .
Figure2.Density profiles multiplied by  2 (left) and images of the dark matter density projected down the  axis (right) for our two lowest-mass dwarf haloes ( ∼ 2 × 10 9 M ⊙ ).The dynamical density profiles obtained from the low resolution snapshots (blue lines) agree very well with both the low and high resolution binned profiles (pink and black points) for most of the radial extent of all the haloes.The largest variations between the dynamical and binned estimates are observed in the outer regions, beyond the virial radius, where large substructures in the outskirts cause spikes in the mass distribution.Any such substructures with mass greater than 3% of the mass of the main halo are shown by brown circles in the halo images, and by corresponding brown arrows in the dynamical profile plots.The panels below the density profiles show the variations in the average radial velocity of the particles contained within concentric shells as a fraction of the virial velocity, which can be used to quantify how close the low resolution halo is to equilibrium.The pink arrows indicate the radius corresponding to 3 times the value of the softening scale of the low resolution simulations (i.e. conv for the low resolution binned profiles).

Figure 3 .
Figure 3. Same as Figure 2 but for the three intermediate-mass dwarf haloes ( ∼ 5 × 10 9 M ⊙ ).Similarly to the other cases, the dynamical density profiles from the low resolution snapshots agree well with both binned profiles.Halo600 is an outlier since it recently had a merger close to the halo's centre which disrupted the equilibrium in the inner regions; as a result the plot of v /v vir shows significant deviations from zero at small radii.Halo624 has a large substructure within its virial radius which will reach the centre of the main halo and merge with it in the next ∼ 500 Myrs.(The structure is found slightly closer to the centre in the high resolution simulation.)The significant disruption caused by this substructure to the halo's equilibrium is also evident in the average radial velocity panel, but our dynamical method nonetheless recovers a sensible 'smoothed' density profile.

Figure 4 .
Figure 4. Same as Figures2 and 3but for the two most massive ( ∼ 10 12 M ⊙ ) out of all seven haloes.Similarly to the other haloes, the dynamical density profiles from the low resolution snapshots agree well with both the low and high resolution binned profiles.For efficiency, the dynamical profiles for these haloes were generated using only a randomly selected fraction (a third) of all the particles within the halo and therefore even smaller errors on the dynamical density profile are achievable in principle.In these examples, all substructures are small (less than 1% of the halo mass) and do not have a visible effect on the density profiles.

Figure 5 .
Figure5.Binned density profile multiplied by  2 obtained by stacking 6 consecutive snapshots from the low resolution simulation of Halo1459 (black points).The dynamical density profile (blue line) and the binned profile obtained from a single snapshot at  = 0 (pink points) are also shown for comparison.Although Poisson noise is mitigated in the stacked profile, the method cannot correct the systematic softening and relaxation errors, and therefore underestimates central densities, unlike the dynamical profile.
Figure7.Dynamical density profile multiplied by  2 (blue line) obtained from the high resolution simulation of Halo1459, compared to the binned density profiles of the high (black points) and ultra-high (green points) resolution snapshots.The binned profile obtained from the low resolution snapshot is shown for reference (pink points).The black arrow indicates the approximate convergence radius of the high resolution binned profile (3 ).The dynamical density profile from the high resolution simulation predicts the ultra-high resolution simulation well, underscoring how the method can be applied at any resolution to extract additional information.

Figure 8 .
Figure8.Same as Figure7but for Halo600.The dynamical profile from the high resolution simulation of this halo shows a steep cusp consistent with the ultra-high resolution binned profile.The high resolution simulation, unlike the low resolution version, did not recently undergo a merger close to the halo's centre.This provides further evidence that the disagreement between the dynamical and binned profiles seen at small radii in the low resolution case is due to disequilibrium caused by the merger event.

Figure 9 .
Figure9.Zoom into the outer regions of the dynamical profile of Halo605 (middle of Figure3) when the reflecting boundary is placed at 100 kpc (blue line) and then moved to 200 kpc (grey line) compared to the low resolution binned profile (pink points).The dynamical profile agrees well with the binned one when the boundary is placed anywhere up to 100 kpc, around twice the virial radius, but differs once contributions from particles out to 200 kpc are included in the calculations.These discrepancies propagate inwards to smaller radii, even below the virial radius (55 kpc, indicated by the vertical dashed line).This behaviour reflects our algorithm's extrapolation of how particles and substructures in the outskirts, while currently unbound, will ultimately fall into the halo at later times in the simulation, altering the density distribution.

Table 1 .
Properties (softening length, particle mass, number of particles, virial radius, virial mass, and brief comments on the density structure) of the seven haloes investigated in this work.The haloes can be grouped into 3 main categories based on their virial mass, from dwarf to Milky Way mass.The number of particles refers to the particles enclosed by each halo's virial radius at  = 0.