{'fraction_non_alphanumeric': 0.0495851244626612, 'fraction_numerical': 0.01542037388783365, 'mean_word_length': 4.345758183032732, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Stochastic Theory of the Hierarchical Clustering I. Halo Mass Function

We present a new theory for the hierarchical clustering of dark matter (DM) halos based on stochastic differential equations, that constitutes a change of perspective with respect to existing frameworks (e.g., the excursion set approach); this work is specifically focused on the halo mass function. First, we present a stochastic differential equation that describes fluctuations in the mass growth of DM halos, as driven by a multiplicative white (Gaussian) noise dependent on the spherical collapse threshold and on the power spectrum of DM perturbations. We demonstrate that such a noise yields an average drift of the halo population toward larger masses, that quantitatively renders the standard hierarchical clustering. Then, we solve the Fokker-Planck equation associated to the stochastic dynamics, and obtain the Press&Schechter mass function as a (stationary) solution. Moreover, generalizing our treatment to a mass-dependent collapse threshold, we obtain an exact analytic solution capable of fitting remarkably well the N-body mass function over a wide range in mass and redshift. All in all, the new perspective offered by the theory presented here can contribute to better understand the gravitational dynamics leading to the formation, evolution and statistics of DM halos across cosmic times.


INTRODUCTION
The halo mass function, namely the statistics describing the number of dark matter (DM) halos per unit comoving volume as a function of halo mass and redshift, is a fundamental quantity in astrophysics and cosmology (see textbooks by Mo et al. 2010 andCimatti et al. 2020). For example, it is a basic ingredient to develop sensible galaxy formation and evolution models (see reviews by Silk &Mamon 2012 andNaab &Ostriker 2017), and it is routinely exploited in cosmological studies relying on the abundance and clustering of collapsed objects and of large-scale structures (see reviews by White 2012 andWechsler &Tinker 2018).
Clearly, the halo mass function can be estimated via high-resolution, large-volume, cosmological N -body simulations (see Sheth & Tormen 1999;Jenkins et al. 2001;Warren et al. 2006;Tinker et al. 2008;Crocce et al. 2010;Bhattacharya et al. 2011;Watson et al. 2013). However, given the natural limits on resolution, computational time, and storing capacity, it can be probed only in limited mass and redshift ranges. Moreover, the results of simulations depend somewhat on the algorithm used to identify collapsed halos (e.g., friend-of-friend vs. spherical overdensity), and on specific parameters related to the identification of isolated objects (e.g., the linking length). On the other hand, to estimate the halo mass function from observations is even more challenging (e.g., Castro et al. 2016;Dong et al. 2019;Sonnenfeld et al. 2019;Li et al. 2020), given the statistical and systematic uncertainties that arise when linking the observable quantities to the halo mass. Therefore, a deep theoretical understanding on how the halo mass function is originated from first principles is of crucial importance.
The modern theoretical framework to address the issue was born with the seminal work by Press & Schechter (1974); these authors were able to compute an analytic expression for the halo mass function by prescribing that a halo would collapse if it resided within a sufficiently overdense region of the initial (Gaussian) perturbation field. Given that the overdensity around a spatial location depends on scale, they recognized that the halo abundance is simply related to the mass fraction in the density field, smoothed on different scales, which is above a critical threshold for collapse. A drawback of this approach is the so called 'cloud-in-cloud' problem, i.e., attention must be paid not to double count overdense regions embedded within a larger collapsing perturbation; in other words, one has to consider only the perturbations that overcome the threshold on a given smoothing scale, but not on a larger one.
The problem was solved with the development of the excursion set framework by Bond et al. (1991), which still nowadays constitutes the standard theory. This envisages the overdensity around a given spatial location to execute a random walk when considered as a function of the smoothing scale; if the smoothing is performed with a sharp filter in Fourier space, the walk is Markovian. The collapse threshold here plays the role of a barrier, and the halo mass function is related to the distribution of first crossing, i.e., the probability that a walk crosses the barrier for the first time on a specific scale. In the original theory, the collapse threshold was gauged on the spherical collapse model of DM perturbations (Gunn & Gott 1972), and as such it was assumed to be independent on halo mass; subsequent developments adopted a mass-dependent threshold, inspired from the ellipsoidal collapse model (see Sheth & Tormen 2002), in order to better reproduce the results of N -body simulations.
Despite this rich literature focused on the theoretical foundations and a number of undoubtable successes in practical applications, the excursion set framework is known to hide some pitfalls and drawbacks: no exact analytic expression of the mass function for a general mass-dependent collapse threshold is known, but for very simple shapes (see Zhang & Hui 2006;Lapi et al. 2013); the merging kernel associated to the excursion set theory is not symmetric, and this causes a mathematical inconsistency or at least an ambiguity in the definition of the merger rates (see Benson et al. 2005;Neistein & Dekel 2008;Zhang et al. 2008); even adopting a massdependent collapse threshold and other refinements, the excursion set formalism struggles to reproduce the halo progenitors' distributions extracted from N -body simulations (e.g., Parkinson et al. 2008;Jiang & van den Bosch 2014); the relation between the probability of first upcrossing and the mass function, which is at the heart of the excursion set framework, is correct on statistical grounds but cannot be strictly true for individual mass elements (see discussion by Mo et al. 2010, their Sect. 7.2.2b).
In this paper we submit a new theory of the hierarchical clustering and halo mass function based on stochastic differential equations in real space, that constitutes a change of perspective with respect to the excursion set formalism. First, we invent a stochastic differential equation that describes fluctuations in the mass growth of DM halos, as driven by a multiplicative white (Gaussian) noise dependent on the spherical collapse threshold and on the power spectrum of DM perturbations; in this approach it is the mass (or mass variance) in a given region of the Universe to perform a (Markovian) random walk as a function of cosmic time. In Sect. 2 we demonstrate that the noise yields an average drift toward larger masses, that quantitatively renders the standard hierarchical clustering. Then, in Sect. 2.1 we solve the Fokker-Planck equation associated to the stochastic dynamics, and obtain as a solution the Press & Schechter mass function; in Sect. 2.2 we point out that the solution is stationary when the original equation is written in convenient variables.
In Sect. 3 we introduce a minimal modification of the stochastic equation in terms of a mass-dependent collapse threshold. Using a parametric shape analogous to that adopted in the excursion set framework, we obtain a closed-form analytical solution of the associated Fokker-Planck equation. Remarkably, such a solution has a shape similar to the empirical fitting formula introduced since Sheth & Tormen (1999); moreover, for specific values of the parameters describing the mass dependence of the collapse threshold, our result reproduces remarkably well the N -body mass function over an extended range of masses and redshifts.
As an aside, in Sect. 4 we explore how to generalize our framework when a colored rather than a white noise is considered, so as to enforce a non-Markovian evolution. Adopting for definiteness a multiplicative Ornstein-Uhlenbeck noise and a constant threshold for collapse, we are able to solve exactly the corresponding Fokker-Planck equation and obtain a closed form solution. With respect to the white noise case, the redshift evolution of the mass function is found to be modified somewhat, in a fashion depending on the correlation time characterizing the noise. Finally, in Sect. 5 we summarize our findings and envisage possible outlooks.
Throughout this work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration 2018) with rounded parameter values: matter density Ω M ≈ 0.3, dark energy density Ω Λ ≈ 0.7, baryon density Ω b ≈ 0.05, Hubble constant H 0 = 100 h km s −1 Mpc −1 with h ≈ 0.7, and mass variance σ 8 ≈ 0.8 on a scale of 8 h −1 Mpc. The most relevant expressions are highlighted with a box.

A STOCHASTIC EQUATION FOR THE
HIERARCHICAL CLUSTERING Our proposal is to capture the essence of the hierarchical clustering via the following nonlinear stochastic differential equation: or equivalently, in terms of mass 5 ( Here η(t) is a Gaussian white noise (physical dimension 1/ √ time) with ensemble-average properties η(t) = 0 and η(t) η(t ) = 2 δ D (t−t ), where δ D is the Dirac deltafunction (the factor 2 is only a convention and clearly it could be reabsorbed into the multiplicative term). The quantity δ c (z) = δ c0 D(0)/D(z) is the critical threshold for collapse extrapolated from linear perturbation theory; in a flat Universe (see Mo et al. 2010;Weinberg 2008), one can use the approximations δ c0 (1), yielding the evolution of the mass variance ln σ (left y−axis) and of the mass M (right y−axis) as a function of cosmic time t; the initial condition M ≈ 10 4 M at z ∼ 10 (corresponding to t ∼ 0.5 Gyr) has been adopted. The grey lines are 30 randomly chosen evolutionary tracks, while the red line and shaded area illustrate the median and the quartiles over 3000 realization of the noise; the dashed blue line shows the evolution of the spherical collapse threshold δc(t).
Finally, σ is the mass variance filtered on the mass scale M : here P (k) is the power spectrum of density fluctuation, andW 2 M (k) is the Fourier transform of a window function whose volume in real space encloses the mass M ; for standard cold dark matter power spectra (e.g., Bardeen et al. 1986), σ(M ) is an inverse, convex, slowly-varying function of M . Note that in the present theory the relation σ(M ) between σ and M is purely deterministic, but both M (t) and σ(t) = σ(M (t)) are to be considered stochastic variables that fluctuate over cosmic time t under the influence of the noise. We stress the change of perspective with respect to the standard excursion set formalism: in the latter the overdensity field δ(σ) executes a random walk as a function of the mass variance σ, which plays the role of a pseudo-time variable; here the mass M (t) or the mass variance σ(M (t)) are themselves stochastic variables, undergoing a Markovian evolution as a function of (real) time t. Note that in the excursion set approach the choice of the filter function in Eq. (3) has a crucial impact, since the random trajectories δ(σ) are Markovian only when a sharp filter in Fourier space is adopted; in the present theory, assuming a different filter function (e.g., Gaussian or top-hat in real space) changes only the deterministic relation σ(M ) but has otherwise no effect on the Markovianity of the stochastic processes M (t) and σ(M (t)).
The rationale naively followed to invent the above Eq. (1) is simple. On the left hand side it appears the time derivative of an adimensional function of the mass that incorporates the power spectrum and the filtering scale; in choosing ln σ we have been inspired by a number of N -body simulations (e.g., Zhao et al. 2009), that suggest the mass growth of halos to be easily described in terms of such a quantity. On the right hand side it appears a stochastic driving η(t), that for dimensional consistency must be multiplied by the (inverse) square root of a characteristic timescale. Since our aim here is to describe the growth of DM perturbations using quantities related to the linear regime, we find it natural to choose the timescale |δ c /δ c | |Ḋ(t)/D(t)|, where the factor D(t) defined below Eq. (2) effectively describes the linear growth of perturbations under gravity in a given cosmological background. Moreover, the adopted noise is 'multiplicative', in that its strength depends on the state of the system, and specifically on the ratio σ/δ c . Regions with σ δ c tend to change their mass more abruptly, while the evolution is slower for σ δ c . In particular, positive variations of M (or in ln σ −1 ) within the filtered region can be reasonably related to mergers among collapsed halos, or mass accretion from the field; negative variations can be interpreted in terms of mass loss due to gravitational interactions with surrounding regions, tidal forces, stripping, and fragmentation. We stress the crucial role played by the multiplicative noise in Eq. (1); as η(t) fluctuates, also the random variable σ and hence the multiplicative factor σ/δ c on the r.h.s. varies, and therefore σ η/δ c is not null even if η is; this noise-induced drift actually makes σ(t) to copy the decrease of δ c (t) with time and, given the inverse convex shape of the deterministic function σ(M ), a net average increase in mass M (t) is enforced. Similar stochastic models with multiplicative noise have been employed to describe a wide range of physical phenomena, from Brownian motion in inhomogeneous media or in close approach to physical barriers, to thermal fluctuations in electronic circuits, to the evolution of stock prices, to the heterogeneous response of biological systems and randomness in gene expression; to our knowledge, this is the first time they are applied to describe the formation of collapsed structure in the Universe.
A more quantitative hint that such a stochastic equation has some value can be derived by performing a simple integration via the Euler-Maruyama method (e.g. Kloeden & Platen 1992; there are other ways to obtain more accurate numerical solutions of stochastic differential equations but these are not needed here); Eq. (1) can be discretized on a time grid t i with i = 0, . . . , n − 1 as follows where w i are random weights extracted from a normal distribution with zero mean and unit variance. In Fig. 1 we show the resulting evolution of ln σ(M (t)) as a function of cosmic time, with initial condition M ∼ 10 4 M or ln σ ∼ 2.2 at z ∼ 10 (reasonable initial conditions do not change significantly the outcome). The grey lines are 30 randomly chosen evolutionary tracks, while the red line and shaded area illustrate the median and the quartiles over 3000 realization of the noise, and the dashed blue line shows the evolution of the spherical collapse threshold δ c (t). It is easily seen that the noise induces a drift of ln σ, making it to decrease, and hence making the mass M to increase. Remarkably, after a burn-in period needed to erase memory of the initial condition, σ(M (t)) tends to copy the evolution of δ c (t); this means that the noise-induced drift effectively renders the standard hierarchical clustering of the halo population, as expressed by the increase with time of the characteristic mass M c (t) set by the condition σ(M c (t)) ∼ δ c (t).

Fokker-Planck equation and the Press & Schechter mass function
We now look for the probability density P(M, t) for a region to be in a state of mass between M and M + dM at time t; this can be found by solving the Fokker-Planck equation associated to Eq. (2), which reads (see Appendix A for details) where we have defined the two quantities D( factorizes the mass and time dependencies. The Fokker-Planck equation may also be written as a pure continuity equation in terms of the average comoving matter densityρ M . Now to solve the Fokker-Planck equation we employ the transformations: Then Eq. (5) turns into which is a standard diffusion equation. In terms of these new variables, the boundary conditions stated below Eq. (5) read lim X→∞ W = 0, W(X, 0) = δ D (X) and (∂ X W)| X=0 = 0. These stem from the following circumstances: (i) for reasonable power spectra σ is a slowly-varying inverse function of M , so that X ∝ σ −1 tends to zero or infinity as M does; (ii) the collapse threshold δ c ∝ D −1 (t) scales inversely with t, so that The solution of this differential problem is standard (it can be easily found via a Fourier transform) and writes Coming back to the original variables we get which after Eqs. (6) and (7) yields the Press & Schechter mass function: Three interesting remarks follow. First, note that the presence of the multiplicative noise in Eqs. (1) and (2) is fundamental in originating a second-order, diffusion-like term in the associated Fokker-Planck equation (see also Appendix A); from the derivation above it is seen that such a term yields the exponential cut-off of the mass function at the high-mass end. Second, one can easily compute the moments (12) For scale-free power spectra M ∝ σ −6/(n+3) holds in terms of the effective spectral index n > −3, so that the above implies M k (t) ∝ δ −6 k/(n+3) c ; this in turn scales as M k (t) ∝ t 4 k/(n+3) in the redshift range where δ c ∝ t −2/3 applies. Third, it is interesting to note that the Fokker-Planck Eq. (5) highlights that the evolution of the probability density function is driven by the two terms on the r.h.s: the first represents its noise-induced drift toward larger masses because of hierarchical collapses, and the second describes its diffusive reshaping at the high-mass end due to the stochasticity in merging/accretion events.

Stationarity
An alternative derivation, that will be useful for the generalization in the next Section, is the following. We again start from Eq. (1) and change variable from σ to ν ≡ δ c (t)/σ(M ): thus now the variable ν is seen to undergo a stochastic Ornstein-Uhlenbeck process. The corresponding Fokker-Planck equation reads We set the boundary conditions lim ν→∞ P = 0 and J | ν=0 = −|δ c /δ c | [ν P + ∂ ν P]| ν=0 = 0, implying ∞ 0 dν P = 1; the former is the natural boundary due to the expected exponential suppression of the mass function for large values of M , which correspond to large ν ∝ σ −1 ; the latter is the no-current boundary condition at the ν = 0 point enforced by the constraint P(ν, t) = 0 for ν < 0 being ν positively defined.
Having incorporated the time dependent quantity δ c (t) into the new variable ν, under the ergodic hypothesis one expects that the relevant solution in terms of this variable should be stationary (e.g., Paul & Baschnagel 2013), i.e., it should satisfy ∂ t P(ν, t) = 0. From Eq. (14), P(ν, t) =P(ν) is determined by where the constant on the r.h.s. must be zero to satisfy the no-current boundary condition J | ν=0 = 0. The solution to this simple ordinary differential equation with normalization ∞ 0 dνP(ν) = 1 isP = 2/π e −ν 2 /2 . The mass function writes as given that |dν/dM | = (δ c /σ 2 ) |dσ/dM |, this is indeed easily recognized to be again the Press & Schechter mass function.
One may wonder whether the general time-dependent solution P(ν, t|ν , t ) of Eq. (14) for a generic initial condition P(ν, t = t |ν , t ) = δ D (ν − ν ) converges to the stationary state, and over which timescale. To this purpose, we note that defining a new time variable τ ∝ − ln δ c (t) brings Eq. (14) into a form that can be easily solved via a Fourier transform; the fundamental solutions (a general result for Gaussian and Markovian variables known as Doob's theorem) read Taking into account the no-current boundary condition, and reverting to the original time variable one obtains P(ν, t|ν , t ) = 1 with δ ≡ δ c (t) and δ ≡ δ c (t ). Plainly, away from any initial condition, for δ substantially lower than δ , this converges to the stationary solution P(ν, t) derived above. Such transitional states could be possibly related to deviation of the mass function from the self-similar shape Eq. (16) in terms of the variable ν. Note, in passing, that the transition probability Eq. (17) cannot be directly related to the halo conditional mass function of the extended Press & Schechter theory; we anticipate that to derive the latter a modification of Eq. (13) is needed, but the issue demands an extended analysis that will be presented in a forthcoming paper.
3. MASS-DEPENDENT THRESHOLD: THE N -BODY MASS FUNCTION It is well known that the halo mass function derived from N -body simulations deviates substantially from the Press & Schechter shape (see reference in Sect. 1 for details); in particular, the former is flatter than the latter both at the high and at the low mass end, and evolves more slowly toward high redshift (see Fig. 2). This mismatch is usually cured by introducing a mass-dependent threshold for collapse: where q, β, and γ are parameters to be set by fitting the N -body outcomes. Such a modified threshold is generally ascribed, though a bit naively, to the fact that perturbations undergo an ellipsoidal rather than a spherical collapse (Sheth & Tormen 2002; see also discussion by Mo et al. 2010). Note, however, that the above shape is quite general in describing a variety of phenomena that can influence the collapse, like tidal torques and angular momentum, cosmological constant, dynamical friction (see Del Popolo 2017 and references therein). In the excursion set framework, the parameters q ≈ 0.707, β ≈ 0.47 and γ ≈ 0.615 are required to fit the N -body mass function; however, there is some degeneracy, in that for example a square-root barrier with q ≈ 0.55, β ≈ 0.5 and γ ≈ 0.5 fits simulations equally well. We stress that values of γ > 1/2 are of some concern in the excursion set framework, because they imply that some walks δ(σ) do not cross the barrier at all, an occurrence thought to represent fragmentation. Notice that no general analytical expression exists for the excursion set mass function, apart   (24)  from very particular barrier shapes (e.g., the constant, linear or square-root barriers; see for example Mahmood et al. 2006;Giocoli et al 2007;Lapi et al. 2013) and in general one must recur to numerical solutions (see Zhang & Hui 2006).
In the framework presented here, we aim to show that a minimal modification of the basic Eq. (1), which incorporates a mass-dependent collapse threshold with shape analogous to the above Eq. (18), will yield a mass function in excellent agreement with N -body simulations. Moreover, we will provide an analytic expression valid for any triples of parameters q, β, and γ < 1/2, that in the limit ν >> 1 (i.e., large masses and/or early times) matches the empirical formula introduced since Sheth & Tormen (1999); this will establish a direct connection between the parameters describing the barrier and the shape of the halo mass function (which is absent in the excursion set approach). We start from Eq. (1) by replacing, in the multiplicative noise term, the standard spherical collapse threshold δ c (t) with the mass-dependent δ c (σ, t) given above; thus now the ratio σ/δ c (σ, t) modulates the noise toward enforcing collapse. We also retain the term |δ c /δ c | 1/2 |Ḋ(t)/D(t)| 1/2 involving the characteristic timescale for the linear growth of perturbations (see discussion in Sect. 2). When formulated in terms of ν, Eq. (1) modified in such a way writes as where δ c (t) is the standard threshold for spherical collapse. The corresponding Fokker-Planck equation reads: with boundary conditions lim ν→∞ P = 0 and J | ν=0 = −|δ c /δ c | [ν P + (ν/B) ∂ ν (ν P/B)]| ν=0 = 0, implying the constraint ∞ 0 dν P = 1. In analogy to the procedure followed in Sect. 2.2, we look for stationary solutions P(ν, t) =P(ν) with ∂ tP = 0; one obtains the equation where the integration constant must be null to satisfy the no-current boundary condition J | ν=0 = 0. The above equation can be easily solved by multiplying both sides by 1/B and recognizing that it becomes separable for the function (ν/B)P; we find where the normalization constant A is determined by the condition ∞ 0 dνP(ν) = 1. We stress that this result holds for any mass-dependent collapse threshold that can be expressed as δ c (σ, t) ≡ δ c (t) B(ν)/ν in terms of the scaled variable ν.
Specializing now to the shape of B(ν) from Eq. (18), we note that for ν → 0 the behavior νP(ν) ∝ B(ν) ∝ ν 1−2γ applies; to satisfy the normalization constraint ∞ 0 dνP = 1 one must require γ < 1/2. Performing explicitly the integration, we get the closed form expression Fig. 3.-Evolution with redshift of the quantity (2 Y ) −1/2 entering the colored noise solution of Sect. 4, that in the white-noise case is simply the spherical collapse threshold δc (dotted line). Colored curves illustrate the result for different values of the volatility ζ and mean reversal rate ω characterizing the noise: blue line is for ω = 10 Gyr −1 , green for 1 Gyr −1 , and red for 0.1 Gyr −1 ; solid lines refer to ζ = ω, while dot-dashed line to ζ = 2 ω, and dashed line to ζ = ω/2 (these last two instances shown only for ω = 1 Gyr −1 ).

The resulting mass function just writes
incidentally, note that the multiplicity function f (ln σ) used in some literature works is just f (ln σ) ≡ νP(ν). We stress that the above is an exact expression, valid for any triple of values q, β, γ < 1/2; the Press & Schechter function is recovered for q = 1 and β = 0. Remarkably, the asymptotic behavior for ν >> 1, corresponding to large masses and/or early cosmic times, is seen to produce a shape akin to the empirical fit of N -body simulations adopted since Sheth & Tormen (1999); however, for finite ν the terms in the exponential are important and must be taken into account. We now use a Levenberg-Marquardt least-squares minimization routine to fit the multiplicity function νP(ν) to the simulation outcomes for FoF halos by Sheth & Tormen (1999), Bhattacharya et al. (2011) and Watson et al. (2013), sampled from their (somewhat different) fitting formulas for −0.6 ≤ ln[σ(M ) D(z)] −1 ≤ 1.3 in mass bins of 0.25 dex. For Sheth & Tormen (1999), we find best fit parameters q ≈ 0.62, β ≈ 0.16 and γ ≈ 0.37, yielding a value of the normalization constant A ≈ 0.63. For Bhattacharya et al. (2011), we obtain q ≈ 0.69, β ≈ 0.09 and γ ≈ 0.42, yielding A ≈ 0.64. For Watson et al. (2013), we get q ≈ 0.69, β ≈ 0.12 and γ ≈ 0.37, yielding A ≈ 0.66. These triples of values, reported for convenience in Table 1, are consistent within the uncertainties in the simulation results and in the fitting procedure.
In Fig. 2 we compare the mass function N (M, t) from Eqs. (23) and (24) to the N -body results by Bhattacharya et al. (2011), finding an excellent agreement over a wide range of masses M ∼ 10 8 − 10 16 M and redshifts z ∼ 0 − 10. In the same Figure we also plot for reference the Press & Schechter mass function; moreover, we show the mass function computed from the excursion set approach (following the numerical algorithm by Zhang & Hui 2006) with the barrier shape of Eq. (20) and the standardly adopted parameters q ≈ 0.707, β ≈ 0.47 and γ ≈ 0.615 (see above).

COLORED NOISE AND NON-MARKOVIAN
WALKS In Nature, white noise is never found to be perfectly realized, but rather constitutes an idealization of the stochastic driving force affecting a physical phenomenon. Thus one may wonder whether the previous treatment can be extended to a colored instead of a white noise; this will correspondingly enforce a non-Markovian evolution of the system. To have a grasp on the impact of colored noise and non-Markovianity on the mass function, we consider a multiplicative stochastic Ornstein-Uhlenbeck process; for simplicity we adopt a mass-independent collapse threshold, so that the endpoint of this computation should be compared with the Press & Schechter mass function. Specifically, we modify our Eq. (2) into the two-dimensional stochastic system where, besides already defined quantities, Γ(t) is an Ornstein-Uhlenbeck noise with average Γ(t) = 0 and a nontrivial correlation between different times controlled by the parameters ω and ζ (both have physical dimension of 1/time); ζ represents the degree of volatility, i.e. the sensitivity of the system to random changes, while ω is the dissipation rate at which the system tends to reverse toward its zero mean. This is perhaps the simplest generalization of the white noise case, since Γ(t) is still Markovian, but of course M (t) is not. In the limit ζ = ω → ∞ one recovers a pure white-noise, since Γ(t) Γ(t ) lim ω→∞ ω e −ω |t−t | = 2 δ D (t − t ). The Fokker-Planck equation regulating the dynamics of the probability density function P(M, Γ, t) for the above system is in analogy with the one-dimensional case, we can write this is as a continuity equation ∂ t P + ∇ · J = 0 in terms of the vectorial differential operator ∇ in terms of: This is the equivalent for colored noise of Eq. (10), which is recovered in the limit ζ = ω → ∞ as Y (t) → dt T 2 (t) = 1/2 δ 2 c . The corresponding expression for the mass function is written as: this is similar to the Press & Schechter shape, but for a modified redshift evolution encoded in Y (t). The quantity (2 Y ) −1/2 , that in the white-noise limit is just δ c , can be regarded as a modified collapse threshold; in Fig. 3 we show how its evolution and absolute value differ from δ c (t), depending on the volatility ζ and mean-reversal rate ω characterizing the noise. We conclude that only values of ζ ∼ ω several Gyr −1 are required not to move far away from the Press & Schechter mass function, and hence from simulations.
5. SUMMARY AND OUTLOOK In this paper we have submitted a new theory of the hierarchical clustering based on stochastic differential equations in real space, that constitutes a change of perspective with respect to the excursion set formalism; this work is specifically focused on the halo mass function.
First, we have invented a stochastic differential equation that describes fluctuations in the mass growth of DM halos, as driven by a multiplicative white (Gaussian) noise dependent on the spherical collapse threshold and on the power spectrum of DM perturbations. By numerically integrating such a stochastic differential equation, in Sect. 2 we have demonstrated that the noise yields an average drift of the halo population toward larger masses, that quantitatively renders the standard hierarchical clustering (see Fig. 1). Then, in Sect. 2.1 we have solved the Fokker-Planck equation associated to the stochastic dynamics, and obtained as a solution the Press & Schechter mass function; in Sect. 2.2 we have pointed out that the solution is stationary when the original equation is written in convenient variables.
Then in Sect. 3 we have introduced a minimal modification of the stochastic equation in terms of a massdependent collapse threshold. Using a parametric shape analogous to that adopted in the excursion set framework, we have obtained a closed-form analytical solution of the associated Fokker-Planck equation. Remarkably, such a solution has a limiting shape for large masses/early times similar to the empirical fitting formula introduced since Sheth & Tormen (1999); in fact, for specific values of the parameters describing the mass dependence of the collapse threshold, our result reproduces extremely well the N -body mass function over a wide range of masses and redshifts (see Fig. 2).
As an aside issue, in Sect. 4 we have generalized our stochastic approach to a colored, instead of a white, noise; in particular, we have investigated the modification to the Press & Schechter mass function when the stochastic dynamics is ruled by a multiplicative Ornstein-Uhlenbeck noise with finite volatility and mean-reversal rate. We have exactly solved the related Fokker-Planck equation, finding that the mass function has shape analogous to the Press & Schechter one when expressed in terms of a modified, effective collapse threshold; the latter may substantially differ from the standard δ c (t) in absolute value and time evolution, depending on the correlation parameters of the noise (see Fig. 3). We conclude that values of such parameters larger than several Gyr −1 are required not to move far away from the Press & Schechter mass function, and hence from simulations.
The next-order development of this work will concern the computation of the conditional mass function, i.e., the mass function of a halo's progenitors. This investigation will naturally extend to merger rates, formation time distributions, and large-scale halo bias. A more detailed comparison of our results with the outcomes of N -body simulations, that includes the specificity of both the numerical experiments as well as of the theory, will be welcome. Other future applications could involve a re-examination of the two-phase mass growth of DM halos, the halo specific angular momentum distribution, the void mass function, and halo statistics in non-standard cosmological frameworks. We very much hope that the new perspective offered by the theory presented here will contribute to a better understanding of the gravitational dynamics leading to the formation and evolution of DM halos and hosted baryonic structures across cosmic times.
We thank our referee for a constructive report, and for the insightful comments and helpful suggestions. We acknowledge Carlo Baccigalupi, Alessandro Bressan, and Giovanni Bussi for enlightening discussions and critical reading. This work has been partially supported by PRIN MIUR 2017 prot. 20173ML3WW 002, 'Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes'. A.L. has taken advantage of the MIUR grant 'Finanziamento annuale individuale attivitá base di ricerca' and of the EU H2020-MSCA-ITN-2019 Project 860744 'BiD4BEST: Big Data applications for Black hole Evolution STudies'.

A. A PRIMER ON THE STOCHASTIC DIFFERENTIAL AND FOKKER-PLANCK EQUATIONS
Given that concepts and techniques related to the stochastic differential and Fokker-Planck equations are not very common among the astrophysics community, for the reader's convenience we present here a short primer, in a modern notation and systematic way. In particular, we focus on the derivation of the Fokker-Planck equation associated to a given stochastic system, in presence of a state-dependent, multiplicative noise; this is extensively used in the main text. More details and applications can be found, e.g., in the book by Risken (1996).
The derivation involves two steps: (i) an expression for the time derivative of the probability density in terms of a Taylor-series of the conditional moments, known as Kramers-Moyal expansion; (ii) the explicit computation of such moments for a random variable satisfying a stochastic differential equation with multiplicative noise. Suppose we are given a system characterized by a physical variable ξ, whose evolution ξ(t) as a function of time t is stochastic. If the evolution is Markovian, by definition the probability density function P(x, t) of finding the system in state ξ(t) = x at time t satisfies: in terms of the conditional (transition) probability P(x, t + τ |x , t) between the times t and t + τ ; in other words, for a Markovian system the transition probability depends only on the value at the next earlier time. We rewrite the integrand as P(x, t + τ |x , t) P(x , t) = P(x + ∆ − ∆, t + τ |x − ∆, t) P(x − ∆, t) in terms of ∆ ≡ x − x and then perform a Taylor expansion in ∆ to obtain Now we insert this expression in Eq. (A1) and perform the integration after changing variable from x to ∆; noting that in the n = 0 term d∆ P(x + ∆, t + τ |x, t) = 1 holds since the conditional probability is normalized, we get where we have defined the conditional moments M n (x, t; τ ) ≡ |ξ(t + τ ) − ξ(t)| n | ξ(t)=x = d∆ ∆ n P(x + ∆, t + τ |x, t) .
Now we Taylor-expand the moments with respect to τ as follows M n (x, t; τ )/n! D n (x, t) τ + ϑ(τ 2 ) ; (A5) note that terms of order τ 0 cannot be present since P (x + ∆, t|x, t) = δ D (∆) by definition and in Eq. (A4) all the conditional moments for n ≥ 1 vanish. For future reference the coefficients D n are defined as All in all, we obtain the so called Kramers-Moyal expansion in terms of the partial differential equation this ends the first step in the derivation. We now compute explicitly the coefficients D n when the variable ξ(t) satisfies a stochastic differential equatioṅ with inital condition ξ(t) = x. Here η(t) is a white (Gaussian) noise with ensemble-average properties η(t) = 0 and η(t) η(t ) = 2 δ D (t − t ); the coefficient 2 in this last expression is only a convenient arbitrary choice, since it can be reabsorbed into the multiplicative function g without loss of generality. We start by transforming the differential into an integral stochastic equation We expand near Then we iterate for ξ(t ) − x in the integrand to get Now taking the ensemble average and using the properties of the white noise yields where in the last term we have used that t t dt 2 δ D (t − t ) g(x, t ) = g(x, t ) since the Dirac-δ operates on an extremal of the integration. Dividing by τ and taking the limit for τ → 0 one immediately recognizes the coefficient D 1 = h(x, t) + g(x, t) ∂ x g(x, t). For higher-order coefficients notice that terms containing the noise are proportional to τ n where n is the number of integrals involved, and vanish for small τ ; actually only one of such terms, containing two integrals and two noises contributes and yields D 2 = lim τ →0 (1/2τ ) t+τ t dt t t dt g(x, t ) g(x, t ) 2 δ D (t − t ) = g 2 (x, t), while D n = 0 for any n ≥ 3. This ends the second step of the derivation.
Putting together the coefficients just derived in the Kramers-Moyal expansion of Eq. (A7), one finds the Fokker-Planck equation corresponding to the original stochastic equation: The quantity g ∂ x g appearing in the coefficient D 1 is a noise-induced drift; this stems from the fact that as η(t) fluctuates, also the random variable ξ(t) and so the function g(ξ(t), t) varies and therefore g(ξ(t), t) η(t) is not null even if η(t) = 0 is. Finally, simple algebra shows that the Fokker-Planck equation may be written as a source-free continuity equation: ∂ t P(x, t) + ∂ x J (x, t) = 0 (A14) in terms of a probability current J (x, t) ≡ h(x, t) P(x, t) − g(x, t) ∂ x [g(x, t) P(x, t)] .
for the marginalized P(M, t) = dΓ P(M, Γ, t) with boundary conditions lim M →∞ P = 0, P(M, 0) = δ D (M ) and ( dΓ Γ P)| M =0 = 0. As a first step, we introduce a new variable X ≡ dM/D(M ) in place of M , and redefine the probability density as W(X, Γ, t) = D(M ) P(M, Γ, t); then the above equation turns into We now perform a two-dimensional Fourier transform W(X, Γ, t) ∝ dk X dk ΓW (k X , k Γ , t) e i (k X X+kΓ Γ) , and obtain the following equation for the Fourier modes Given the boundary conditions, it is convenient to look for solutions with shapẽ where Σ XX (t), Σ XΣ (t), Σ ΣΣ (t) are only functions of time. Inserting this ansantz into the previous equation yields the following ordinary differential equations These can be straightforwardly solved as Inverting the Fourier transform in Eq. (B3) one finds out the solution where ||Σ|| = Σ XX Σ ΓΓ − Σ 2 XΓ is the determinant of the 2 × 2 symmetric matrix constructed with the Σs; this can be easily checked to satisfy the desired boundary conditions, that in terms of the variable X and function W read lim X→∞ W = 0, W(X, 0) = δ D (X), and ( dΓ Γ W)| X=0 = 0.
Finally, marginalizing over Γ and coming back to the original variables, one obtains in terms of the quantities: note that the correct white-noise limit is recovered for ζ = ω → ∞ since Y (t) → dt T 2 (t) = 1/2 δ 2 c (t), and that actually (2 Y ) −1/2 constitute an effective collapse threshold, dependent on the parameters of the colored noise.