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. Interpreting the scaling radius as the mean free path length of the particles, a dark matter halo can be constructed as a core of free moving particles, transitioning into a homogeneous dark matter fluid. The mean radius of the specific particle distribution fixes the power-law index in the core. The fluid 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 definition of the halo boundary, such that two typical limiting cases arise that explain the $r^{-3}$-slope for galaxy-cluster outskirts and the $r^{-4}$-slope for galactic outskirts. Hence, the different power-law slopes of the individual regions allow to put limits on the properties of the dark matter particles like their number density and mean free path length. 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) and (10) 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 (27) 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 (23) and (9), and the shape of individual luminous matter distributions, like (8) or (13). 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 powerlaw 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 (21), showing that the mass density profiles of dark matter halos at one time seem to be independent of their merger history. 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. 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 Section 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 Section 3, we set up a maximum-likelihood approach to derive the power-law index for the ensemble of dark matter particles. In Section 4, limiting approximations of the general case set up in Section 3 are considered. We derive the power-law index for a dilute dark matter gas, a homogeneous dark matter fluid, and describe the boundary regions for galaxy-scale and galaxy-cluster-scale dark matter halos. In Section 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 Section 2. We conclude in Section 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 (28). 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 of-ten 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.
2.1 Common mass density profiles for dark matter halos The three-dimensional Navarro-Frenk-White (NFW) mass density profile, as introduced in (23) and (24), 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 (10) 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, Equation 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 .
Instead of characterising the NFW profile by ρ s and r s , it is common to use its mass given by Fitting the NFW profile to equilibrium mass densities in simulations, (24) 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 (21) by ρ M (r) = ρ s (r/r s ) 3/2 1 + (r/r s ) 3/2 .
As for Equation 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 (9). 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, (19), (18), (11), (12), (25).

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 (14) as ρ PIEMD (r) = ρ s 1 + (r/r s ) 2 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 (14), 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. (15), (22), (6). While (14), (22), and (16) use Equation 7 to describe the dark matter halos of elliptical galaxies, (15) and (6), 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 (13) ρ J (r) = ρ s (r/r s ) 2 (1 + r/r s ) 2 .
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, Equation 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, Equation 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 Section 5 and pursue the statistical approach to describe a dark matter halo as an ensemble of particles. The prerequisites for the calculations performed in Sections 3 and 4 are stated in Section 3.1. Subsequently, we introduce the probability distribution to describe the location of an individual particle in a cosmic structure in Section 3.2. The joint probability density distribution for an ensemble of identically distributed particles is set up in Section 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 Section 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. (26). 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 in contrast to the standard approaches mentioned in Section 1. Yet, it is a reasonable starting point in our approach because any correlation between the particles requires additional interactions -evidence for additional dark matter particle interactions apart from gravitation is still to be observed -or it 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.

3.2
The localisation probability of a particle within a structure As already stated in Section 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 σ < ∞), whose size remains to be specified by the physical application later. 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 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 ), i.e.
Due to the finite core of the probability distribution in Equation 9, we can set r min = 0. This is the main reason for choosing the Lomax power-law distribution (Equation 9) instead of the standard power-law distributioñ with its singular core, which requires r min > 0 forp(r) to remain finite and thus welldefined. Equation 11 can be employed instead of Equation 9 as soon as r min r σ . Solving Equation 10 for r min = 0 as detailed in Appendix A we arrive at in which, using the abbreviation t max ≡ 1 + r max /r σ , Hence, the normalisation constant is a function of α, r σ , and r max . Combining Equations 9 and 12, we obtain the normalised probability distribution for a particle being localised at radius r within the structure.
In the same way, we deriveÑ for Equation 11 in Appendix A so that the powerlaw distribution for r min r σ reads with, using x ≡ r/r σ , 3.3 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 Equation 17, we obtain as derived in detail in Appendix B. For general t max , d f (α, t max )/dα/ f (α, t max ) is given in Appendix B. Analogously, we obtain dL/dα for Equation 15 replacing f (α, t max ) by g(α, x max , x min ) and ln(1 + r j /r σ ) by ln(r j /r σ ). As mentioned in Section 3.1, we find that Equation 18 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 Equation 17. As this factor is a constant, it does not occur in the derivative of L(α) in Equation 18.
In fact, Equation 18 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 Equation 18 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 Section 4, we therefore consider different actual distributions for the second term and investigate different approximations to Equation 18 to determine α for the most common physically motivated limiting cases. For each of them, we also calculate, which extremum of the likelihood they represent.

Power-law indices for different approximations
Equation 18 depends on the extension of the structure, given by the maximum radius, r max , 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 Equation 18 are obtained by relating r max , r σ , and the r j with respect to each other and changing the number of particles belonging to the structure. In a simulation, the latter amounts to changing the resolution of the structures.
As the mean free path length between the particles within the structure is an important length scale to characterise the particle ensemble and potentially occurring phase transitions of the structure, we interpret the scale radius r σ as the mean free path length. 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", treats r j r max r σ with the limit of r max → r σ , i.e. the mean free path length is larger than the boundaries of the structure. In this case, the individual particles move freely in the range between r min and r max . It requires the usage of Equation 9 as one-particle probability density distribution because r min < r σ is not necessarily larger than zero. 2. The second case is the homogeneous fluid case, in which the particles are homogeneously distributed in the volume between r min and r max . This case requires large particle numbers up to r max , so that we consider n p → ∞ here. It can be described with the standard power-law ansatz of Equation 15, if we assume that r min r σ , such that r min > 0 is guaranteed.
3. In the third case, we call the "thinning outskirts case", we assume a homogeneous fluid for which r max goes 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 a homogeneous fluid for which r max is chosen too small, such that the particle distribution assigned to that structure spills over the predefined maximum extension.
Our definition of the mean free path length is understood as the length scale at which the particles transform from individually moving particles into a collectively moving fluid. Otherwise, defining the mean free path length for the particles as the average distance travelled between two particle collisions could be considered inconsistent with the prerequisite that the particles are collionless, i.e. independent.

Core case
At first, we consider the approximation in which the mean free path length is much larger than the extent of the structure. It is called the "core case" because this approximation is valid in the centre of structures close to the resolution limit of the simulations or observations. As derived in detail in Appendix A, Equation 13 can be approximated by which also fixes the first term in Equation 18. Using x j ≡ r j /r σ , the second term in Equation 18 is approximated as In the last step, we set the arithmetic mean value of all radial particle positions to k x max with some 0 < k ≤ 1 depending on the distribution of the n p particles. Carefully taking account of the absolute value when deriving f (α, x max ), Equation 18 reads which is solved by For the limit that the maximum radius goes towards the mean free path of the particles and that the radii of the particles are uniformly distributed, such that k = 1/2, we arrive at lim such that the particles in the core can have a power-law distribution that scales with r −1 under the given prerequisites. Assuming that it is not the radius which is uniformly distributed but the three-dimensional particle position in the volume of the halo, we obtain k = 3/4, which leads to and consequently a scaling of r −5/3 for the power-law distribution. Hence, the specific particle distribution directly determines α. Profiles steeper than r −1 in the core up to r −2 can be attained if the particle distribution has an arithmetic mean radial value with k > 1/2, shallower ones for 1/3 < k < 1/2. (For k ≤ 1/3, our power-law ansatz to effectively describe the gravitational interaction in the particle ensemble is not valid anymore.) Consequently, if the particles are actually distributed according to a power-law, we expect shallower slopes in Equation 22 than for the homogeneous distribution, when taking the limit x max → 1.
Any value of 1/3 < k ≤ 1 leads to a likelihood minimum for α given by Equation 22.

Homogeneous fluid case
If the particles form a homogeneous fluid, the mean free path length between particles goes to zero, i.e. r σ → 0. Hence, we can assume that 1 + r/r σ ≈ r/r σ (in particular that t max ≈ r max /r σ 1). At the same time, the number of particles needs to be large, n p 1 to consider them as a fluid. Assuming t max 1, we can approximate f (α, t max ) of Equation 13 by its leading order term and obtain such As derived in Appendix B, the same result is obtained directly using the power-law distribution of Equation 15, which demonstrates that Equation 15 is a good approximation to Equation 9 for r r σ . Inserting the result of Equation 27 into Equation 18 and setting the latter to zero, we obtain Solving for α yields Given these approximations, the particles in a homogeneous fluid are located at positions This implies that the power-law index of the most likely 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.

Inserting the positions of Equation 30 into Equation 29
, we can approximate the sum in the denominator by Stirling's formula as detailed in Appendix C to arrive at lim n p →∞ α(n p ) = lim which belongs to a likelihood maximum. Inserting this power-law index into Equation 14 or Equation 15, 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 Equation 29, which means that the particles are still considered as a fluid. Instead of assuming a homogeneous distribution for the sum term in the denominator, we assume that, for all j = 1, ..., n p This approximation, ln(t 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 Equation 32 into Equation 29, we obtain .
Letting r max approach the range of the gravitational interaction length, we obtain 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 Equation 29 again and assuming that the volume of the structure is set too small, such that r σ r max ≤ r j for all j = 1, ..., n p , i.e. the particle distribution, as a homogeneous fluid, spills 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 (35) 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 Equation 29 by using Equation 35 and t j /t max ≈ r j /r max In the last step, we defined 0 < k ≤ 1 in the same way as in Section 4.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 Equation 35, but instead assume that r j r max for all j = 1, ..., n p , the sum in Equation 36 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 Equation 34 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 Equation 35, we can read off Equation 37 that steeper one-particle power-law distributions than r −4 are obtained, if 0 < δr < r max . Hence, analogously to the core case (Section 4.1), the specific distribution of particles with respect to r max fixes α in the boundary regions of a structure.

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 Section 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 Equation 9 (or Equation 11) as in which n p enters due to the normalisation condition mentioned above. Using Equation 39, it is straighforward to determine the mass density of a structure consisting of identical particles with mass m that are distributed according to the number density n(r) as ρ(r) = mn(r) = mn p N(α, r σ , r max ) 1 + r In the last step, we introduced a scaling density ρ σ in analogy to the scaling density ρ s for the mass densitiy profiles of Section 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 Figure 1 out of an ensemble of dark matter particles fulfilling the prerequisites of Section 3.1. Fig. 1 Schematic composition of the most likely dark matter halo shape based on the probabilistic derivations and approximations detailed in Section 3 (here shown with "thinning outskirts").
The core: The innermost part of the halo is given by radial positions that are smaller than the mean free path length of the particles. There, dark matter forms an ideal gas, whose slope of the 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 Section 4.1). Shallower mass density profiles are expected for a mean radius closer to the centre, steeper profiles for a mean radius closer to the mean free path length of the particles.
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 Section 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 Section 4.2 based on probability theory is consistent with the physical one based on statistical mechanics and thermodynamics.
The boundary region ("outskirts"): The fluid approximation holds on scales that are larger than the mean free path length of the particles. Thus, it is also valid when we consider the decrease of ρ(r) far away from the centre. 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 Sections 4.3 and 4.4). 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. Table 1 Integration of all models of Section 2 into our approach using the approximations as detailed in Section 4 to Equation 18. Section 5.1 contains further details. mass density core outskirts comments profile α(k) = 2 − 1/k NFW k = 1/2 thinning uniform radial distribution of particles in the core MP k > 1/2 thinning steeper radial distribution of particles in the core EP (β = 0) k = 1 no outskirts this case is an isothermal sphere PIEMD no power-law overflowing inner part is not a power-law distribution JP k = 1 overflowing core is an isothermal fluid 5.1 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 Equation 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 Section 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 Equation 18 as summarised in Table 1. These profiles thus seem to be constructed based on the physically motivated prerequisites stated in Section 3.1, in particular on the assumption that the scale-free gravitational interaction can be described by a power-law probability density. The varying αs between these models account for the fact that different approximations to Equation 18 lead to slightly varying extremal power-law indices, as detailed in Section 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 the core belongs to a likelihood minimum and the same applies to the case of thinning outskirts. 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 Equation 35.
For the EP and the PIEMD, the integration into our approach is different, as they are not based on the power law of Equation 9. In order to find their extremal α, the calculations of Sections 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 Equation 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 Equation 6 yields ρ E (r) = ρ s (r/r s ) 2 . (42) 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 Equation 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 Section 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-Lematre constant (see e.g. (20) and references therein for recent investigations of potentially arising biases in the modelling) or to infer galaxy cluster masses (see e.g. (17) 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 Section 4.
The constraining equation on the power-law index reveals that the scaling radius is more than an arbitrary length scale used to arrive at dimensionless equations, as has been assumed so far, see e.g. (23). The scaling radius marks the length scale at which an ensemble of freely moving individual dark matter particles transfers into a dark matter fluid. Thus, it becomes clear, why conventional approaches based on fluid dynamics have failed to explain the halo profiles so far.
The limits of the general constraining equation of the power-law index mentioned above are also determined with respect to this length scale. The mathematical results are given a physical interpretation by deriving the respective mass density profile from a joint, specific spatial configuration of the dark matter particles. As a result, we find that the mass density profile belonging to a set of freely moving particles confined to a volume much smaller than the scaling length shows the same behaviour as assumed in the core region of the heuristic power-law mass density profiles introduced in Section 2. The fact that the extremum power-law index belongs to a likelihood minimum, such that we expect the ensemble configuration to evolve away from it, 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 Section 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 Section 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. (28) 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.