Dark matter halo mass functions and density profiles from mass and energy cascade

Halo abundance and structure play a central role for modeling structure formation and evolution. Without relying on a spherical or ellipsoidal collapse model, we analytically derive the halo mass function and cuspy halo density (inner slope of −4/3) based on the mass and energy cascade theory in dark matter flow. The hierarchical halo structure formation leads to halo or particle random walk with a position-dependent waiting time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _g$$\end{document}τg. First, the inverse mass cascade from small to large scales leads to the halo random walk in mass space with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _g\propto m_h^{-\lambda }$$\end{document}τg∝mh-λ, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_h$$\end{document}mh is the halo mass and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}λ is a halo geometry parameter with predicted value of 2/3. The corresponding Fokker-Planck solution for halo random walk in mass space gives rise to the halo mass function with a power-law behavior on small scale and exponential decay on large scale. This can be further improved by considering two different \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}λ for haloes below and above a critical mass scale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_h^*$$\end{document}mh∗, i.e. a double-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}λ halo mass function. Second, a double-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document}γ density profile can be derived based on the particle random walk in 3D space with a position-dependent waiting time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _g \propto \Phi (r)^{-1} \propto r^{-\gamma }$$\end{document}τg∝Φ(r)-1∝r-γ, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi$$\end{document}Φ is the gravitational potential and r is the particle distance to halo center. Theory predicts \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =2/3$$\end{document}γ=2/3 that leads to a cuspy density profile with an inner slope of −4/3, consistent with the predicted scaling laws from energy cascade. The Press-Schechter mass function and Einasto density profile are just special cases of proposed models. The small scale permanence can be identified due to the scale-independent rate of mass and energy cascade, where density profiles of different halo masses and redshifts converge to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-4/3$$\end{document}-4/3 scaling law (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _h \propto r^{-4/3}$$\end{document}ρh∝r-4/3) on small scales. Theory predicts the halo number density scales with halo mass as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto m_h^{-1.9}$$\end{document}∝mh-1.9, while the halo mass density scales as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto m_h^{4/9}$$\end{document}∝mh4/9. Results were compared against the Illustris simulations. This new perspective provides a theory for nearly universal halo mass functions and density profiles.


INTRODUCTION
Within the standard ΛCDM (cold dark matter) cosmology [1,2,3,4], the formation of structures proceeds hierarchically with small structures coalescing into large structures in a "bottom-up" fashion.For systems involving long-range interaction, the formation of haloes of different sizes is necessary to maximize system entropy [5].Therefore, highly localized halo structures and their evolution are major features of ΛCDM model [6,7].As a counterpart of "eddies" in hydrodynamic turbulence, "haloes" are the building blocks in the flow of dark matter [8,9,10].Halo abundance and internal structure play a central role for modeling structure formation and evolution.These two quantities are also critical to understand the small scale challenges for ΛCDM when comparing model with observations [11,12,13,14].However, despite having been extensively studied over many decades, our understanding is still not entirely satisfactory.
First, the abundance of dark matter haloes is described by a halo mass function.The seminal Press-Schechter (PS) model allows one to predict the shape and evolution of mass function based on a density peak approach [15].This model relies on a threshold value of overdensity (  ) that can be obtained from the nonlinear collapse of a spherical over-density [16,17].Bond et al. provided an alternative derivation using an excursion set approach (EPS) that puts the theory on a firmer footing by removing the fudge factor in original PS model [18], which was further extend to excursion set with correlated steps [19,20,21].The PS model was further improved by Jedamzik with a formalism explicitly counting all cosmic materials to address the so-called "cloud-in-cloud" problem in density peak approach [22].Lee and Shandarin adopted Zeldovich approximation and extended the PS formalism to a non-spherical dynamical model [23].Other developments include combination of the peak and excursion set approaches [20], a moving barrier as a better density threshold [24], and more recent efforts on developing emulators of halo mass functions for a range of different cosmologies [25].
However, when compared to N-body simulations, both PS and EPS models overestimate the number of low-mass haloes and underestimate the number of massive haloes.There are also significant errors at high redshifts [26].Further improvement was achieved by computing the density threshold   for ellipsoidal collapse [27,28].In contrast to the spherical collapse where   is independent of halo mass, the ellipsoidal collapse leads to a mass-dependent overdensity threshold   .This modification (hereafter ST) considerably complicates the derivation but provides a better agreement with simulations.
Because of its simplicity, the PS-EPS-ST mass functions are still a very popular analytic model.However, the theoretical basis of this approach is at best heuristic.First, the derivation requires a threshold overdensity from a simplified (if not over simplified) collapse model (either spherical or ellipsoidal).Second, the linear density field is required to identify collapsed structures that is deeply in the non-linear regime.In principle, halo mass function should be an objective intrinsic property of self-gravitating collisionless system that is independent of any simplified (spherical or ellipsoidal) collapse models.In this paper, a different approach is taken to derive the halo mass function without resorting to any simplified models.This approach is based on the random walk of haloes in mass space, which is a direct result of inverse mass cascade in dark matter flow [10].
Next, the structure of haloes is described by the halo density profile that can be studied both analytically and numerically with N-body simulations [29,30].Since the seminal work of spherical collapse [17], the power-law density profile was derived under the self-similar approximation.The secondary in-fall model suggests a power-law density dependent on the initial density of the region that collapsed [31,32].High-resolution N-body simulations have shown nearly universal profile with a cuspy density shallower than isothermal profile at smaller radius and steeper at larger radius [33,34].For the cuspy inner density from N-body simulations, there seems no consensus on the exact value of the asymptotic logarithmic density slope .Since the first prediction of  = −1.0 in NFW profile [33], the inner density slope of simulated haloes have different values from  > −1.0 [35] to  = −1.2[36], and  ≈ −1.3 [37,38,39].In addition, there still lacks a complete understanding for the origin of nearly universal density profile [7].In this paper, similar to the halo random walk in mass space for halo mass function, a new approach is presented based on the particle random walk in real space, which provides a possible theory for nearly universal halo structures and density profiles.

EXISTING HALO MASS FUNCTIONS
For comparison with our mass function model, a brief overview of existing mass functions is presented here.The exact definition of mass function varies widely in the literature.The two widely used mass functions are defined as where ( ℎ , ) is the number density of haloes,  0 is the background density.Here   ( ℎ ) is the density fluctuation when density field is smoothed at mass scale  ℎ , which can be computed from the density power spectrum.When a normalized variable  =  2  / 2  ( ℎ ) is used, the third definition  () can be introduced such that the multiplicity mass function f(  , ) = 2  ().In this definition, the PS mass function reads The modified PS model (ST model) can be compactly written as: where the normalization condition requires: The best fitted parameters from simulation is  = 0.3222,  = 0.707, and  = 0.3 (hereafter ST1), while  = 0.3222,  = 0.75, and  = 0.3 was suggested by Sheth and Tormen [40] (hereafter ST2).Both models satisfy the normalization condition ∫ ∞ 0  ()  = 1.Many empirical mass functions were also proposed by fitting to the high-resolution simulation data.For example, a universal mass function by Jenkins etc. (hereafter JK) covers a wide range of different cosmologies and redshifts that is written as [41], where the threshold density   = 1.6865.Using a similar form of mass function to ST, Warren proposed (hereafter WR) [42]    () = 0.7234 It should be noted that these empirical mass functions might not satisfy the normalization constraint and can be difficult to extrapolate beyond the range of fit.
The other widely used empirical mass function by Tinker etc. was also calibrated from numerical simulations with haloes identified as isolated spherical overdensity masses.The range of halo mass is between 10 11 and 10 15 ℎ −1  ⊙ with redshift  ⩽ 2 [43].TK mass function reads or equivalently, where best fitted parameters  = 0.186,  = 1.47,  = 2.57 and  = 1.19 for haloes with a critical density ratio Δ  = 200.Table 1 summarizes different halo mass functions f(  , ) in Eq. ( 1).The double- mass function is analytically derived in Section 4.

MASS AND ENERGY CASCADE BETWEEN HALOES
To derive the halo mass function and density profiles, we first introduce the relevant context and background.In CDM cosmology, haloes are continuously merging with small structures (mass accretion).This facilitates an inverse mass cascade in halo mass space, i.e. a continuous mass transfer from small to large mass scales ("inverse") to allow hierarchical structure formation (see Fig. 1).To explain this, we first identify all haloes in entire system and then group them according to their mass  ℎ .In simulation, a clear definition of halo is required to identify these haloes.This definition is usually related to a critical density   from a simplified collapse model.At this step, we just treat haloes as existing objects without triggering a specific halo definition.In Fig. 1, halo of mass  ℎ merging with a single merger of mass  results in a new halo of mass  ℎ + .This causes a continuous mass flux from small to large scales along the chain of merging, i.e. an inverse mass cascade at a rate of   .
Next, the mass of entire halo group (  ) including all haloes of the same mass  ℎ is   =  ℎ  ℎ , where  ℎ is the number of haloes in that group.Now let's consider the most dominant and frequent merging, i.e. the merging with a single merger (or a single particle  of mass ) in Fig. 1, where  ℎ is the average waiting time of a given halo group, i.e. the average time interval between two subsequent merging events involving single mergers with any one halo in the same group.Therefore, the rate of mass transfer (or cascade) from mass scale  ℎ to scale  ℎ +  during the time interval  ℎ should be i.e. the entire halo mass  ℎ is transferred to a larger scale in a time interval  ℎ .This equals the rate of change for total mass in all haloes greater than  ℎ .Here  ℎ () is the total mass in all haloes,   ( ℎ , ) =   / 0 (see Eq. ( 1)) is the probability distribution of total halo mass  ℎ with respect to  ℎ .The integration gives the total mass in all haloes greater than scale  ℎ .The 'minus' sign stands for the "inverse" cascade from small to large scales.When self-gravitating collisionless system reaches a statistically steady state, this rate of mass transfer must be scale independent (i.e.  is independent of  ℎ ).If this is not the case, there would be a net accumulation of mass at some intermediate mass scale below  * ℎ .We exclude this possibility because we require statistical structures of haloes to be self-similar and scale free for haloes smaller than  * ℎ .This leads to the rate of mass cascade   independent of mass scale  ℎ up to a critical mass  * ℎ [10].Therefore, taking the derivative of Eq. ( 8) with respect to  ℎ leads to where   =  ℎ  ℎ is the halo group mass,   is mass of a single particle (mass resolution in N-body simulation).
Here the scale-independent   requires the halo group mass   ( ℎ , ) ≡   ( ℎ ) to be independent of time, i.e. a "small scale permanence" where the group mass   of different halo masses  ℎ and different redshifts  should collapse on to a common scaling law (Eq.(10) and Fig. 2).Once the statistically steady state is established, the rate of mass cascade   becomes scale-independent.The halo group mass   in propagation range becomes time independent due to scale-independent   .Mass is simply injected at the smallest scale (scale of single mergers), propagated to larger scales in propagation range ( ℎ <  * ℎ ), and consumed to grow haloes in deposition range ( ℎ >  * ℎ ).Halo group mass   ( ℎ ) is constant in time for haloes  ℎ <  * ℎ , and grows with time for haloes  ℎ >  * ℎ .Similarly, due to scale-independent energy cascade, the "small scale permanence" for halo density profile will be identified in Section 5 (Fig. 10).
To validate this concept, Fig. 2 presents results from large scale cosmological Illustris simulation (Illustris-1-Dark) [44].Illustris is a suite of large volume cosmological DM-only and hydrodynamical simulations.The selected Illustris-1-Dark is the DM-only simulation of 106.5Mpc 3 cosmological volume with 1820 3 DM particles for the highest resolution.Each DM particle has a mass around 7.6 × 10 6  ⊙ .The gravitational softening length is around 1.4kpc.Haloes in simulation were identified by a standard friends-of-friends (FoF) algorithm with linking length parameter b = 0.2 and halo center placed at the minimum of the gravitational potential of entire halo.Simulation has cosmological parameters of a total matter density Ω  = 0.2726, dark energy density Ω  = 0.7274 at  = 0, and a dimensionless Hubble constant ℎ = 0.704.
Next, if we focus on a given halo in a halo group, the waiting time   for that particular halo to merge with a single merger should be different and much greater than  ℎ (the waiting time for entire group).Here   is expected to be inversely proportional to the surface area of that halo.The larger surface area  ℎ , the more likely for Halo mass m h (M sun )  10)) at small mass scale (propagation range) with halo geometry parameter  ≈ 0.88.that halo to merge with a single merger, and the smaller waiting time   .Therefore, for haloes with a given mass  ℎ , this waiting time where  is a key halo geometry parameter.Intuitively,  ≈ 2/3 for large haloes (i.e. ℎ ∝  2/3 ℎ ).This is also equivalent to the waiting time   ∝ Φ −1 , where Φ ∝  ℎ / ℎ is the gravitational potential and  ℎ ∝  1/3 ℎ is the size of halo.The greater halo gravitational potential Φ, the larger velocity dispersion  2 from virial theorem (or halo temperature), the smaller waiting time   , and the more frequently halo merging with single mergers.Particle waiting time is dependent on its local potential.This will be used for deriving halo density profile in Section 6.
Depending on the number of haloes  ℎ in a given halo group, the two waiting times   and  ℎ are related to each other as Again, due to scale-independent rate of mass cascade   (not varying with  ℎ in propagation range), Eq. ( 10) requires the number of haloes for any given mass  ℎ , or equivalently a power-law group mass   =  ℎ  ℎ ∝  − ℎ at small mass scales, i.e. the small scale permanence in Fig. 2. In the same figure, we obtain  ≈ 0.88 for Illustris simulation and number of haloes in halo group  ℎ ∝  −1.9 ℎ that is in good agreement with other work [45].
To summarize, the mass cascade at statistically steady state involves two ranges, the propagation and deposition range.The propagation range for haloes with mass  ℎ <  * ℎ involves a sequence of merging with single mergers (the smallest structure) to simply propagate mass to larger scales.In this range, the rate of mass transfer   is independent of halo mass  ℎ and halo group mass   is constant in time.The deposition range ( ℎ >  * ℎ ) involves the consumption (deposition) of mass cascaded from scales below  * ℎ to grow haloes above  * ℎ (Fig. 1).Therefore, the inverse mass cascade can be described as: "Little halos have big halos, That feed on their mass; And big halos have greater halos, And so on to growth." In addition, haloes possess finite kinetic and potential energy.
Accompanied by the mass cascade, there exists a simultaneous energy cascade across haloes of different masses [46,47].The rate of energy cascade   ∝    2 / ℎ ∝ −  2 , where  2 is the mean kinetic energy of all particles in all haloes.The specific rate of energy cascade per unit mass (  < 0 for inverse energy cascade) can be estimated from the time variation of velocity dispersion  2 0 for all dark matter particles, where  0 ≈ 350/ from N-body simulation and  0 is the current age of universe [9].Therefore, similar to the mass cascade in propagation range, there exist an inverse (kinetic) energy cascade from small to large scales with a constant rate   .In this range of scales, the small scale structures evolve so fast and do not feel the slowly evolving large scale structures directly except through constant rate   .This description indicates that relevant quantities in this range of scales should be determined by and only by   ( 2 / 3 ), gravitational constant  ( 3 / •  2 ), and the relevant length scale r.By a simple dimensional analysis, the halo mass enclosed within  and corresponding halo density should follow the scaling [9] i.e. the 5/3 law and -4/3 law.These results can be demonstrated and confirmed by both N-body simulations (Figs. 12 to 15) and halo density profiles from random walk in Section 6 (Eq.( 30)).

DOUBLE-𝜆 HALO MASS FUNCTION
To derive halo mass function, the inverse mass cascade can be transformed into a halo random walk in mass space that mimics the random work of particles for diffusion problem.Just similar to the particle diffusion, we can derive the relevant Fokker-Planck equation and corresponding solution, from which halo mass function can be analytically solved.This is not just mathematically convenient, but reveals some fundamental aspects of halo mass function as an intrinsic property of self-gravitating collisionless system.As shown in Fig. 1, haloes are continuously migrating in mass space from one scale ( ℎ ) to neighboring scale ( ℎ + ) by merging with single mergers.This leads to a probability distribution to find a halo at a given mass.The waiting time (or jumping frequency) for a given halo to migrate from a given mass  ℎ to neighboring mass  ℎ +  is   in Eq. (10).Different from the standard random walk with a constant waiting time, the halo waiting time   is dependent on the mass of halo, i.e. a position-dependent   (Eq.( 10)).For halo with a given mass  ℎ , the waiting time   ∝  − ℎ , where  is a key halo geometry parameter we discussed.
First, the random walk of haloes in mass space describes the stochastic variation of the mass of a given halo due to continuous merging with single mergers of mass .Following the Langevin equation, we can write a stochastic equation for halo mass  ℎ [10]  ℎ () where /  represents the average rate of mass change.For a powerlaw waiting time   ∝  − ℎ , we find the position-dependent diffusivity should take the form of Here  0 () is a proportional constant for diffusivity   .The white Gaussian noise  () satisfies the covariance ⟨ () ( with a zero mean ⟨ ()⟩ = 0. Equation ( 13) describes the stochastic evolution of halo mass  ℎ with a waiting time   ( ℎ ) ∝  − ℎ .Second, in Stratonovich interpretation [48], the Langevin equation (Eq.( 13)) yields to a distribution function  ℎ ( ℎ , ) satisfying the Fokker-Planck equation (resembling particle diffusion) which describes the evolution of probability function  ℎ for halo mass  ℎ in mass space.Obviously, the halo mass function   ( ℎ , ) is exactly the distribution function  ℎ , i.e.   ≡  ℎ .Finally, solution to Eq. ( 15), i.e. the halo mass function, is a stretched Gaussian with an exponential cut-off for large  ℎ and a power-law behavior for small  ℎ , The mean square displacement in mass space is where  * ℎ () is the critical mass scale and  0 is just a proportional constant.With the exponent of 1/(1 − ) ⩾ 1 in Eq. ( 17), it is clear that the random walk of haloes in mass space is of a super-diffusion nature.Now   ( ℎ , ) (Eq.( 16)) can be rewritten in terms of where the dimensionless constant The time dependence of   is absorbed into  * ℎ .Intuitively,  ≈ 2/3 for large haloes in deposition range with low concentration, whose central structures are still dynamically adjusted due to fast mass accretion.While for small haloes with high concentration (propagation range), the mass accretion is slow and inner structure is stable [49].These small haloes can be treated as fractal objects with a fractal surface dimension  ℎ ⩽ 3. The geometry parameter  =  ℎ /3 can be greater than 2/3 (see Fig. 2).These high concentration low mass haloes are usually found in denser environments [50].The denser environment might lead to a rougher halo surface and higher surface fractal dimension  ℎ .Therefore, two different  (i.e.double-) are required for two ranges (propagation range with  ℎ <  * ℎ and deposition range with  ℎ >  * ℎ ) due to different halo properties and surrounding environments.The single- halo mass function in Eq. ( 18) can be naturally generalized to a double- halo mass function with  1 and  2 for propagation and deposition ranges, respectively.Therefore, the double- mass function reads By introducing variable  = ( ℎ / * ℎ ) 2/3 , the three parameter double- mass function can be finally written as, where model parameters  and  have clear physical meaning.Both are related to halo geometry parameters  1 and  2 as, Clearly, Eq. ( 21) reduces to the Press-Schechter (PS) mass function if  1 =  2 = 2/3 and  0 = 1/2.However, the derivation of double- mass function does not rely on any collapse model (spherical or ellipsoidal).The critical overdensity   from collapse model is not required in this formulation.In simulation, haloes are usually defined using the critical overdensity   to compute the halo mass function.The derivation of double- mass function of Eq. ( 21) does not depend on the exact definition of halo.Different definitions of halo in simulation might affect both halo mass  ℎ and the critical mass  * ℎ , but not the ratio  = ( ℎ / * ℎ ) 2/3 , and therefore not the double- halo mass function.More importantly,  1 =  2 = 2/3 or  =  = 1 is a natural result of current theory.This formulation reveals that the halo mass function in the form of Eq. ( 21) is an intrinsic property of self-gravitating collisionless dark matter system that is independent of spherical or ellipsoidal collapse models.
The halo geometry exponent  has a fundamental meaning to relate halo surface area (or effective mass accretion area) to its mass.The cosmology and redshift dependence of  1 and  2 can be systematically studied by fitting the model to the simulation data of different cosmologies, similar to the study in Bocquet et al. [25] and Euclid Collaboration et al. [51].
Alternatively, similar to the scale radius   for halo density where logarithmic density slope is -2, we may introduce a scale mass  ℎ where logarithmic slope  ln(   )/ ln( ℎ ) = −1 such that  ℎ =  (2 0 ) 3/(2 )  * ℎ from Eq. (20).With a new scaled variable ν = ( ℎ / ℎ ) 2/3 , mass function in Eq. ( 21) can be further simplified with  and  as the only two parameters To validate the derived double- mass function, we presents results from Illustris simulation (Illustris-1-Dark) [44].Figure 3 presents the halo mass function f(  , ) in Eq. ( 1).The best fit of double- mass function to the simulation data at all z gives values of  0 = 1.162,  = 0.365, and  = 1.185 (Fig. 3), which leads to  1 = 0.856 and  2 = 0.605 from Eq. ( 22) for the propagation and deposition ranges, respectively.This leads to a slope of − 1 − 1 ≈ −1.9 for halo number density ( ℎ , ) ∝  −1.9 ℎ (Eq.( 10)), in very good agreement with Fig. 2 and other work [45].Compared to predicted value of  = 2/3 for matter dominant universe, the effect of dark energy in Illustris simulations seems to enhance the value of  1 and decrease the value of  2 , reflecting the changes in environments and halo properties due to the presence of dark energy and accelerated expansion.
The PS mass function overestimate the mass in small haloes and underestimates the mass in large haloes.The JK mass function matches simulation for large mass haloes with large deviation for small haloes.The fitted WR mass function does not satisfy the normalization condition, where ∫ ∞ 0    () diverges.The WR mass function also deviates at small mass with a finite limit f( −   1) with simulation results at  = 0, 4, 8, and12, as a function of halo mass  ℎ .Relative errors of different mass functions when compared to binned simulation data are also presented in the bottom plots.Similar conclusions can be obtained from these plots, where WR, ST, TK and double- mass functions agree with simulation at lower redshift.Double- mass function is slightly better at higher redshifts  = 8 and 12.

MASS SCALE 𝑚 *
ℎ AND SMALL SCALE PERMANENCE The inverse mass cascade and halo mass function (Eq.( 20)) require a critical halo mass scale  * ℎ that can be related to halo velocity dispersions from virial theorem where  2  ( ℎ ) is the velocity dispersion of all DM particles in a halo with a given mass  ℎ , which represents the temperature of that halo.Here ⟨⟩ represents the average for all haloes in the same group with same mass  ℎ .In addition,  2 ℎ =  ( ℎ ) is the dispersion (variance) of halo velocity  ℎ (the mean velocity of all particles in the same halo) for all haloes in the same group, where  2 ℎ represents the temperature of halo group that is relatively independent of halo mass  ℎ [5,47].Figure 8 presents an example of the variation of ⟨   21)) predicts less mass in larger haloes and slightly better agrees with the simulation.Halo mass m h (M sun )  24)).
compute the critical mass  * ℎ for other redshifts.The variation of  * ℎ with the scale factor  is presented in Fig. 9.In linear regime,  * ℎ ∝  3 is expected, while in nonlinear regime  * ℎ ∝  3/2 [10].With halo mass function in Eq. ( 18) and the small scale permanence for   in Eqs. ( 9), (10), and Fig. 2, the halo group mass   =  ℎ   (  is particle mass) should satisfy such that the total mass in all haloes  ℎ () ∝  * ℎ 1− when statistically steady state is established in the nonlinear regime.With  = 2/3 for  ℎ =  * ℎ ,  ℎ () ∝  1/2 is expected.The time variation of total halo mass  ℎ is also presented in Fig. 9.
Next, similar to the small scale permanence for group mass   in Fig. 2, we will present the small scale permanence for halo density profile.From the scaling laws due to energy cascade, the density scaling   ∝  −4/3 is proposed in Eq. ( 12), which already hints the small scale permanence.To demonstrate this concept, the density profiles for haloes with a critical mass  * ℎ at different redshifts are studied first.In Illustris-1-Dark simulation, all haloes with mass between 10 ±Δ  * ℎ are identified at different redshifts  with Δ = 0.1.The spherical averaged density profile is computed for every halo.The density profile for haloes with critical mass  * ℎ is computed as the average density profile for all haloes with mass between 10 ±Δ  * ℎ .Figure 10 presents the time evolution of halo density profiles for haloes with critical mass  * ℎ ().The small scale permanence from energy cascade can be clearly demonstrated as the density profiles for haloes with critical mass at different redshifts all collapse onto the predicted density scaling (blue solid line  ℎ ∝  −4/3 ) on small scales.Finally, if gravity is the only interaction and dark matter is fully collisionless and cold, extending the established scaling in Fig. 10 to the smallest length scale and and the earliest time (or the highest ) might be able to identify dark matter particle mass, size, lifetime, and many other properties [46].ℎ at different redshifts  collapse at small scale  onto the predicted density scaling (-4/3 law with  ℎ ∝  −4/3 ) from the theory of energy cascade (solid blue line from Eq. ( 12)).

DOUBLE-𝛾 HALO DENSITY PROFILE
The halo density profile can be analytically derived based on a similar idea as deriving halo mass function.Within CDM paradigm, the formation of structures starts from the gravitational collapse of small scale density fluctuations and proceeds hierarchically such that small structures coalesce into large structures in a "bottom-up" fashion.The halo structure is formed hierarchically through a series merging with smaller structures (dominantly with single mergers in Fig. 1).Now let us follow the mass accretion history of a given halo in Fig. 11, where halo mass   ≡   () (or halo size  ≡  (), the radius enclosing mass   ) continuously varies with time from 0 to   (or size from 0 to  ).The mean waiting time of every merging with a single merger  has a simple scaling as   ∝  −  ∝ Φ −1 , where  is a geometry parameter (see Eq. ( 10)) and Φ() ∝   / is the gravitational potential at .In 3D space, halo size  can be related to the position   of merger  as  = √   •   .Since both halo mass   () and Φ() can be related to size  (), the waiting time   should also be a function of  (), which means a varying waiting time dependent on the particle distance  to halo center where  is an exponent for -dependence of waiting time   , which can be related to the slope of density profile (see Eq. ( 30)).Since haloes are formed by sequential merging, every DM particle in any halo was a single merger at the time they joined that halo.That particle starts to continuously perform a 3D random walk with a position-dependent waiting time   dependent on its local potential Φ or  (Eq.( 26)) right after the merging, where Φ() is determined by the total enclosed mass within .In this regard, halo random walk in mass space is consistent with the particle random walk in 3D space.The random walk of DM particles has a position dependent waiting time   ∝ Φ() −1 ∝  − , where  = √   •   is the distance to halo center.The waiting time is also dependent on the local potential Φ(), or from virial theorem, the velocity dispersion  2 that represents the local temperature.Since energy cascade theory predicts the 5/3 law for mass scaling   ∝  5/3 for the inner region of virialized haloes (see Eq. ( 12)), we have potential Φ() ∝   / ∝  2/3 such that  = 2/3 from Eq. (26).A position dependent waiting time   () is an important feature for hierarchical formation of halo structure.A longer waiting time   () at small  means a more stable core region than the outer region.
Finally, the particle distribution resulting from this positiondependent random walk in 3D space gives rise to the halo density, as shown in Fig. 11.Therefore, to find the halo density profile, we need to derive the particle distribution function due to the random walk in 3D space with   () ∝  − .The 3D particle random walk can be described by a Langevin equation for particle position   (similar to Eq. ( 13) for halo random walk in mass space), Due to position-dependent waiting time   (), the positiondependent diffusivity reads where  0 () is a proportional constant.The smaller , the smaller diffusivity or longer waiting time, and the higher particle density.In Itô convention, the 3D Fokker-Planck equation in Cartesian coordinate can be directly obtained for particle distribution function   ( , ) ( = 1, 2, 3 for Cartesian coordinates), The corresponding solution of Eq. ( 29) in spherical coordinate is Since the distribution function   (, ) is equivalent to halo density, we find that the parameter  is half of the density slope at small .
From this insight, assume  is unknown, we can predict the value of  as follows: Since the waiting time   ∝ Φ() −1 ∝  − , halo density should scale as   ∝  −2 from Eq. (30).The halo mass enclosed in  scales as   ∝    3 ∝  3−2 .The local potential at  should scale as Φ() ∝   / ∝  3−2−1 .The waiting time of particle at  should satisfy Eq. ( 26) that requires 3 − 2 − 1 =  such that  = 2/3 and the density slope 2 = 4/3.It should be noted that the random walk theory for halo structure formation confirms the -4/3 law (  ∝  −4/3 ) predicted by the energy cascade theory in Eq. (12).Predictions are tested against simulations in Figs. 12 to 15.Similar to halo mass function (Eq.( 20)), the exponent  can be different in two different ranges, i.e. the power law below the scale radius   and the exponential decay above   .Using two different  for -dependence of waiting time   () ∝  − , i.e.  1 and  2 for two different ranges, based on the single- distribution in Eq. ( 30), the double- distribution reads Introducing the conventional scale radius   () where the logarithmic slope of   (, ) equals -2, we should have Substituting Eq. ( 32) into Eq.( 31) and introducing a dimensionless spatial-temporal variable  = /  (), distribution function reads Finally, the two parameter particle distribution function can be written as (with a similar form as mass function in Eq. ( 23)) where two dimensionless parameters  and  are -2 10 -1 10 0 10 1 10 2 10 3 10 4 r (kpc)  12)) for halo density is presented as the solid blue line.The double- density model (Eq.was also plotted for all haloes as dashed lines. The time variation of the distribution function is absorbed into the scale radius   ().The double- distribution function reduces to the Einasto profile with  = 2.The cumulative distribution in spherical coordinate can be easily obtained as, where Γ(, ) is an upper incomplete gamma function.So far we provide physical interpretation and a possible theory for halo density.The general density profile can be finally written as where   () is the density at scale radius   .Simulated haloes were found to have different density slopes in different simulations as discussed in Section 1.This might be due to the different radial flow and mass accretion rate in these haloes, whose density profile can be modelled by the general solution in Eq. ( 37) [9].On small scale, virialized haloes are incompressible with vanishing (proper) radial flow [54].For fully virialized haloes with vanishing radial flow, we would expect -4/3 law for inner density with 2 1 = 4/3, which is consistent with the limiting density slope in Eq. ( 12).Combining Eq. ( 37) with / = 2/3 leads to density profile that is consistent with the prediction from energy cascade in Eq. ( 12), The small scale permanence for halo density in Fig. 10 becomes where   is an amplitude parameter of halo density,   =  is a shape parameter of density profile, and   is the scale radius.
To validate the proposed density model in Eq. ( 38), spherical averaged density profile was first obtained for all haloes with given mass in a range of 10 ±Δ  ℎ at different redshifts .Next, we obtained the average halo density profile for all haloes in the same range at same redshift.The radial flow in these haloes might be cancelled out  12)) for halo density is presented as the solid blue line.for comparison, the double- density model (Eq.( 38)) was also plotted as dashed lines.Model fits better for halo density at higher redshift.The asymptotic density slope −4/3 at small  can be identified.12)) for halo density is presented as the solid blue line.The double- density model (Eq.( 38)) was also plotted as dashed lines.The asymptotic density slope −4/3 at small  can be identified.
after this averaging such that the averaged halo density can be better described by Eq. ( 38) with an inner slope of 2 1 = 4/3.Figures 12 to 15 present the halo density profiles of different halo mass  ℎ at different redshifts  from Illustris dark matter only simulations: Illustris-1-Dark (solid lines), where Δ is selected to be 0.1.The double- density model (Eq.( 38)) was also used to fit all haloes and plotted as dashed lines in these figures.The best-fit model parameters   ,   and   can be obtained for different halo mass  ℎ and redshifts  (as presented in Figs.16 to 18).The double- density model provides a reasonably well fit to all haloes at all redshifts, with slightly better fit at higher redshift in a matter-dominant universe.
Figure 16 presents the variation of amplitude parameter   with the dimensionless parameter  defined in Eq. (24).As expected, the amplitude parameter   ∝  2/3 increases with halo mass  ℎ at  12)) for halo density is presented as the solid blue line.The double- density model (Eq.( 38)) was also plotted as dashed lines for comparison.
Figure 17 presents the variation of shape parameter   with .The shape parameter   is relatively independent of parameter  at low redshift .It varies in a small range between 1 and 3 and slightly decreases with halo mass  ℎ , which corresponds to a range of  2 = 2/3 for large haloes and  2 = 0 for small haloes with  1 = 2/3 (see Eq. ( 35)).In the range  >   , the potential Φ is relatively independent of  due to exponential decay of density.Therefore, the waiting time becomes less dependent on  in this range with  2 ⩽  1 .Table 2 lists relevant values of  and  in different ranges.Figure 18 presents the variation of the best fitted scale radius   with  at different redshifts , where   increases with  with an approximate scaling of   ∝  1/2 .In summary, the amplitude parameter   is related to the rate of cascade  in haloes (Eq.( 40)), while the shape parameter   is related to the parameter  (Eq.( 35)), i.e. the position dependence of waiting time   ∝  − .
It would be also interesting to compare the density profile obtained in this work with the Einasto and NFW profiles.Figure 19 presents the comparison for small (10 8.5  ⊙ ) and large haloes (10 13  ⊙ ) at redshift  = 0 (haloes in Fig. 12).These density profiles include: 1) the general double- profile in Eq. (37) with  and  being independent; 2) the Einasto profile with  = 2 in Eq. (37); 3) the double- profile with  = 2/3 in Eq. (37) (or Eq. ( 38)) for fully virialized haloes; and 4) the standard NFW profile.Bottom plots present the relative errors between these density profiles and simulation results.As expected, the general double- profile provides the best fit of simulated halo density, compared to NFW profile.The double- profile with  = 2/3 (Eq.( 38)) provides a slight better fit than Einasto profile for small haloes, and a much better fit for large haloes.Finally, additional tests for different halo definitions and cosmologies should be very helpful to include data from simulations other than Illustris series.In this case, parameters in halo mass function and density models (Eqs.( 21) and ( 37)) need to be fitted for different cosmologies.From this study, we can find how model parameters (halo parameters  and ) vary with different cosmologies, which will require extensive work in future study.Here a quick test of double- density for some simulated haloes in the literature was presented.Figure 20 provides the best fit by the general model in Eq. ( 37) for these simulated haloes.Since the analytically derived double- profile reduces to Einasto profile for  = 2/3, the general double- profile is expected to provide a better fit than Einasto profile for all simulated haloes.

CONCLUSION
In this paper, a simple theory was presented for halo mass function and density profile.The small scale permanence is proposed for halo group mass   and halo density profile  ℎ due to scale-independent rate of mass and energy cascade (Figs. 2 and 10).Both halo mass function and halo density profile can be analytically derived based on this simple theory.The position-dependent waiting time   ∝  − ℎ leads to an analytical mass function modelled by a stretched Gaussian with a power-law behavior on small scale and exponential decay on large scale (Eq.( 18)).This can be further improved by considering two different values of  in propagation and deposition ranges, i.e. a double- mass function in Eq. (21).Similarly, a double- halo density profile is proposed based on the particle random walk in 3D space with a position-dependent waiting time   ∝  − (Eq.( 37)).The predicted value of  = 2/3 leads to a cuspy density profile with an inner slope of -4/3, consistent with the energy cascade theory (Eq.( 12)).The Press-Schechter mass function and Einasto profile are just special cases of the proposed model.Models were compared and validated against the Illustris simulations.Future work will involve additional tests for proposed models in different cosmologies.

1 Introduction 2
Existing halo mass functions 3 Mass and energy cascade between haloes 4 Double- halo mass function 5 Mass scale  * ℎ and small scale permanence 6 Double- halo density profile 7 Conclusion

) ℎ − 1 11 −Figure 1 .
Figure 1.Schematic plot of the inverse mass cascade for hierarchical structure formation.Halo of mass  ℎ merges with single merger (free DM particles of mass ) to cause the mass flux into haloes on larger scales  ℎ +  and the next merging along the chain.This facilitates a continuous mass cascade from small to large scales.A scale-independent mass flux   is expected for haloes in the mass propagation range (<  * ℎ ).Mass cascaded from small scales is simply propagated in the propagation range and consumed to grow haloes with mass >  * ℎ in the deposition range.

Figure 3 .
Figure 3.Comparison between different halo mass functions f(   , ) and simulation at different redshift z.The PS mass function overestimate the mass in small haloes and underestimates the mass in large haloes.The fitted JK mass function matches simulation only in a given range with large deviation for small mass haloes.The WR mass function deviates at small mass with a limit f(  −1  → 0, ) = −1.695.The double- mass function (Eq.(21)) with best fitting parameters  0 = 1.162,  = 0.365, and  = 1.185 (or  1 = 0.856 and  2 = 0.605) matches the simulation and is slightly better than ST mass functions at large halo mass.

Figure 4 .
Figure 4. Comparison of mass functions with Illustris-1-Dark simulation (solid blue) at z=0.The PS mass function overestimates mass in small haloes.The fitted JK mass function matches simulation only in a given range.The double- mass function (Eq.(21)) matches both simulation and the ST and WR mass functions at z=0.Bottom plot presents the relative errors between simulation and different mass functions.

Figure 5 .
Figure 5.Comparison of mass functions with Illustris-1-Dark (solid blue) at z=4.The simulation results agree with all mass functions except PS.Double- mass function (Eq.(21)) predicts a slightly lower mass in larger haloes.

Figure 7 .
Figure 7.Comparison of mass functions with Illustris-1-Dark simulation (solid blue) at z=12.Compared to other mass functions, the double- mass function (Eq.(21)) predicts less mass in larger haloes and slightly better agrees with the simulation.

2 Figure 9 .
Figure 9.The variation of critical halo mass  * ℎ and total mass  ℎ in all haloes with scale factor . Two regimes can be identified.In the linear regime  * ℎ ∝  3 .In nonlinear regime  * ℎ ∝  3/2 and  ℎ ∝  1/2 , where statistically steady state is established with a scale-independent rate of cascade.Density profiles of haloes with critical mass  * ℎ are presented in Fig. 10.

Figure 10 .
Figure 10.The evolution of halo density profiles for haloes with critical mass  * ℎ ().Figure demonstrates the small scale permanence, i.e. the density profiles for haloes with mass  *ℎ at different redshifts  collapse at small scale  onto the predicted density scaling (-4/3 law with  ℎ ∝  −4/3 ) from the theory of energy cascade (solid blue line from Eq. (12)).

Figure 11 .
Figure 11.Schematic plot of the growth of a given halo in both mass   and size  via continuous merging with single mergers, where the waiting time   (  ) ∝  −  .Every merging event corresponds to a single move of particle  in a random walk process, where the waiting time   ( ) ∝  − .Single mergers continuously join halo and perform 3D random walk.Particle distribution from 3D random walk gives rise to the halo density.

Figure 14 .
Figure14.density profiles for different halo mass  ℎ at  = 8 (solid lines).The predicted scaling law (Eq.(12)) for halo density is presented as the solid blue line.The double- density model (Eq.(38)) was also plotted as dashed lines.The asymptotic density slope −4/3 at small  can be identified.

2 Figure 18 .
Figure 18.The variation of scale radius   for halo density with  at different redshifts .The scale radius increases with  as   ∝  1/2 .

Table 1 .
Different Halo Mass Functions f(   , ) [47]halo velocity dispersions ⟨  2  ⟩ ( ℎ ) and  2 ℎ at  = 8 from Illustris-1-Dark simulation.The two velocity dispersions represent the temperature of haloes and temperature of halo groups[47].The large fluctuation at large mass scale is due to fewer massive haloes. Her ⟨  2  ⟩ ∝  is relatively independent of  ℎ .The critical halo mass  * ℎ ( = 8) = 9 × 10 10  ⊙ is found by setting ⟨  2  ⟩ ( * ℎ ) =  2 ℎ (Eq.( Figure Halo density profiles for different halo mass  ℎ at  = 0 (solid lines).The predicted scaling law (Eq.( Halo density profiles for different halo mass  ℎ at  = 4 (solid lines).The predicted scaling law (Eq.( Halo density profiles for different halo mass  ℎ at  = 12 (solid lines).The predicted scaling law (Eq.( Figure 16.The variation of amplitude parameter   for halo density with the dimensionless parameter  at different redshifts .In principle,   increases with halo mass  ℎ .This is related to the waiting time   ∝  − ℎ .fixed redshift or decreases with time at fixed mass  ℎ .The mass cascade across haloes is accompanied by a simultaneous energy cascade across haloes.The rate of cascade is independent of mass scale for group of haloes of the same mass.For individual haloes with mass  ℎ <  * ℎ , the rate of energy cascade  in these haloes is smaller due to the longer waiting time   ∝  − ℎ .The effective rate of energy cascade  in individual haloes is inversely proportional to   , ( ℎ , ) =  ℎ / * ℎ    =  3/2   .

Table 2 .
The variation of shape parameter   for halo density with  at different redshifts .The shape parameter   varies in a small range between 1 and 3 and slightly decreases with halo mass  ℎ .Halo parameters  and  from theory and simulation