Cosmic structures from a mathematical perspective 1: dark matter halo mass density profiles

The shapes of individual self-gravitating structures of an ensemble of identical, collisionless particles have remained elusive for decades. In particular, a reason why mass density profiles like the Navarro–Frenk–White or the Einasto profile are good fits to simulation- and observation-based dark matter halos has not been found. Given the class of three dimensional, spherically symmetric power-law probability density distributions to locate individual particles in the ensemble mentioned above, we derive the constraining equation for the power-law index for the most and least likely joint ensemble configurations. We find that any dark matter halo can be partitioned into three regions: a core, an intermediate part, and an outskirts part up to boundary radius rmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_\mathrm {max}$$\end{document}. The power-law index of the core is determined by the mean radius of the particle distribution within the core. The intermediate region becomes isothermal in the limit of infinitely many particles. The slope of the mass density profile far from the centre is determined by the choice of rmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_\mathrm {max}$$\end{document} with respect to the outmost halo particle, such that two typical limiting cases arise that explain the r-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{-3}$$\end{document}-slope for galaxy-cluster outskirts and the r-4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{-4}$$\end{document}-slope for galactic outskirts. Hence, we succeed in deriving the mass density profiles of common fitting functions from a general viewpoint. These results also allow to find a simple explanation for the cusp-core-problem and to separate the halo description from its dynamics.


Introduction
During the past decades, large-scale N-body simulations have successfully reconstructed cosmic structure formation with increasing resolution and complexity, as observations corroborate, see e.g. [5,9,11], or [18] for summarising overviews. Complementary efforts to search for a closed-form description have arrived at a hydrodynamical theory that explains cosmic structure evolution up to the non-linear regime, see [30] for a recent overview. A new approach to describe cosmic structure formation based on a kinetic field theory can derive an analytic, parameter-free equation for the non-linear cosmic power spectrum, see [3] for an introduction.
All of these approaches describe cosmic structures as evolving from initial overdensity seeds in the power spectrum that are propagated forward in cosmic time. To yield the non-linear agglomerations we observe today, these approaches thus require an initial configuration of overdensity seeds in phase-space together with a dynamical model of its evolution over cosmic time. As a result, they obtain the statistical properties of mass density perturbations of the observable universe as a whole based on the evolution of the phase-space-volume of the initial overdensity seeds. Subsequent attempts to derive the shape of individual mass density profiles of simulated and observed galaxies or galaxy clusters from these approaches have not been successful yet. A potential reason that the mass density profiles have not been explained by one of these approaches so far is the missing physical link between fluctuations in the power spectrum of the mass overdensity and the definition of an individual cosmic structure, like a (dark matter) halo.
Several heuristic fitting functions have been developed that capture the shape of individual dark matter halos, among others the Navarro-Frenk-White profile [26] and the Einasto profile [10], and the shape of individual luminous matter distributions, like the de-Vaucouleurs profile [8] or the Jaffe profile [15]. Without a deeper understanding how these profiles can be derived from more fundamental principles, it still remains an open question why these formulae are good fits to the simulated and observed mass density distributions.
To further investigate this question, we develop a statistical approach to derive the shape of dark matter halo density profiles from the framework of probability theory. The approach is based on an ensemble of independent, identical dark matter particles that are spatially spread according to generic, spherically symmetric power-law probability density distributions which are typical for scale-free interactions like gravitation. Although we only consider ensembles of dark matter particles, the approach is generic and can be applied to any ensemble of particles that shares the same prerequisites.
Contrary to the approaches mentioned above, we split cosmic structure evolution into the description of structures at one instant in cosmic time and a second part describing their dynamical evolution. In this work, we focus on the former part and leave out the latter. This is feasible because we do not track the change of initial structures across cosmic time to model a current one. Instead, at one instant in cosmic time, we define a structure as an agglomeration of particles on a predefined length scale of interest. This approach is supported by the simulations of [24], showing that the mass density profiles of dark matter halos at one time seem to be independent of their merger history. [14] arrive at similar results, when setting up their derivation of dark matter halo mass density profiles. They find that approaches to explain dark matter halo mass density profiles based on statistical mechanics, including their own, deal with final equilibrium states of the halo mass densities. Hence, the resulting dark matter halo mass density profiles, like their "DARKexp" final equilibrium state, are independent of any dynamical relaxation process.
Defining a (cosmic) structure as a particle ensemble at a specific scale enables us to set up a general characterisation of a dark matter halo that is solely subject to gravitational forces acting on this scale. Our effective approach is thus independent of the processes happening on smaller scales and of the dynamics on any scale leading to the state under consideration, similar to [14], but without the need to introduce an energy-state-space or a phase-space. We find that the parametric fitting functions mentioned above can be derived from this general characterisation as combinations of limiting cases for the halo core, the isothermal middle part, and the boundary regions, taking into account the finiteness of the observation or simulation to be described, meaning their finite resolution, finite volume, and finite amount of dark matter particles.
The work is organised as follows: In Sect. 2, we briefly summarise the parametric mass density profiles for the most commonly used dark matter halo profiles of galaxies and galaxy clusters. Given a small set of assumptions as prerequisites, detailed in Sect. 3, we set up a maximum-likelihood approach to derive the power-law index for the ensemble of dark matter particles. In Sect. 4, limiting approximations of the general case set up in Sect. 3 are considered. We derive the power-law index for the core, the isothermal, intermediate region, and describe the boundary regions for galaxy-scale and galaxy-cluster-scale dark matter halos. In Sect. 5, we set up a mass density profile from the particle ensemble, such that the mathematical results derived in the previous sections can be assigned a physically meaningful interpretation. In this way, we can explain the shape of the most commonly employed halo mass density profiles of Sect. 2. We conclude in Sect. 6 with a summary of the approach, its current explanatory power, and an outlook to the next parts of this series.

Halo mass density profiles from a physical viewpoint
Our characterisation of cosmic structures decisively depends on the definition thereof. The most important property of a structure is that its constituents are gravitationally bound to each other. The structure should be of finite extent and separable from similar structures in its vicinity, unless structures merge. Furthermore, a structure requires a characteristic (length) scale. Mathematically, a characteristic (length) scale serves as a scaling (radius) to obtain dimensionless variables. Physically, the scaling (radius) is important to mark characteristic scales of potentially occurring phase transitions.
Tracing all individual constituents within cosmic structures and taking into account their mass to set up a mass density distribution of these structures is computationally infeasible for luminous matter and principally impossible for dark matter. Therefore, we describe the mass density distribution of gravitationally bound, separable structures of finite extent by a continuous mass density ρ(r), averaging the distribution of particles over scales that are larger than their mutual mean free path length and that are smaller than the extent of the structure, see [4] for details.
We briefly discuss the most common mass density profiles for galaxy-and galaxycluster-sized dark matter halos. An extended list can be found e.g. in [7] and [31]. All models are designed to characterise the mass density profile of a dark matter halo with a small number of parameters that have a physical interpretation. We restrict our analysis to spherically symmetric mass density profiles, i.e. we neglect the perturbations by interactions with the environment. Small perturbations are usually modelled by introducing elliptical isocontours, larger deviations from spherical symmetry often require to use more complex mass density profiles. The assumption of spherical symmetry is a good approximation to the inner parts of galaxy clusters and heavy, elliptical galaxies. Making this simplification, we can derive the most-likely parameter values analytically as a proof of concept. These derivations can be extended to cases with less symmetry in a future work.

Navarro-Frenk-White profile
The three-dimensional Navarro-Frenk-White (NFW) mass density profile, as introduced in [26] and [27], reads in which r denotes the radius, r s a scale radius. ρ s is the scale density. It is four times the mass density at r s , i.e. ρ s = 4ρ NFW (r s ). As stated in [11] and references therein, the physical origins of this profile have not yet been fully understood. Usually, r s is defined as the radius for which the logarithmic slope of the mass density profile is -2. For a general limiting radius r max , the NFW-halo mass is given as Note that M(0) = 0, despite the singular core of the mass density profile at r = 0. Usually, Eq. 1 is integrated to r max = r 200 . This radius is determined by assuming that the mass of the dark matter halo up to r 200 equals the mass of a sphere with the same radius but a constant density which equals 200 times the critical density of the homogeneous cosmic background density, ρ 200 . This radius is chosen as an approximation to the radius at which the halo is in virial equilibrium, called r vir .
Instead of characterising the NFW profile by ρ s and r s , it is common to use its mass given by Eq. 2 and the concentration parameter c With the definition of c and ρ 200 , Eq. 2 for r max = r 200 ≈ r vir becomes Fitting the NFW profile to equilibrium mass densities in simulations, [27] found that it can describe these mass density profiles well over two decades in radius independent of the halo mass, the initial power spectrum, and the values of the cosmological parameters. Depending on the resolution of the data set, other mass density profiles may yield better fits. Yet, the NFW profile is still one of the most commonly employed models of sufficient fitting quality for a lot of applications.

Moore profile
Another dark matter halo profile commonly discussed in the literature, is the Moore profile (MP, for short). It is defined in [24] by As for Eq. 1, r s denotes a scale radius, usually chosen as mentioned above. While its behaviour for large radii is the same as for the NFW profile, it was defined to be a modification of the NFW profile with a steeper slope in the central part of the halo.

Einasto profile
A third common three-dimensional dark matter halo profile form is given by the Einasto profile (EP), as defined in [10]. The profile has the form in which the scaling radius r s is usually chosen as mentioned above. The additional third parameter of the profile, β, is called the shape parameter. As β allows for a change in shape over the radius, which is not possible in the NFW or Moore profile, the Einasto profile has been found to better fit highly resolved N-body simulations of dark matter halos, [12,13,21,22,28].

Pseudo-isothermal elliptical mass profile
The fourth commonly used dark matter halo profile is the pseudo-isothermal elliptical mass distribution (PIEMD). It is introduced in [16] as with a scaling density ρ s and two scaling radii, r s and r c , denoting the scale radius of the halo core and of a cut-off radius to truncate the halo mass density, respectively. In [16], the radius r is assumed to constrain an elliptical isocontour, however, we do not use the ellipticity here and assume spherically symmetric isocontours for the halo mass density. The exponent of the second term in the denominator is usually set to n = 1. The behaviour in the centre of the halo depends on the size of r s relative to r c . For r r s and r r c , and n = 1, the density falls of as r −4 . This profile has been successfully used to constrain projected halo mass profiles from observations of multiple images of gravitationally lensed background sources, see e.g. [6,17,25]. While [16,19], and [25] use Eq. 7 to describe the dark matter halos of elliptical galaxies, [6] and [17], among many others, also use it for galaxy-cluster scale mass density profiles.

Jaffe profile
The fifth commonly used dark matter halo profile, most often employed for galaxyscale dark matter halos, is the Jaffe profile (JP), as introduced in [15] As for the other profiles, ρ s and r s denote the scaling density and scale radius, respectively. For r r s , the density in the halo centre scales like r −2 , i.e. even steeper than the MP, and it falls off like r −4 for r r s , similar to the PIEMD. Yet, contrary to the latter, Eq. 8 has a cuspy core, given by the first bracket in the denominator, and only one scale radius r s . Projecting the JP on the two-dimensional sky, the de-Vaucouleurs profile, [8], is obtained. Indeed, Eq. 8 is motivated as a corresponding three-dimensional mass density profile for the the projected, observed light distributions of elliptical galaxies and bulges of spiral galaxies that are heuristically fitted very well by de-Vaucouleurs profiles.

Halo mass density profiles from a statistical viewpoint
Having introduced the physical mass density profiles, we now leave the mass-density approach until Sect. 5 and pursue the statistical approach to describe a dark matter halo as an ensemble of particles. The prerequisites for the calculations performed in Sects. 3 and 4 are stated in Sect. 3.1. Subsequently, we introduce the probability distribution to describe the location of an individual particle in a cosmic structure in Sect. 3.2. The joint probability density distribution for an ensemble of identically distributed particles is set up in Sect. 3.3. We obtain the general equation to determine the power-law index for the one-particle probability distributions from the extrema of the joint distribution of the ensemble.

Prerequisites
As a general starting point, we consider n p independent, identical, and identically distributed particles of mass m, located at positions r j , j = 1, ..., n p in a gravitationally bound structure. We restrict ourselves to the case of a spherically symmetric structure (see Sect. 2). Furthermore, we assume that the particles are only interacting with each other by weak-field gravity. By the latter, we mean any gravitational law that can be effectively described as a scale-free interaction on a constant number of individual particles, which is also the basis of the hydrodynamical approach, the simulations, and the kinetic field theory of [3] because observations of the large-scale matter power spectrum support this assumption, see e.g. [29]. A good example fulfilling these requirements would be an ensemble of dark matter particles that interact with each other by Newtonian gravity to form a dark matter halo.
Here, the meaning of "particle" depends on the scale that we want to describe the structure on. In simulations, for instance, "particles" can be dark matter agglomerations of several thousand solar masses. In observations, the smallest "particles" that can be described are given by the resolution at which the dark matter halo is constrained given all measurement imprecisions and uncertainties of the modelling.
The assumption of independent particles implies that they are uncorrelated. This is similar to approaches mentioned in [14]. It is a reasonable starting point because any correlation between the particles requires additional interactions and evidence for additional dark matter particle interactions apart from gravitation is still to be observed. Including correlations requires additional knowledge of the dynamical status of the structure such that correlations are traced back to a previous state of the structure -which we do not take into account in this work. As we will show in the subsequent parts of the series, this assumption seems to be a good approximation for current simulations and observations.
Whether the particles in the ensemble are distinguishable or not, is of no concern. We will see in the course of the calculations that our approach treats distinguishable and indistinguishable particles alike. The equations to constrain the power-law index are invariant with respect to this property.

The localisation probability of a particle within a structure
As already stated in Sect. 3.1, following the principle of indifference, we assume that the identical particles are identically distributed. Given that the interaction which is supposed to generate the halo shape is scale-free, let be the probability of a single particle to be located at radial position r in a three dimensional, spherically symmetric mass density distribution. To obtain a dimensionless radius, we scale r by a scale radius r σ (0 < r σ < ∞). From the mathematical point of view, the scaling to a dimensionless radius is required for the log-likelihood of the joint probability distribution of the entire ensemble. α is the power-law index. We require α > −1 to ensure that the probability density distribution monotonically decreases with increasing distance from the halo centre.
We are free to choose r σ such that r > r σ for all particles. This allows us to drop the 1 in Eq. 9 because, r > r σ avoids the singularity of a standard power-law at r = 0. p(r ) thus remains finite and well-defined. We determine the normalisation constant N from the requirement that the particle belonging to the structure is located in the volume that is delimited by the minimum radius, r min , and maximum radius of the structure, r max (r max > r min > r σ ), i.e.
Solving Eq. 10 as detailed in Appendix 1 we arrive at in which, using the abbreviation x max ≡ r max /r σ (and analogously for x min ), Hence, the normalisation constant is a function of α, r σ , r min , and r max . Combining Eqs. 9 and 11, we obtain the normalised probability distribution for a particle being localised at radius r within the structure.

The power-law index of the most likely ensemble configuration
Next, we determine the log-likelihood for the n p independent and identically distributed particles as The maximum-likelihood estimator for α is then determined by solving dL/dα = 0. From Eq. 14, we obtain as derived in detail in Appendix 1. Using Eq. 12, we obtain for the first term in Eq. 15 Assuming that x max x min for α = 2, we arrive at As mentioned in Sect. 3.1, we find that Eq. 15 is the same for distinguishable and indistinguishable particles because indistiguishable particles only require a global factor of 1/n p ! multiplied to the product of probabilities in Eq. 14. As this factor is a constant, it does not occur in the derivative of L(α) in Eq. 15.
In fact, Eq. 15 yields an extremum of the likelihood and not necessarily its maximum. It thus remains to be determined whether the solution found is a likelihood minimum, saddle point, or maximum.
We note that α enters Eq. 15 over the normalisation N , i.e. the first term, which is related to the predefined extent of the structure under consideration and does not depend on the distribution of the n p particles. This distribution is contained in the second term. Although we described the probability of one particle to be at a certain location by a power law, the actual distribution of an ensemble of simulated or observed particles may deviate from this assumption. In Sect. 4, we therefore consider different actual distributions for the second term and investigate different approximations to Eq. 15 to determine α for the most common physically motivated limiting cases. For each of them, we also calculate, which extremum of the likelihood they represent.
Analogously, we can derive Eq. 14 with respect to r σ , r min and r max to determine extremum configurations for the scaling radius and the boundaries of the structure. The derivations can also be found in Appendix 1. For r σ , we obtain, This result implies there is no preferred scale radius r σ , which is expected by construction, assuming that r σ was introduced as an auxiliary parameter to obtain dimensionless variables and assuming r j > r σ for all particles j = 1, ..., n p . The calculation for r max yields (20) we conclude that, for arbitrary halo boundaries, the derivative of the likelihood with respect to r max is always zero for α = 2. Since r max > r σ > 0, the last term in Eq. 20 only approaches zero for r max → ∞. This result implies that, given a finite r max , like r 200 , the only extremum is at α = 2. The calculation for r min is analogous and discussed in Appendix 1 because r min can simply be identified with the radial position of the innermost particle of the ensemble. We can summarise the results so far that, as expected from a scale-free ansatz for p(r ), no preferred values for the scaling radius r σ or the halo boundaries r min and r max can be derived. Therefore, we continue our analysis with Eq. 15 to find extrema for the power-law index α. Equation 15 depends on the extension of the structure, given by the maximum and minimum radius, r max and r min , the scale radius r σ , the number of particles belonging to the structure, n p , and their radial positions r j , j = 1, ..., n p . Thus, different approximations to Eq. 15 are obtained depending on the relationship between the pre-defined halo extent, r max and r min , the scaling radius r σ , and the radial particle positions r j with respect to each other and the number of particles belonging to the structure. In a simulation, the latter amounts to changing the resolution of the structures.

Power-law indices for different approximations
Having chosen r σ < r min ≤ r j , we can then distinguish four different limiting cases. These four cases, subsequently discussed in individual sections, are: 1. The first case, we call the "core case", considers the central halo part from r min to the boundary of the core, called r core . In this case, we assume r min < r j < r core r max . We furthemore introduce n c as the number of particles in the core and assume that it is small compared to n p , the total amount of particles in the entire halo, so that the individual particles move freely in the range between r min and r core . 2. The second case is the homogeneous fluid case, in which the particles are homogeneously distributed in the volume between r core and r max . This case requires large particle numbers up to r max , so that we consider n p → ∞ here. 3. In the third case, we call the "thinning outskirts case", we let r max go to large extensions, r max → ∞, while the particles assigned to the structure are assumed to remain at a much smaller radius, r j ∞, j = 1, ..., n p . This case ignores particles in the simulated or observed particle ensemble that are gravitationally bound to a structure but are located far from its centre. 4. The fourth case is called the "overflowing outskirts case", in which we assume r max is chosen too small, such that the particle distribution assigned to that structure spills over the predefined maximum extension.

Core case
At first, we consider the central part of the halo from r min up to a radius r core that we call the boundary of the halo core. Inspired by simulations, we take r min to be the inner-most radius that is available, r core to be the radius at which the simulation has converged and r max to be r 200 . For the simulations of [28], r min = 1.5 · 10 −4 r 200 and r core = 0.6 · 10 −3 r 200 , such that r min /r core = 2.7 · 10 −2 . We thus assume that r min r core , so that we can use Eq. 15 and insert Eq. 17 with r core now taking the role of r max as the outer boundary of the central halo part to obtain 1 which, using that x j /x core = r j /r core , can be simplified as The extremum values of α are thus independent of r σ and only subject to the number of particles in the core, n c , and their distribution. The accuracy of the approximation of the logarithm on the right-hand side increases, the closer r j is to r core . Since we assume that 0 < r j < r core , convergence of the full Taylor expansion is fulfilled. For a homogeneous distribution of particle radii between r min and r core , we have such that Assuming that it is not the radius which is uniformly distributed but the threedimensional particle position in the volume of the halo, we obtain such that α = −2, which is already smaller than the smallest possible α = −1. Hence, the specific particle distribution directly determines α. For the uniform distribution of particle radii according to Eq. 24, we find p(r ) ∝ r −1 . Shifting the mean radius of the particle distribution towards r core leads to shallower profiles, until we reach p(r ) ∝ r 0 . In simulations, p(r ) is constant per bin, such that p(r ) ∝ r 0 occurs naturally as soon as r core = r min , when the core only covers the first bin filled with particles. Similarly, probability densities steeper than p(r ) ∝ r −1 , i.e. α > 0 are obtained when the mean radius of the particle distribution is shifted away from r core . Any value of α < 2 given by Eq. 24 leads to a likelihood maximum.

Homogeneous fluid case
Next, we describe intermediate part of the halo with the particle distribution from r core up to the halo boundary r max . Since r min and r core are both much smaller than r max , for the sake of simplicity, we consider the entire halo from r min to r max .
As in Sect. 4.1, we insert Eq. 17 into Eq. 15 and solve for α to obtain With 0 < r min ≤ r j ≤ r max < ∞, all terms in the denominator are negative yielding α < 2.
Assuming that n p 1, the particle distribution becomes a homogeneous fluid and the particles are located at positions This implies that the power-law index of the extremum particle distribution becomes independent of r max and only depends on n p . Assuming that we narrow the radius range of a power law and cut off the strongly decreasing tail, we can understand the homogeneous fluid distribution as a detail of a power law distribution, so that the actual particle distribution can still be a power law. We approximate the sum in the denominator by Stirling's formula as detailed in Appendix 1 to arrive at lim n p →∞ α(n p ) = lim which belongs to a likelihood maximum. Inserting this power-law index into Eq. 13, we arrive at a power-law distribution for the individual particles that scales with r −2 .

Thinning outskirts case
For the "thinning outskirts case", we reuse Eq. 27. Instead of assuming a homogeneous distribution for the sum term in the denominator, we assume that in which we have to reintroduce x max = r max /r σ again to make the argument of the logarithm dimensionless. This approximation, ln(x j ) ≈ 0 for all j = 1, ..., n p , implies that we do not take into account particles in the tail of the power-law distribution. These particles are gravitationally bound to the structure but may not be assigned to this structure because they are far away from its centre. In simulations or observations, they may also be included into neighbouring structures to which they are closer.
Inserting Eq. 30 into Eq. 27, we obtain Letting r max approach the range of the gravitational interaction length, we obtain lim which results in a likelihood minimum and a one-particle power-law distribution that scales with r −3 .

Overflowing outskirts case
Lastly, we treat the "overflowing outskirts case" by reusing Eq. 27 again and assuming that the volume of the structure is set too small, such that r max ≤ r j for all j = 1, ..., n p , i.e. the particle distribution spills completely over the defined extension of the structure. Furthermore, we assume that the r j and r max are of the same order of magnitude, such that r = r max + δr , with δr ≤ r max (33) holds for the arithmetic mean value r of the radial positions r j and the arithmetic mean δr of their deviations δr j from r max . Given these assumptions, further motivated below, we can now approximate the sum in the denominator of Eq. 27 by using Eq. 33 In the last step, we defined 0 < k ≤ 1. Thus, we obtain Letting k approach the maximum value, we obtain which results in a likelihood maximum and a one-particle power-law distribution that scales with r −4 . If we do not impose the assumptions stated in Eq. 33, but instead assume that r j r max for all j = 1, ..., n p , the sum in Eq. 34 can still be set equal to n p k, now for k 1. Taking the limit of k → ∞ implies that the arithmetic mean of the actual particle distribution goes to infinity for a finite r max . We then obtain the same limit as in Eq. 32 again. This case may occur, if it is very hard to separate the structure of interest from its environment and its volume is assumed too small compared to the actual extend of the structure. Comparing this case with the thinning outskirts case, we thus arrive at consistent results, as we do not account for particles far away from the centre which should be assigned to the structure in both cases. If we can set r max at least at the right order of magnitude with the second constraint in Eq. 33, we can read off Eq. 35 that steeper one-particle power-law distributions than r −4 are obtained, if 0 < δr < r max . Hence, analogously to the core case (Sect. 4.1), the specific distribution of particles with respect to r max fixes α in the boundary regions of a structure.
Since α = 3 > 2, the approximation of Eq. 17 using ln(r max ) seems inappropriate to use here. Yet, requiring that the particle distribution is located at radii all larger than r max , the condition that r max > r min used in this approximation is also inverted. We can thus continue to employ this approximation, as it remains valid for this case.

From particle ensembles to mass density profiles
In statistical mechanics, the transfer from the positions and momenta of an ensemble of n p particles to continuous spatial number densities and densities in momentum space is induced by defining a probability distribution function in phase-space. This function assigns each element of the phase-space volume a probability that the configuration of the positions and momenta of the particle ensemble currently under consideration lies in this phase-space volume-element, see e.g. [2] for an introduction. The size of the infinitesimal phase-space volume-element thus sets a scale below which individual particle properties are averaged over. For a collisionless system of particles, the individual one-particle probability density functions in phase-space are uncorrelated, which implies that the probability distribution function of the ensemble is separable and amounts to a product of one-particle probability density functions, see [4]. Thus the prerequisites of collisionless systems are consistent with our prerequisites in Sect. 3.1, especially the assumption that the particles are independent, and hence allow us to link both approaches.
In the framework of phase-space probabilities, the number density of particles at a position r , n(r ), is given as the probability density function in phase-space for one particle integrated over all its potential momenta and normalised to obtain all n p particles when integrating over the entire volume. Consequently, the number density of particles can be determined from the one-particle probability density function of Eq. 9 as in which n p enters due to the normalisation condition mentioned above. Using Eq. 37, it is straighforward to determine the mass density of a structure consisting of identical particles with mass m that are distributed according to the Fig. 1 Schematic composition of the most likely dark matter halo shape based on the probabilistic derivations and approximations detailed in Sect. 3 (here shown with "thinning outskirts") number density n(r ) as In the last step, we introduced a scaling density ρ σ in analogy to the scaling density ρ s for the mass densitiy profiles of Sect. 2. Integrating ρ(r ), we obtain the mass of the structure up to radius r M(r ) = mn p 4π r r min drr 2 p(r ).
As a result from this transfer, we find that the mass density profile of a structure ρ(r ) consisting of an ensemble of collisionless particles is directly proportional to the one-particle probability distribution p(r ) as introduced in this section. Consequently, we can assemble a dark matter halo in its most likely shape as depicted in Fig. 1 out of an ensemble of dark matter particles fulfilling the prerequisites of Sect. 3.1. The core: The innermost part of the halo is given by radial positions that are smaller than some chosen core radius, called r core . There, dark matter forms an ideal gas of collisionless freely moving particles. The slope of the corresponding mass density profile is determined by the specific distribution of the particles. The arithmetic mean radius around which the particle distribution is concentrated can be employed as its characteristic length. If the radii of the particles are uniformly distributed, a ρ(r ) ∝ 1/r -behaviour is observed (see Sect. 4.1). Shallower mass density profiles are expected for a mean radius closer to r core , steeper profiles for a mean radius closer to the halo centre. The homogeneous fluid region: With increasing radial distance from the halo centre, the dark matter particles become more numerous and dense to form a homogeneous fluid for the limiting case of infinitely many particles (see Sect. 4.2). The same limiting ρ(r ) ∝ 1/r 2 -behaviour as for this region is also derived for the mass density profile of a singular isothermal sphere using the equation of hydrostatic equilibrium for an ideal gas density, see e.g. [4] for details. The former arises as the most likely configuration that any self-gravitating system of independent particles strives towards. The latter is the maximum-entropy configuration of the gas in phase space, [2]. Thus, the ideal isothermal gas with its smooth mass density profile has the same properties as our fluid approximation of the collisionless dark matter particle ensemble. In this sense, our mathematical derivation of α from the most likely spatial configuration of an ensemble of particles in Sect. 4.2 based on probability theory is consistent with the physical one based on statistical mechanics and thermodynamics.
The boundary region ("outskirts"): How the structure is separated from others or from a background, strongly depends on the definition of the structure and its specific environment. As no permanent thermodynamical equilibrium states exist due to the infinite range of the gravitational interaction, the definition of a maximum radius to separate structures is arbitrary compared to standard definitions based on equilibrium states. Consequently, depending on the location at which we set r max with respect to the ensemble mean of the particle distribution that forms the halo of interest, we obtain different slopes for the decreasing power-laws of ρ(r ) (see Sects. 4.3 and 4.4). As Eq. 19 reveals, ρ(r ) ∝ r −3 is preferred to obtain an extremum of the likelihood for any chosen r max < ∞. In any of those cases, the derivations show that the respective power-law slope is set by discarding weakly-bound particles far from the halo centre.

Explanation for the most commonly used mass density profiles
As a summarising result of this section, we disentangled a dark matter halo into a core, an intermediate, and an outskirts part, based on a mass density profile that is proportional to the power law in Eq. 9. Comparing the extremal power-law indices α for these halo parts, we find that they vary with the distance from the halo centre and decisively depend on the predefined halo extensions and given morphological symmetry.
In addition, α depends on the underlying distribution of dark matter particles. Hence, we can now comprehend why the most commonly used dark matter halo mass density profiles of Sect. 2 are good fits to most of the simulated dark matter halos and those inferred from observations: All heuristic models are able to produce the varying α over the extension of the halo. Starting at the core, the profile slope steepens with increasing radius, such that α monotonically increases with increasing radius. For the NFW profile, the MP, and the JP, the mass density profile consists of a product of two power laws which can be directly related to the approximating cases of Eq. 15 as summarised in Table 1. These profiles thus seem to be constructed based on the physically motivated prerequisites stated in Sect. 3.1, in particular on the assumption that the scale-free gravitational interaction can be described by a powerlaw probability density. The varying αs between these models account for the fact that different approximations to Eq. 15 lead to slightly varying extremal power-law indices, as detailed in Sect. 4. Simulation-and observation-based examples for these approximations will be shown in the next parts of the series.
Potential reasons for the debates about the profile slope in the core and in the outskirts part can now be explained by the result that the extremal power-law index for these parts are highly dependent on the specific particle distribution. In addition, the extremum of the thinning outskirts case is a likelihood minimum. Contrary to that, all density profiles agree that halos show an isothermal part between the core and the outskirts, which belongs to the most likely spatial configuration of a many-particle ensemble and is therefore ubiquitous in simulated and observed cosmic structures. The overflowing outskirts is also a likelihood maximum. Its decrease of the mass density profile is often found for galaxy-scale masses. On that scale, the extension of the structure is characterised by the half-light radius, which seems to be a good estimate for the order of magnitude on which the mass is distributed as well and which motivated the approximations stated in Eq. 33.
For the EP and the PIEMD, the integration into our approach is different, as they are not based on the power law of Eq. 9. In order to find their extremal α, the calculations of Sects. 3 and 4 need to be repeated using their specific probability density distributions. Instead of doing this, we show that the two density profiles can approximate the power law of Eq. 9, such that the EP and the PIEMD have a similar behaviour as the other profiles. Consequently, this corroborates their good fits to simulations and observations of dark matter dominated structures. Further examples for successful fits to simulations and observational data will follow alongside the ones mentioned above in the next parts of this series.
EPs for β 1 approximate an isothermal sphere because the limit of β → 0 for Eq. 6 yields For β = 0 and r r s , the exponential function can be Taylor-expanded to yield such that the core of an EP is similar to a power law. For the PIEMD, the profile is a product of two different power laws constructed to yield the same behaviour as the power-law profiles mentioned above for large r . For r → 0, it avoids any singularity and given r r s ≤ r c , we obtain to leading order which is similar to the next-to-leading-order Taylor expansion of Eq. 9 without the linear term in r .

Conclusion
The main aim behind this new approach is to complement current models of structure growth that intertwine the description of a cosmic structure with its dynamics and thus require to track initial seed structures through cosmic time in order to explain the non-linear mass agglomerations we observe today. In the first part of this series on structure evolution over cosmic time, we focussed on separating the description of a cosmic structure from its dynamics and reducing the former to the necessary amount of prerequisites and characterising equations. Being able to find a description of any structure at an arbitrary point in cosmic time that can be used as a starting point for any theory of structure growth is a first step to making our standard approaches more efficient and to envisioning tests of alternatives that have been computationally intractable so far. Our goal in this first part was to derive a general description of dark matter halo shapes from more fundamental principles than using heuristic models like the ones introduced in Sect. 2. Our newly gained understanding of the morphology of dark matter halos also allows us to support the currently employed gravitational lensing models to infer the Hubble-Lemaître constant (see e.g. [23] and references therein for recent investigations of potentially arising biases in the modelling) or to infer galaxy cluster masses (see e.g. [20] and references therein for recent investigations on the quality of mass density reconstructions and potentially arising biases in the modelling).
To describe the morphology of a dark matter halo, we employed the physical prerequisites of a set of collisionless, identical, only gravitationally interacting particles that are confined to a spherically limited volume called the dark matter halo. Translating these requirements into mathematical terms yields an ensemble of independent, identical particles that are identically distributed within a finite sphere according to a power law probability density distribution. We set up the joint probability density distribution of the entire ensemble and determined the power-law indices that belong to the most and least likely spatial configurations of the ensemble, given certain physically motivated limits on the confining volume and the number of particles in the ensemble, as detailed in Sect. 4.
In agreement with previous works, e.g. [26], we find that the scaling radius r σ introduced to obtain a dimensionless radial variable drops out of the equation that constrains the power-law index. Due to the scale-free approach, there is no preferred r σ (see Eq. 18). Nevertheless, introducing r −2 , the radius at which the slope of the mass density profile is −2, as scaling radius, the concentration parameter c = r −2 /r 200 yields valuable information about the scale of isothermality with respect to the scale of halo relaxation. Since isothermality is reached with an infinite number of dark matter particles in our approach, while all other parts contain a finite amount of freely moving particles, we will further investigate whether a r σ exists as a phase-transition scale between a dark matter fluid and freely moving particles. A succeeding investigation can characterise the mean free path of (potentially colliding) dark matter particles and allows to interpret the halo concentration parameter as the Knudsen number for dark matter halos.
Beyond existing approaches, we arrive at a derivation of the power-law index of the mass density profiles without the need to employ any equilibrium considerations and the velocity distribution of the particle distribution. Avoiding the combination of dynamics and structure description, our approach clearly shows the connection between a joint, spatial configuration of dark matter particles on the microscopic scale and their respective mass density profile as their averaged macroscopic representation. As a result, we find that the mass density profile belonging to a set of freely moving particles confined to a small volume in the inner-most halo part shows the same behaviour as assumed in the core region of the heuristic power-law mass density profiles introduced in Sect. 2. The fact that the extremum power-law index strongly depends on the specific distribution of the particles in the core with respect to the chosen maximum core radius, could resolve the cusp-core-debate (see [5] for an introduction) without the need to introduce any effects caused by luminous matter, as further detailed in Sect. 5.1. Analogously, we succeeded in retrieving the isothermal intermediate part between core and outskirts in our framework.
When determining the power-law indices for the outskirts on galaxy-and galaxycluster scale (see also Sect. 5.1 for further details), we found that the power-law index is set by the relation between the predefined halo boundary and the maximum radius of the particle ensemble under consideration. This result enables us to understand why the galaxy-and galaxy-cluster scale halos show a different decrease in the mass density profile far from the centre. It emerges due to the two different definitions of halo boundaries. On galaxy scale, we observe the half-light radius of the luminous part of the galaxy and assume that the respective dark matter halo is of approximately the same order of magnitude in its extensions. For galaxy-clusters, which we consider the largest cosmic structures that form the nodes in the cosmic web, we assume that, in principle, all particles irrespective of their distance to the cluster centre are gravitationally bound to the halo. At the same time, we divide the total ensemble of particles in a simulation or observation on large scales into several neighbouring clusters, such that the mass density profile for an individual cluster arises as the trade-off between the long-range effect of gravitational attraction and the assignment of particles to their nearest local minimum of the gravitational potential.
On the whole, our approach shows that an individual dark matter halo can be considered an emergent structure, i.e. an assignment of individual particles to an ensemble of particles close to a local minimum in the gravitational potential according to predefined boundary and symmetry constraints. Hence, our picture of a halo is a transient state of a set of particles agglomerating around a salient point in space. Over cosmic time, different particles can be assigned to the halo, such that its shape can be maintained with a different set of particles or its morphology can change as a consequence of the dynamics, for instance allowing for merger events.
Since we aim at describing individual halos as required by the observational applications like gravitational lensing or kinematical characterisations, our halo model is currently not linked to the power spectrum of density fluctuations, which is the basis of the approach outlined in [3]. Furthermore, it is different from the current standard model based on fluid dynamics, as already mentioned in the introduction. In this picture, halos are entities of a specified set of particles that is tracked in their collective motion over cosmic time, such that the particle dynamics introduces changes in the halo morphology.
Given its success at explaining the dark matter halo shapes in a very clear and simple way, the subsequent parts of this series will deal with a thorough consistency check of the proposed approach by dark-matter-only simulations and by observations of low surface brightness galaxies, which are assumed to mainly consist of dark matter. Beyond that, we will also apply our approach to corroborate the common mass density models for gravitational lenses, see e.g. [31] for an overview, and investigate whether it can shed light onto the bulge-halo-conspiracy, [1]. Thereby, we intend to integrate our approach into the landscape of existing methods to improve our understanding of structure growth and enhance our tool set for doing so. x min As g(α, r σ , r max , r min ) becomes singular for α = 2, we check that the limit of Eq. 47 exists. It exists and yields g(2, x max , x min ) = ln(x max ) − ln(x min ).