Culminating the Peak Cusp to Descry the Dark Side of Halos

The ConflUent System of Peak trajectories (CUSP) is a rigorous formalism in the framework of the peak theory that allows one to derive from first principles and no free parameters the typical halo properties from the statistics of peaks in the filtered Gaussian random field of density perturbations. The predicted halo mass function, spherically averaged density, velocity dispersion, velocity anisotropy, ellipticity, prolateness, and potential profiles, as well as the abundance and number density profiles of accreted and stripped subhalos and diffuse dark matter, accurately recover the results of cosmological N-body simulations. CUSP is thus a powerful tool for the calculation, in any desired hierarchical cosmology with Gaussian perturbations, of halo properties beyond the mass, redshift, and radial ranges covered by simulations. More importantly, CUSP unravels the origin of the characteristic features of those properties. In this paper, we culminate its construction. We show that all halo properties but those related to subhalo stripping are independent of the assembly history of those objects, and that the Gaussian is the only smoothing window able to find the finite collapsing patches while properly accounting for the entropy increase produced in major mergers.


INTRODUCTION
One fundamental limitation for the theory of structure formation is the lack of an exact analytic treatment in non-linear regime. Indeed, such a treatment is available for monolithic spherical collapse. But the typical seeds of dark matter (DM) halos are ellipsoidal, and so is also their collapse. Moreover, halos suffer major mergers, so their collapse is rather lumpy. This is the reason that halo properties have been traditionally studied by means of numerical simulations (see e.g. the review by Frenk & White 2012).
But simulations are not without problems. They are very CPU-time consuming, and hence, only deal with moderately large numbers of particles. As a consequence, their dynamic range and spatial resolution are limited. In addition, their analysis involves complex selection procedures similar to those used in observations which may bias the results. Last but not least, even though they are very useful to determine halo properties, they are not well-suited to explain the origin of those properties.
To shed light on that 'dark side' of halos, structure formation has also been studied by analytic means. Press & Schechter (1974) used the top-hat spherical collapse e.salvador@ub.edu approximation to infer the halo mass function (MF) from the statistics of perturbations in the initial Gaussian random density field. This approach was refined by correcting for cloud-in-cloud configurations (Bond et al. 1991), extended to calculate conditional MFs and merger rates (Bower 1991;Lacey & Cole 1993), and modified to account for ellipsoidal collapse (Monaco 1995;Sheth et al. 2001;Sheth & Tormen 2002). Some authors (Colafrancesco et al. 1989;Peacock & Heavens 1990;Appel & Jones 1990;Bond & Myers 1996;Paranjape & Sheth 2012) went a step further and took into account that halos evolve from maxima (peaks) in the filtered initial density field (Doroshkevich 1970;Bardeen et al. 1986, hereafter BBKS). One particularly elaborate treatment along this line was the so-called ConflUent System of Peak trajectories (CUSP) formalism developed by Manrique & Salvador-Solé (1995, which assumed the existence of a halo-peak correspondence, calibrated by means of the results of simulations. Other authors concentrated in deriving the density profile for halos that result from the spherical collapse of the homogeneous mass distribution around a density perturbation, the so-called secondary-infall model (Gunn & Gott 1972;Fillmore & Goldreich 1984;Bertschinger 1985), or from the spherical collapse of the perturbation itself (Avila-Reese et al. 1998;Manrique et al. 2003;Salvador-Solé et al. 2007). This approach was also pursued within the peak theory (Del Popolo et al. 2000;Ascasibar et al. 2004;MacMillan et al. 2006;Salvador-Solé et al. 2012a,b).
However, the analytic approach faces the following apparently insurmountable fundamental difficulties (Ds): • D i) Even though halos seem to arise from the collapse of patches with negative total energy corresponding to peaks in the smothed initial density field (Hahn & Paranjape 2014), the relation between halos and protohalos in simulations is not the expected one (Ludlow & Porciani 2011).
• D ii) Moreover, there seems to be no one-to-one correspondence between halos and peaks. Peaks of a given density contrast are overcrowded (Appel & Jones 1990), and halo seeds often split in a few disjoint nodes (Porciani et al. 2002).
• D iii) The edge of a virialized halo 1 is a fuzzy concept, there being many possible halo mass definitions (e.g. Davis et al. 1985).
• D iv) There is no clear argument in favor of any particular smoothing window (e.g. Bond 1988) and the mass encompassed by a window of a given scale is unknown. 2 • D v) Peaks are triaxial (Doroshkevich 1970) and the time of their ellipsoidal collapse, along one axis first (pancakes), then another axis (filaments) and the third axis in the end (Lin et al. 1965;Zel'Dovich 1970) is unkown.
• D vi) As DM is collisionless, monolithic collapse is followed by shell-crossing (Hénon 1964), with no analytic treatment not even assuming spherical collapse.
• D vii) In addition, the real collapse of halos is lumpy, i.e. they suffer major mergers producing a violent relaxation with no analytic treatment and a poorly known final state (Lynden-Bell 1967;Shu 1978).
• D viii) The material assembled in halos is a mixture of diffuse DM (dDM) and other halos (Angulo & White 2010), which become subhalos tidally stripped by the host potential well. The situation is thus very complex (e.g. Ghigna et al. 1998).
• D ix) In addition, massive subhalos suffer dynamical friction and eventually merge at the halo center, whereas the exact treatment of dynamical fric-tion is only available for infinite homogeneous systems (Chandrasekhar 1943).
Yet, using CUSP, Salvador-Solé et al. (2012a,b) were able to derive not only the typical density profile, but also the velocity dispersion and anisotropy profiles and even the ellipticity and prolateness profiles of halos in CDM cosmologies (see also Viñas et al. 2012 for WDM cosmologies). Juan et al. (2014a,b) formally proved the basic hypothesis of CUSP that there is a one-to-one correspondence between halos and peaks and rederived the halo MF and density profile from first principles and with no single free parameter. Lastly, Salvador-Solé, Manrique, & Botella (2021a,b,c) have recently derived the properties of halo substructure.
All CUSP predictions are in very good agreement with the results of simulations, so this formalism can be used to extend the halo properties beyond the radius, mass, and redshift ranges and cosmologies covered by simulations. More importantly, CUSP unravels the origin of all those properties and their characteristic features. Yet CUSP has not attracted much attention possibly because of its non-linear construction over many years, the lengthy mathematical developments included in its development, and the fact that it has never been clearly shown how CUSP solves or avoids the above mentioned Ds for an accurate analytic treatment of structure formation so that the prejudice that this is not possible remains.
More importantly, CUSP has two weak points. First, the derivation of all halo properties is made assuming monolithic collapse, while real halos suffer major mergers. The fact that CUSP recovers, nonetheless, the results of simulations seems to imply that the inner properties of halos do not depend on their assembly history perhaps because in major mergers halos lose the memory of their past history (Salvador-Solé et al. 2012a). Some results of simulations support, indeed, this idea (Huss et al. 1999;Hansen et al. 2006;Wang & White 2009;Barber et al. 2012), but others point to the opposite conclusion (Gottlöber et al. 2001(Gottlöber et al. , 2002Sheth & Tormen 2004;Fakhouri & Ma 2009Hahn et al. 2009). Second, the properties of peaks from which all halo properties follow depend on the particular smoothing window used to filter the initial density field. CUSP employs the Gaussian window for practical reasons (it greatly simplifies the calculations; BBKS), but the reason why this particular window leads to such good results is poorly understood.
In this paper we address these issues. We prove that the properties of halos do not depend on their assembly history (except for the properties related to subhalo stripping), and show that the use of a Gaussian filter is mandatory for the filtering of the initial density field to properly trace DM clustering. With these results we culminate the construction of CUSP.
Taking advantage of the opportunity, we provide a compact orderly review of CUSP skipping all detailed mathematical developments and showing, instead, how CUSP solves or circumvents the above mentioned Ds and clarifies some intriguing questions (Qs) risen by simulations, listed in the end of the Paper. To this end, each time one D is solved or one Q is answered it is indicated in parentheses.
The layout of the Paper is as follows. In Section 2 we establish the halo-peak correspondence at the base of CUSP. In Section 3 we derive the halo and subhalo MFs from peak counts. The inner properties of halo seeds are derived in Section 4 and those of halos following from them are given in Section 5. The role of major mergers and the Gaussian window is addressed in Section 6. The implications of our results are discussed in Section 7. Throughout the Paper we provide some Figures in order to illustrate the goodness of the CUSP predictions, calculated for the LCDM WMAP7 cosmology and using the BBKS power spectrum with Sugiyama (1995) shape parameter.

HALOS AND PEAKS
In hierarchical cosmologies with Gaussian density perturbations, DM clustering is fully determined by the linear power spectrum. Thus, by filtering the density field at any arbitrary small time t i with varying graining scales so as to uncover the collapsing patches of different masses, it should be possible to reconstruct the growth, abundance, and inner properties of halos at any time t > t i . But is there really any smoothing window able to do that?
Given the isotropy of the Universe, the window must be spherical. On the other hand, as collapsing patches are finite, the window must be of compact support. Lastly, even though we do not understand why (see Sec. 6 for the explanation), simulations with finite resolution converge, so the window must also be of compact support in Fourier space. The Gaussian is the only window satisfying these necessary conditions, and as shown next it works, indeed.

Ellipsoidal vs. Spherical Collapse
In top-hat spherical collapse, all peaks with a fixed positive density contrast δ th at any scale R th in the density field at t i collapse at the same time t, and give rise to halos with mass M satisfying (e.g. Peebles 1980) In equations (1) and (2) ρ c (t) is the cosmic mean density at t, δ th c (t) and D(t) are the critical linearly extrapolated density contrast for spherical collapse at t, and the linear growth factor (equal to 1.686 and to the scale factor a(t), respectively, in the Einstein-de Sitter cosmology). The relation (2) is equivalent to the following one in terms of the 0th-order top-hat spectral moment σ th 0 , which ensures that the seed of one halo at different t i is the same evolving perturbation with height ν th (M, t) ≡ δ th (t, t i )/σ th 0 (M, t i ) = δ th c (t)/σ th 0 (M, t). However, in Gaussian density fields, peaks are triaxial (for whatever filter), so the patches they encompass undergo ellipsoidal collapse, which invalidates the previous relations. Indeed, the time of collapse (along all three axes) t c depends in this case not only on the mass and size, but also on the triaxial shape and central concentration of the patch, i.e. on the density contrast δ, scale R, ellipticity e, prolateness p, and curvature x (defined in eq. [13]) of the associated peak: t c = t c (R, δ, e, p, x). It is thus not unsurprising that peaks with identical δ th (t) but different R th (M ) collapse at different times, and that the masses of spherical patches with R th (M ) differ from the masses M of halos resulting from their collapse (D i).
Nonetheless, as the e and p probability distributions functions (PDFs) of Gaussian-filtered peaks with δ and x at scale R are very sharply peaked at their maximum values like the x PDF itself (BBKS), we can take them with fixed values at the respective maxima, e max , p max , x max . Thus, provided that there is some relation between the mass M of halos at the cosmic time t (or the corresponding patches at t i ) and the scale R of their associated peaks as provided by any given halo mass definition (see below), all ellipsoidal patches with different masses M identified by triaxial peaks with a positive δ at the scale R(M, t) will collapse essentially at the same cosmic time t = t i + t c [R(M, t), δ, e max , p max , x max ], with t a monotonous decreasing function of δ as in top-hat spherical collapse (D v). Certainly, the small scatter ∆ in the e, p, and x PDFs around the respective maxima will translate into a scatter in the time of collapse, which will propagate into the MF of halos at t, as well as on the radius r of the ancestor at t of the halo with final mass M , implying (in halos growing insideout; see below), a scatter in the value at r of any halo profile, ∆ξ(r) = dξ dr To calculate those scatters we need the partial derivatives of t c with respect to e, p, and x), which can be estimated through numerical experiments. But the purpose of the present Paper is not to derive such scatters, but the mean halo properties, and this can be done analytically.

Halo-Peak Correspondence and Halo Mass Definition
Thus, for any given halo mass definition, Gaussian ellipsoidal collapse leads to a one-to-one correspondence between halos with M at t and peaks in the density field at t i filtered on scale R similar to that found in top-hat spherical collapse. Moreover, to guarantee the arbitrariness in t i the functions δ(t) and R(M, t) defining that correspondence must depend on t i just as in top-hat spherical collapse, or, using the 0th-order Gaussian spectral moment instead of the scale R, Hereafter σ j stands for the jth Gaussian spectral moment. Note that the halo-peak correspondence does not depend on the small scale mass distribution within each patch determining whether the collapse is monolithic or lumpy (i.e. Q ii). Specifically, in all cosmologies, the consistency conditions that all the DM in the Universe must be locked inside halos of different masses and that the mass M of a halo must be equal to the volume-integral of its density profile allow one to find (Juan et al. 2014a) the functions r δ (t) and r R (M, t) or r σ (M, t) setting the correspondence associated with every specific halo mass definition (D iii). 3 In all cases of interest, these functions are well-fitted by the analytic expressions, with S(t) defined as  Komatsu et al. (2011) with the values of coefficients d, s 0 , s 1 , s 2 and A for some relevant cases, quoted in Table 1. (See App. A for an approximate expression of r R (M, t).) Strictly speaking, the previous halo-peak correspondence does not take into account that, due to cloudin-cloud configurations, some nested peaks do not accomplish their collapse at t because, before that, they are captured by the more massive halo corresponding to the host peak and become subhalos. Thus, the previous halo-peak correspondence must be corrected for failed halos and the corresponding nested peaks in order to have a strict one-to-one correspondence between real virialized halos and non-nested peak. In turn, similar one-to-one correspondences are foreseen between subhalos and nested peaks at different levels (see below).

Halo Growth and Peak Trajectories
In hierarchical cosmologies, low-mass halos are much more abundant than high-mass ones, so minor mergers are extremely frequent, and give rise to a smooth mass increase of the massive partner called accretion. Accretion does not alter the quasi-equilibrium state of the accreting object; there is just a gentle virialization of the accreted matter through shell-crossing. In turn, accreted halos survive as subhalos. In contrast, major mergers are rare events causing the destruction of all partners (usually two) and the virialization ex novo of the whole system through violent relaxation. This relaxation takes a few crossing times, so a substantial fraction of all halos at any time t (usually the most massive ones, with longer crossing times and formed more recently; e.g. Raig et al. 2001) are not fully virialized. It is thus unsurprising that the seeds at t i of ∼ 15 − 20% of them still have multiple nodes (Porciani et al. 2002) (D ii). Thus, when calculating the halo MF at t, we will include all collapsed halos at that time, regardless of whether they are fully virialized or not, but when calculating their typical (steady) properties we will consider only virialized objects.
Taking into account the relation ∂δ(r, R) ∂R = R∇ 2 δ(r, R) ≡ −x(r, R)σ 2 (R)R (13) satisfied by a Gaussian window, the Taylor expansion of δ(r, R) at scale R + ∆R (∆R R) of a peak implies there can only be one peak within a distance ∆R from another at scale R (Manrique & Salvador-Solé 1995). We can thus readily identify the two peaks tracing any individual accreting halo. All peaks with varying R tracing the same accreting halo describe a continuous δ(R) trajectory in the δ-R plane. Note that the peaks along any of those trajectories are not anchored to a fixed point; their position r slightly sloshes around. However, the series of (non-concentric) patches of different masses they encompass as δ diminishes are embedded within each other because the separation of contiguous peaks at R and R + ∆R is at most ∆R.
The continuous peak trajectory δ(R) traced by an accreting halo satisfies the differential equation Similarly, the peak trajectory traced by a halo accreting at the mean instantaneous rate satisfies the same equation but with the curvature x(δ, R) replaced by the inverse of the mean inverse curvature of peaks with δ at R (Manrique & Salvador-Solé 1995) or, given the very sharply peaked x-distribution, simply the mean peak curvaturex(R, δ), calculated in BBKS. For moderately high peaks as it corresponds to halos of galactic scales,x(R, δ) is well-approximated by γν, where γ is σ 2 1 /(σ 0 σ 2 ) (BBKS), and equation (14) leads to Approximating the power spectrum by a power-law, P (k) ≈ Ck n , we have d ln δ/d ln R ≈ m and δ(R) ∝ R m , with m = −[(n+3)/2] 3/2 , showing that the form of those peak trajectories is nearly universal in all realistic cosmologies. This result is what causes (App. A) all scaled halo profiles to be nearly universal (Q i). Another interesting consequence is that, in the limit of vanishing R, peak trajectories converge to a finite value. Note that, contrarily to what happens in the excursion set (ES) formalism (Bond et al. 1991), where the δ(R) trajectories describe random walks, the δ(R) trajectories in CUSP are continous monotonously decreasing functions of R (eq. [14]). This is well-understood: ES follows the variation in the mass of patches around fixed points when the scale R of the smoothing k-sharp window increases, whereas CUSP follows the variation in the mass of patches around moving peaks when the scale R of the smoothing Gaussian window increases. In the former case the density contrast at the point may increase or decrease, while in the latter their curvature of the peak automatically determines how much δ decreases. (The curvature of the peak is indeed connected with the accretion rate of the corresponding halo; see eq. [14].) Far from being a problem, the monotonous decrease of δ(R) with increasing R better reflects the monotonous increase of halo masses with increasing cosmic time. This difference between the two formalisms is what allows CUSP to monitor the growth of individual halos, while ES can only monitor the growth of halos in a statistical way (see also Sec. 6).
But, when an accreting halo suffers a major merger, its associated peak trajectory is interrupted. No peak at R + ∆R can be identified to the peak at R (in fact the two merging peaks become saddle points at the new scale; Cadioul et al. 2020) and a 'new' peak with the same δ on a scale R substantially larger than those of the progenitors appears (it cannot be identified to any peak at the previous scale R − ∆R), which traces the newly formed halo. The rates at which peak trajectories disappear or appear due to major mergers can be used to calculate the halo destruction and formation rates (via these events) as well as their formation and destruction time PDFs (Manrique et al. 1998). The fact that, due to major mergers, the number density of (non-nested) peak trajectories in the δ-R plane progressively declines with increasing R is precisely what motivated this formalism to be dubbed ConflUent System of Peak trajectories.

(SUB)HALO MASS FUNCTIONS AND PEAK ABUNDANCES
When a halo is accreted by another halo or merges with it, it becomes a subhalo or is destroyed, respectively. But its own first-level subhalos are never destroyed. Either they become second-level subhalos in the former case or they are passed as first-level subhalos in the new halo in the latter case. And so on for higher-level subhalos.
Given the one-to-one halo-peak correspondence, this halo and subhalo behaviour automatically leads to a parallel behaviour of the corresponding peaks. The result is a complex nesting network of peaks with any fixed density contrast at different scales (Appel & Jones 1990). But this network is not a problem: it simply tells us that to count halos at any t we must count non-nested peaks with the corresponding δ and to count subhalos of any desired level we must count nested peaks of the same level (D ii).

Halo MF and Abundance of Non-Nested Peaks
The number density of peaks with δ at scales between R and R + dR is the number density of peaks at scale R with density contrastδ greater than δ that cross such a density contrast when the scale is increased to R + dR or, equivalently, with δ satisfying the condition Therefore, it is equal to the integral of the density of peaks with heightν =δ/σ 0 (R) and curvaturex per infinitesimalν and dx, calculated by BBKS, overδ in the range delimited by the inequality (16) and overx in the range of all possible values. The result is (Manrique & Salvador-Solé 1995) (17) where ν = δ/σ 0 (R) is the peak height, and R is defined as √ 3 σ 1 /σ 2 . As mentioned, the mean peak curvature, x(R, δ), is separable in a first approximation, so is also N pk (R, δ).
But the number density (17) refers to all peaks with δ between R and R+dR, while virialized halos correspond to non-nested peaks only. The number density of nonnested peaks with δ between R and R + dR, N (R, δ)dR is the solution of the Volterra integral equation where N nst pk (R, δ|R , δ) is the conditional density of peaks with δ at R subject to being nested in a peak with δ at R given below (eqs. [21]-[25]). The comoving density (18) of peaks with δ at R at t i is thus equal to the comoving density of halos with M at t or MF, N (M, t), but for the change of variable from R to M (eq. [8]).
The new change of variable from M to σ 0 leads to the halo multiplicity function at t, defined as And, using the variable (see Fig. 1), with the Gaussian height ν approximately equal to [a(t)/D(t)] times the top-hat one, ν th (App. A). As shown in Juan et al. (2014b), the FoF mass with linking length b = 0.2, FoF(0.2), is equivalent to the SO mass with overdensity ∆ vir relative to the mean cosmic density, SO(∆ vir ). This explains that the halo multiplicity function for FoF(0.2) is privileged (Q iii) in the sense that it has a roughly universal shape in top-hat filtering (Juan et al. 2014b).

Subhalo MF and Abundance of Nested Peaks
Following the same procedure above but from the conditional density of peaks with infinitesimalν andx at R subject to lying at a distance r (in units of the tophat filtering radius associated with R) of a peak with ν at R , per infinitesimalν andx calculated by BBKS, we can also compute the conditional number density of peaks with δ at scales between R and R + dR subject to being at a distance r (same units) from such a peak (21) wherex(R, δ, R , δ , r) is the mean curvature of peaks with δ and R at a distance r from another peak with δ at R , and ν (r) stands for δ (r)/σ 0 (R )g(r, R ), where we have used the notation: e(r) = 1 − 2 (r), where δ (r) and ∆δ (r) are the mean and rms density contrast, respectively, at r from a peak (BBKS).
Thus, the conditional number density of peaks with δ at scales R to R + dR subject to being nested in a peak with δ at R is (Manrique et al. 1998) where factor with S equal to the mean non-nested peak separation (same units) drawn from their mean density (eq. [18]), is to correct for the overcounting of host peaks as some of them are also nested. But, if we are interested in the number of first-level subhalos in a halo (from now on, numbers are denoted by a calligraphic N so as to distinguish them from number densities), the conditional number density (22) is not enough. We need the conditional number of peaks with δ per infinitesimal scale around R subject to being nested it at first-level in a non-nested peak with δ at R . The result is where N nst (R, δ|R , δ) is the conditional number density of peaks with the same characteristics subject to being nested in non-nested peaks with identical δ at R (eq. [22]) corrected for nesting at any intermediate scale R between R and R , given by the solution of the Volterra equation Note that, in the integral on the right of equation (25), we have used the conditional number density of peaks with δ per infinitesimal scale around R subject to being directly nested within peaks at R (i.e. without being nested in any smaller scale peak), N d nst pk (R , δ|R , δ), calculated in Salvador-Solé, Manrique, & Botella (2021a) to avoid overcorrection for intermediate nesting (peaks can be nested in more than one intermediate scale peak).
We remark thatx(R, δ, R , δ, r), likex(R, δ), is very nearly separable, in the subhalo mass range, in a function of R and another function of the remaining arguments (Juan et al. 2014b). This separability then propagates (see Salvador-Solé, Manrique, & Botella 2021b for details) to the conditional number of nested peaks, N (R, δ|R , δ) (eq. [24]), which has the important consequences mentioned below.
Given the one-to-one correspondence between firstlevel subhalos and first-level nested peaks, the MF or number of accreted subhalos per infinitesimal mass in halos with M h at t h , N acc (M s ), coincides, but for the change of variable from scale R to subhalo mass M s , with the number of peaks per infinitesimal scale around R(M s ) and δ(t h ) that are directly nested into a peak with the same density contrast at the larger scale R(M h ). The corresponding cumulative subhalo MF N acc (> M s ) (see Fig. 2) is essentially proportional to M −1 s (Q x), and very nearly universal when M s is scaled to M h (Q i). The derivation of the MF of subhalos tidally stripped by the halo potential well is postponed to Section 5.3.2.

PROTOHALO PROPERTIES
The typical halo properties, hereafter denoted by subscript h, arise from those of the corresponding collapsing patches or protohalos, hereafter denoted by subscript p. (Quantities with no subscript refer to both objects indistinctly.) Specifically, since in linear regime the velocity and density fields are tightly related, protohalos are fully characterized by their spherically averaged density, ellipticity, and prolateness profiles, that is we need no kinematic information. To obtain those structural protohalo profiles we should deconvolve the height, ellipticity, and prolateness of the corresponding peaks. But this is not possible because the Gaussian window used to find those peaks yields the loss of information at large wavenumbers. Fortunately, as shown next, for purely accreting halos we can use the information on the peaks along the full δ(R) trajectory to achieve the desired deconvolutions. (The peak trajectories of ordinary halos suffering major mergers are interrupted at those events, so we have no access to that information.) For this reason we concentrate, until Section 6, in halos evolving by pure accretion (PA), the only ones for which we can determine the properties of the corresponding protohalos. Then, monitoring the collapse and virialization of those protohalos, we can derive the properties of the final purely accreting halos. This procedure can be followed for individual objects (provided their δ(R) trajectory is known) or for ensembles of halos with the same mass and redshift so as to derive their typical properties. In the present study we concentrate in this latter application of the method. Specifically, we derive the mean spherical averaged radial profiles and global quantities of (proto)halos with M h at t h .

Density Profile
Taking the origin of the coordinate system on the peak at scale R, the density contrast δ at r p = 0 is nothing but the convolution with the Gaussian window of that radius of the intrinsic density contrast field δ p (r p ) in the protohalo. We thus have after integrating over the polar angles, where δ p (r p ) is the spherical average of δ(r). Equation (26) shows that we can infer the peak trajectory δ(R) traced by a purely accreting halo from the density profile of its protohalo. Conversely, given the peak trajectory δ(R) of a purely accreting halo, we can readily solve the Fredholm integral equation of first kind (26) for δ p (r p ) (Salvador-Solé et al. 2012a). This can be done for individual halos or for halos with M h at t h accreting at the mean instantaneous rate, leading to the mean spherically averaged density profile δ p (r) of the corresponding protohalos. Note that, as in the limit of vanishing R, the peak trajectories converge to a finite value (null asymptotic slope), so do the unconvolved mean protohalo density profiles in the limit of vanishing r p .

Eccentricity Profiles
Similarly, reorienting the Cartesian axes j along the main axes of the triaxial peak with δ at R, the squared semiaxes of the peak, A 2 j (δ, R), equal to the second order spatial derivatives ∂ 2 j of the centered density contrast field δ p (r p ) scaled to the Laplacian (A 2 1 + A 2 2 + A 2 3 = 1) is given by (Salvador-Solé et al. 2012b) Writing δ p (r p ) in terms of δ p (r p ) and the protohalo axis profiles a p j(r p ) (App. C), and integrating over the polar angles, we arrive at where we have taken into account that the axes A j of peaks are are well-known functions of their curvature (BBKS), and introduced the function G p (r p ) defined in Appendix C. Equation (28) shows that the A j (R) trajectories of peaks can be obtained from the semiaxes profiles of their protohalos. Conversely, given a peak trajectory δ(R) determining the A j (R) trajectories, we can solve Fredholm integral equation of the first kind (28) for the quantities a 2 p1 (r p )/a 2 p j(r p ), in the same way as equation (26). Then, from the quantities a 2 p1 (r p )/a 2 p j(r p ), we can infer the ellipticity and prolateness profiles, or equivalently, the primary and secondary eccentricity profiles, p (r p ) and ε p (r p ), defined in Appendix C. We can thus invert equations (28) from δ(R) solution of equation (15) (with the partial derivative ∂δ/∂R equal to the total derivative) fixing the A 2 j (dδ/dR, R) trajectories, and determine the profiles a 2 p1 (r p )/[a 2 p j(r p )G p (r p )] for the three orientations j, leading to the primary and secondary eccentricity profiles, p (r p ) and ε p (r p ), for individual protohalos or for halos accreting at the mean instantaneous rate.

HALO PROPERTIES
The properties of purely accreting halos follow from those of the corresponding protohalos thanks to the conservation in ellipsoidal collapse and virialization of a few quantities specified below. In fact, as the latter properties follow from the peak trajectories traced by such halos, the properties of purely accreting halos ultimately follow from the characteristics of those peak trajectories.
The link between the properties of protohalos and halos is provided by the conservation of three quantities during ellipsoidal collapse and virialization: the mass inside the corresponding radii, the total energy inside those radii, and the volume delimited by them. The latter two quantities are not conserved in an absolute way (the system expands and contracts, and loses energy through shell-crossing during virialization), but relative to those of the spherically symmetrized system.

Density Profile
The one-to-one correspondence between halos and peaks guarantees the conservation of the mass in the triaxial protohalo at t i to that of the triaxial halo at t, M h = M p . As shown in Appendix A, this equality and the halo mass definition implicit in that correspondence fully determine the (roughly universal) mean spherically averaged halo density profile in purely accreting halos. However, to understand what is behind that formal derivation it is convenient to monitor in detail the ellipsoidal collapse and virialization of the system. For the reasons explained in Appendix B, this can be done assuming spherical symmetry.
In linear regime, homeoids expand radially by the factor D(t) (hence, homothetically) with neither shellcrossing nor energy exchange between different regions. After reaching turnaround, shells collapse and bounce, which causes them to cross other shells. When a shell moving inwards crosses another one moving outwards some (potential) energy is transferred from the latter to the former due to the change in the inner mass they see. And when they again cross in the opposite direction the energy exchange is reversed. But the amount of energy exchanged depends on the radius of the crossing: near pericenter, they exchange more energy than near apocenter. This difference causes a net energy flux outwards, i.e. from shells collapsed earlier, with smaller apocenters, to shells collapsed later, with larger apocenters. This energy loss of shells causes their apocenter loci to shrink so shells cross increasingly nearer the center, and exchange less energy. On the other hand, the correlation between their phases is slowly lost. Thus, the Lagrangian energy outflow diminishes. In the end, the energy exchange ceases, and shell apocenters stop contracting.
Note that shell-crossing proceeds with no apocenter crossing. Indeed, if particles in two different accreting shells did coincide at apocenter (hence with null radial velocity), they will coincide along the whole trajectory, so they would not belong to different accreting shells. Therefore, the particle apocenter loci contract, due to the energy loss through shell-crossing, in an orderly manner until they stop. In other words, accreting halos grow from the inside out (D vi). This particular growth allows one to relate the radius r h in the final object in equilibrium to the protohalo properties within the corresponding r p .
The virial radius of the accreting halo increases with increasing time as the object grows inside-out. We can thus virtually move one after the other the shells reaching turnaround without any crossing so as to match their apocenters in the final halo. 4 The energy profile of the resulting 'toy object' inside any radius r h , E h (r h ) will differ, of course, from that of the real halo E h (r h ) because the latter accounts for the energy loss through shell-crossing, while the toy object is built with no shell-crossing. In fact, E h (r h ) equals the total energy of the system at turnaround or directly at the initial time, E h (r h ) = E p (r p ). The toy object is of course not in equilibrium, but this can be fixed. As E h (r h )− W h (r h ), where W h (r h ) is the potential energy of a homogeneous sphere with M h (r h ) at r h , must be positive, 5 we can expand virtually every inner shell, avoiding shell-crossing, so as to end up with a uniform density equal to the mean density of the real halo inside r h , and still have an excess of spherical kinetic energy inside that radius. This kinetic energy can then be redistributed over the sphere, exchanging the radial and tangential components of the spherical velocity variance so as to satisfy the spherical virial relation with null spherical radial velocity variance. We are then led to a steady homogeneous toy object with the same mass M h (r h ) and radius r h as the real halo, though with potential energy W h (M h ) = −3/5 GM 2 h /r h (M h ), total energy E h (M h ) = E p (M p ), and null surface term. It thus satisfies the virial relation where , which is given, in the parametric form, by with ρ p (r p ) = ρ c (t i )[1 + δ p (r p )] calculated in Section 4. H i is the Hubble constant at t i , and is the peculiar velocity caused by the mass excess within r p , where we have neglected the velocity dispersion of DM particles at t i , and taken the cosmic virial factor f (Ω) ≈ Ω 0.1 at t i equal to one. Note that equation (29) explains why the virial radius of halos is recovered in spherical top-hat collapse assuming energy conservation (Q v), as if the system had not lost energy through shell crossing. This result thus justifies the definition of the halo virial radius Bryan & Norman 1998. The inversion of r h (M h ) given by Equation (29) leads to the mass profile M h (r h ), and by differentiating it at the spherically averaged density profile ρ h (r h ) (see Fig. 3). Conversely, given a purely accreting halo with mass profile M h (r h ), we can calculate the total energy of its seed, E p (M p = M h ) (eq. [29]), and determine ρ p (r p ) from (eqs. [30] and [31]). Thus there is in PA a one-toone correspondence between the density profile for halos and protohalos. These previous results hold for the profiles of individual halo-protohalo pairs as well as for mean profiles.
In addition, the inner asymptotic protohalo density profile, . Equation (29) then leads to ρ h (r h ) ∝ r α h , and, since α is null (see Sec. 4.1), we conclude that there is strictly no cusp in the mean spherically averaged density profile for CDM halos (Q iv). Nevertheless, the finite central value is approached slightly more slowly than in the Einasto profile.

Eccentricity and Kinematic Profiles
The kinematics and triaxial shape of halos are determined by the triaxial shape of protohalos. Since this relation is not included in the one-to-one halo-peak correspondence, there is no shortcut in this case; the only way to derive it is by taking into account the above mentioned conservation relations during ellipsoidal collapse and virialization. Moreover, for obvious reasons, we cannot assume spherical symmetry in this case. However, we will take into account that, to leading order in the the departure from spherical symmetry, hereafter simply the 'asphericity', all physical properties of evolving triaxial systems coincide with those of their spherically symmetrized counterparts (see App. B).
To leading order in the asphericity, the volume of ellipsoids coincides at any time with that of the corresponding spheres in the spherically symmetrized system. That is, the ratio between the two volumes must be conserved, to leading order in the asphericity, over the ellipsoidal collapse and virialization of the system. We thus have where a 1 , a 2 , and a 3 are the semiaxes of the ellipsoids and r is their equivalent radius (see the definition in App. C), equal to the radius of the corresponding sphere in the spherically symmetrized system. Similarly, to leading order in the asphericity, the total energy within radius r of the triaxial system, E(r), coincides at any time with that E(r) in the spherically symmetrized system. Note that the energy lost by a sphere of initial radius r p through shell-crossing during virialization is indeed accounted for in E(r). What is not accounted for is the energy exchange between the sphere and the rest of the system in the triaxial system, so that energy loss is included in the residual δE(r). Consequently, the ratio between E(r) and E(r) must be conserved to leading order in the asphericity. This conservation has two consequences. First, as E(r) = E(r)+δE(r), it implies with δE(r), given in Appendix B, dependent on the eccentricities of the halo. Second, as E(r) is the same in the spherically symmetrized system, the transfer from the radial to the tangential kinetic energy due to the non-radial motion produced in the non-linear evolution of the triaxial system must go in parallel to a departure of the potential energy so as not to alter the total spherical energy. In other words, the fractional velocity variance transferred from the radial to the tangential direction (equal to half the fractional 1-D tangential velocity variance generated) must be equal to half the typical fractional deviation of the potential, Equation (35) relates the velocity anisotropy of halos with their eccentricity (Q vi). Those two conservation relations determine the triaxial shape and kinematics of the halo (App. C). Indeed, the halo eccentricity profiles, h (r) and ε h (r), arise from those of protohalos, p (r) and ε p (r) (eqs. [C23] and [C25]) through a relation (eq. [C26]) that involves the halo velocity variance σ 2 h (r). This relation is enough to determine the halo eccentricities near the center, where they coincide with those of protohalos with well-known sharply peaked PDFs (BBKS). But to solve the problem at any r h the closure relation (35) is required.
Writing the anisotropy β h (r h ) as a function of the scaled potential variance (eq. [35]), the latter as a function of the halo eccentricities (App. C), and the eccentricities as functions of the halo velocity variance and the protohalo eccentricities (eqs. [C24] and [C25]), the anisotropic equilibrium equation to leading order in the leads to a differential equation for σ 2 h (r h ). Its solution for the usual boundary condition of null velocity at infinity leads to a pseudo phase-space density profile ρ h (r h )/σ 3 h (r h ) of the same power-law form as in Bertschinger's (1985) spherical self-similar collapse model. This characteristic pseudo phase-space density profile found in simulations (see Fig. 4) is thus due to the (quasi-)self-similar collapse of a (quasi-) homogeneous mass distribution, and the phase-mixing produced by shell-crossing like in Bertschinger (1985) model (Q vii). The only difference, apart from the non-strict-similarity of the real collapse, is that Bertschinger assumes spherical symmetry, while CUSP takes into account the ellipsoidal collapse leading to a partial conversion of the radial velocities in linear regime to tangential velocities (eq. [35]).
Once σ h (r h ) is known, we can compute the halo scaled density-potential covariance profile from that in the protohalo (eq. [C24]), and then the scaled density and potential variance profiles (App. C). The latter leads, through equation (35), to the velocity anisotropy profile β h (r h ), which adopts the universal form (Hansen et al. 2006) found in simulations (Q viii).
These kinematic profiles lead to the eccentricity profiles of the halo from those of the protohalo (eqs [C24]-[C26]). In peaks with small and moderately large average curvature, the typical e p and p p values peak at 1.7 and 1.3, respectively, regardless of the values of δ and R (BBKS), implying that p and ε p are about 0.81 and 0.64, respectively. Since U (r h ) is close to one (see Fig. 2 in Salvador-Solé et al. 2012b), the predicted mean h (r h ) and ε h (r h ) profiles at small and moderately large radii are constant and close to 0.9 and 0.8, respectively, consistently with the results of simulations (Q ix). At large radii, the rms density fluctuation profile tends to a power-law of index ∼ −0.1, the eccentricities are outwards decreasing, and the isopotential contours become more spherical than the isodensity contours (Salvador-Solé et al. 2012b) as found in simulations (Q ix).
Once again, not only are the kinematic and eccentricity profiles of halos determined by the eccentricity profiles of the protohalo, but the latter can be obtained from the former (eqs. [C24]-[C25]). And all the previous results hold for the profiles of individual halo-protohalo pairs as well as for their mean profiles.

Accreted Subhalos
According to the results of Section 3.2 and the insideout growth of halos in PA, the abundance of (first-level) subhalos per infinitesimal non-truncated mass around M s inside the radius r h in halos with current mass M h and virial radius R h at t h , N acc (< r h , M s ), coincides, after the appropriate change of variables, with the mean number per infinitesimal scale of peaks with δ[t(r h )] at R(M s ) nested (at first-level) in the peaks with scale R{δ[t(r h )]} found along the mean peak trajectory of those accreting halos, with mean (sperically averaged) mass profile equal to M h (r h ). We thus have, Form now on a bar on a function of r means its mean value inside r. The dDM mass fraction at r h in the halo with M h at t h directly accreted from the inter-halo medium is essentially equal to the dDM mass fraction f acc dDM (t) in that medium (t > t min ) at the time t = t(r h ), being N (M, t) the halo MF at t (Sec. 3). Indeed, halo clustering starts at some finite time t min with a minimum halo mass M min equal to the free-streaming mass associated with the DM particle in the real Universe or the resolution mass in cosmological simulations. Consequently, all the DM that at t min should lie in halos with masses M < M min remains in the form of diffuse DM (dDM) that is progressively accreted onto halos (D viii) (see Fig. 6). As can be seen from equation (37), the above mentioned separability of N [R s , δ|R(δ), δ] in a function of R s and another function of the remaining arguments propagates (with an extra factor of R s arising from the δ-derivative and the change of variable from R s to M s ) to the separability of the same kind of N acc (r h , M s ).
The mean spherically averaged number density profile for accreted subhalos with M s , n acc (r h , M s ), defined as N acc (r h , M s )/(4πr 2 h ), scaled to the corresponding total mean number density,n acc (R h , M s ) = 3N acc (M s )/4πR 3 h , is a function of r only. Taking into account the relation (38), it can be written as Equation (39) shows that the mean spherically averaged scaled number density profiles of accreted subalos for all masses M s overlap in one only profile, which follows the scaled density profile of the halo. The ultimate reason for this result is thus the above mentioned separability of the conditional number of nested peaks (Q xi).

Stripped Subhalos
Soon after being accreted, subhalos begin to be tidally stripped (or even disrupted) by the host potential well so that they liberate other stripped subhalos and dDM previously locked inside them, and end up with smaller truncated masses, M tr s . Note that, as halos undergoing PA grow inside-out, once subhalos are accreted their orbits are kept unaltered, and since their initial conditions at accretion are independent of mass, the orbits are also. Certainly, subhalos suffer two-body interactions between themselves and with dDM. But these interactions do not affect the very numerous low-mass subhalos because they are already accounted for in shellcrossing (see the discussion in Salvador-Solé, Manrique, & Botella 2021a). Only rare massive subhalos suffer an extra effect, the so-called dynamical friction. However, the typical properties of substructure are defined from large subhalo ensembles, so rare massive objects are excluded from the present treatment (D ix).
The abundance N stp (r h , M tr s ) of stripped subhalos per infinitesimal radius and truncated subhalo mass around r h and M tr s is then given by where N stp [M,t(r h )] (r, M tr s ) is the differential abundance of stripped subsubhalos with M tr s at r in a subhalo with mass M when the host halo had mass M (r h ). In equation (40), N acc (v, r h , M ) is the abundance of accreted subhalos per infinitesimal mass, radius, and tangential velocity around M s , r h , and v, respectively, which factorizes in N acc (r h , M ) (eq. [39]) times the massindependent tangential velocity distribution of the halo at r h , and R(r h , M ) and R tr (v, r h , M ) are the original and truncated radii of subhalos with original mass M and tangential velocity v at r h , respectively. Double angular brackets denote average over the tangential velocity v of subhalos at their apocentric radius r h .
The first term on the right of equation (40) gives the abundance of subhalos accreted at r h with original nontruncated mass M s that give rise by direct stripping to subhalos with truncated mass M tr s , and the second term gives the number of subhalos with truncated mass M tr s that are released at r h from the stripping of more massive subhalos.
Indeed, after subhalos are accreted (their apocenters are stabilized) at r h , they begin to be repeatedly stripped and heated at each passage by the pericenter at r per (v, r h ), causing them to reach a new equilibrium configuration at next apocenter. Taking into account the truncation undergone in every passage under the impulse approximation (González-Casado et al. 1994), the subhalo ends up, after p passages by pericenter (from t(r h ) to t h ), with a truncated-to-original mass ratio given by (Salvador-Solé, Manrique (Q xii), where the functions m(r h ) and q(r h ) can be calculated from the mass and radius those subhalos would have had they not suffered any stripping during virialization and Q i (v, r h ) is the ratio of subhalo radii after and before truncation at the i-th pericentric passage. The ratio Q i (v, r h ) is the solution of the equation where Q(v, r h ) stands for r per (v, r h )/r h and f (c) for ln(1 + c) − c/(1 + c) in the case of that the density profile of halos (or subhalos) is adjusted by the NFW (Navarro et al. 1996) analytic profile (see Juan et al. 2014b for the Einasto profile). In equation (42), c(r h ) = r h c(R h )/R h is the concentrations of the halo and c i (v, r h ) is the concentration of the subhalo when it reaches equilibrium at apocenter after the i-th pericentric passage is the solution of the recursive equation where h(c) stands for f (c)(1 + c)/{c 3/2 [3/2 − s 2 (c)] 1/2 }, being s 2 (c) the isotropic 3D velocity variance of objects with mass M and radius R scaled to cf (c)GM/R, and u(v, r h ) is the heated-to-original subhalo energy ratio going together with tidal truncation, assumed of the form For K = 0.77 and β = −0.5 very good agreement is found with the results of numerical experiments (Salvador-Solé, Manrique, & Botella 2021b). These two parameters are the only ones appearing in the whole theory. They are necessary because subhalo stripping and heating do not arise from the collapse and virialization of peaks, and hence, must be modeled apart. According to the previous relations the truncated-tooriginal subhalo mass ratio (eq. [41]) turns out to be independent of the original mass of subhalos as a consequence of their similar concentration when they are accreted and the way they are truncated a pericenter (D xi) (Salvador-Solé, Manrique, & Botella (2021b)).
After some algebra, the solution of the integral equation (40) takes the compact form (Salvador-Solé, Manrique, & Botella 2021b) where µ(r h ) is the v-averaged mass-independent truncated-to-original subhalo mass ratio profile and f rel (r h ) is the proportion of stripped subhalos with M tr s at r h previously locked as subsubhalos which have been released in the intra-halo medium when their subhalo hosts have been stripped. The quantity f rel , which can be found by solving a Fredholm equation given in Salvador-Solé, Manrique, & Botella (2021b), never exceeds 6 %, so it can be neglected in a first approximation. Equation (45) implies that the mean spherically averaged scaled number density profiles for stripped subhalos with M tr s , n stp (r h , M tr s ) defined as usual in terms of the corresponding abundance N stp (r h , M tr s ) takes the form with and f acc dDM (r h ) are the total dDM mass fraction and the dDM mass fraction of accreted subhalos, respectively, at the radius r h of the halo with M h at t h (D viii). The fraction f stp dDM (r h ) is given by and a similar relation holds for f acc dDM (r h ) in terms of N acc (r h , M s ).
Equation (46) shows that the mean spherically averaged scaled number density of stripped subhalos of different truncated masses M tr s overlap in one only profile as a consequence of the separability of the abundance of accreted subhalos through equation (45). But this profile is bent with respect to ρ h (r) because the dDM mass fraction increases towards the halo center due to the more marked stripping there (see Fig. 5).
Taking into account the relation (48) that follows from equations (47) and (45) and the separability of N stp (r h , M tr s ) and N acc (r h , M tr s ), the integration over r h of equation (45) leads to the relation between the respective MFs. On the other hand, from the f stp dDM (r h ) and f acc dDM (r h ) dDM mass fractions profiles (Fig. 6), we can compute the global total and accreted dDM mass fractions in halos. We find that about ∼ 90 % of the mass of current Milky Way-mass halos is in the form of dDM (Q xiii).
Notice that all halo profiles can be expressed in terms of ρ h (r), f acc dDM (r), µ(r) and f rel (r). Since ρ h is roughly universal as a function of the scaled radius, r/R h , and the f acc dDM profile following from f acc dDM (t) and t(r/R h ) is also universal (see App. A), µ and f rel are too (Salvador-Solé, Manrique, & Botella 2021b). Thus, all properties regarding accreted and stripped subhalos in purely accreting halos (or having suffered the last major merger a very long time ago) are roughly universal (Q i).

MAJOR MERGERS AND GAUSSIAN WINDOW
The imprints of tidal stripping on subhalos and dDM are never erased, so these components retain the memory of the halo past history. Thus, the abundance and radial distribution of stripped subhalos and the associated dDM in halos depend on the halo assembly history (see Salvador-Solé, Manrique, & Botella 2021c for their derivation from those of purely accreting halos).
But all the remaining halo properties, including those regarding accreted subhalos and dDM, are independent of the halo assembly history, and can be obtained assuming PA (D vii). This surprising result, shown next, explains the empirical fact that all relaxed halos have similar properties regardless of the time of their last major merger (see also the discussion in Sec. 7).
In PA all halo properties evolve in a continuous manner. When a homeoid (or ellipsoidal shell) collapses, it crosses the homeoids having previously collapsed, with similar triaxial shapes and orientations. The energy Figure 6. Total dDM mass fraction as a function of radius (upper ticks) and cosmic time when the halo reached that radius (lower ticks) predicted by CUSP for the same halos as in Figure 5 in a 100 GeV WIMP universe (solid red line) and an Aquarius-like (Springel et al. 2008)  outflow (in Lagrangian coordinates) set in such shellcrossings causes the particle apocenter loci to contract (and change their triaxial shape) until the correlation between the orbital phases of particles in neighboring shells is lost. Virialization is thus achieved through the randomization of orbital phases that goes together with the contraction of the system. However, the order of the particle apocenter loci is kept unaltered, implying that the radial mapping of the system is essentially conserved. It is thus unsurprising that the radial profiles and triaxial shape and kinematics of halos can be inferred from those of protohalos, and vice-versa. Consequently, there is no increase of entropy, meaning that the randomization of particle velocities through shellcrossing, which affects their orbital phases, is not a real relaxation of the system.
On the contrary, there is no continuity of halo properties interrupted in major mergers. Particles previously assembled in each of the progenitors and orbiting within it suddenly falls towards the new center of mass of the system where they mix up. This mixing yields a new phase of less symmetric shell-crossing around the center of mass and energy pumping outwards until the parti-cle orbital phases are randomized, and the system stops contracting like in accretion. However, in this case, the particle apocenter order within each progenitor is not translated into the new structure, so there is a randomization of both particle phases and apocenters themselves, which causes the scrambling of particles in the new system. Such a scrambling destroys and smears out the progenitors so that the system becomes triaxial at all scales, from the largest one corresponding to the peak tracing the new halo, downwards. In other words, the memory of the initial multipolar mass distribution at small scales is lost, and the resulting mass distribution cannot be inverted anymore. Consequently, there is an increase of entropy, meaning that the randomization of particles velocities, which affects in this case the order of their apocentric radii and not just their orbital phases, is a real (violent) relaxation of the system. Shortly, the chaotic orbits followed by particles during the virialization of halos in major mergers cause the entropy increase in the two processes (Beraldo et al. 2019).
Of course, the structural profiles (mass density, number density of accreted subhalos, and the accreted dDM density), kinematic profiles, and eccentricity profiles of merged halos can be inverted as if they had been set by PA. The result will be the structural profiles (mass density and number density of nested peaks) and eccentricity profiles of fake protohalos evolving into the merged halos by PA. The mean profiles of real merged halos are thus equal to the mean profiles arising from their respective fake protohalos. But the latter mean profiles coincide with the mean profiles of real halos evolving by PA from protohalos because the curvature of a peak does not depend on the small-scale mass distribution inside (and outside) the protohalo, so the mean curvaturex of fake and real peaks evolving by PA with given density contrast and scale are the same. Thus the differential equation (14) governing the respective mean peak trajectories is also the same, and, since the two kinds of peaks have identical characteristics at the scale R and time t of the final halos, the two mean δ(R) trajectories coincide, just as the respective mean (R) and ε(R) trajectories following from them. Therefore, even though the small-scale mass distributions in protohalos of merged and purely accreting halos with given mass, triaxial shape, and kinematics at t are very different, their respective mean unconvolved profiles coincide.
Moreover, not only do the mean profiles of both kinds of halos coincide, but the profiles of every individual purely accreting halo coincide with those of one merged halo, and vice-versa. Indeed, given a merged halo, we can reshuffle its particles in the protohalo, keeping the convolved density profile and triaxial shape of the peak unaltered. By doing this the real protohalo is converted into a fake protohalo evolving by PA into the merged halo. Conversely, given a purely accreting halo, we can concentrate its particles along the radial direction in a few small scale undulations, keeping the convolved density profile of the peak unaltered, then along their homeoids, keeping the triaxial shape of the peak unaltered, so as to build up a few massive clumps (i.e. peaks on smaller scales) within the large-scale peak without altering its convolved properties. This way the real protohalo of the purely accreting halo is converted into a protohalo of a merged halo with the same profiles as the original purely accreting halo.
Thus both halo populations are indistinguishable, globally (the whole ensemble) as well as individually (every single realization). The only difference between them is in the entropy accompanying their collapse and virialization: in halos evolving by PA there is no entropy increase because the radial mapping of the protohalo is conserved, whereas in major mergers the entropy increases due to the reshuffling of particles produced in the merger at t, or in the protohalo at t i causing the new protohalo to evolve by PA to the same final halo. But the entropy increase in major mergers can only be detected provided the initial configuration is known. This is what happens in simulations where it is seen, indeed, that major mergers cause an entropy increase with respect to the case of accretion compatible with no increase (Obreschkow et al. 2020).
Thus halo properties do not depend on their assembly history (Q xiv). This result has in turn the following important consequence for the smoothing window.
As the halo assembly process carries the loss of information on small scales of collapsing patches, the smoothing window that identifies those patches must filter out those small scales, i.e. it must be of compact support in Fourier space. Any other window with compact support in physical space as required to have finite collapsing patches, such as the top-hat one, has a Fourier transform with non-null values out to infinity. And any other window with compact support in Fourier space, such as the k-sharp window, does not encompass finite collapsing patches because its inverse Fourier transform has non-null values out to infinity. The way this problem is circumvented in the ES formalism is by using the ksharp window with varying 0th-order spectral moment to identify the patches of different masses which should collapse according to the top-hat window with identical 0th-order spectral moments. Unfortunately, there is no reason for the masses seen by the two different windows with identical 0th-order spectral moment to coincide.
Thus the Gaussian is the only smoothing window that accounts for the entropy increase in halo clustering (due to major mergers), and finds in a fully consistent manner the finite patches with known time of collapse and virialization (D iv). In this respect we remark that CUSP does provide the relation between the 0th-order spectral moments of the Gaussian and the top-hat windows that see collapsing patches of the same mass (eqs. [9] and [11]).

SUMMARY AND DISCUSSION
CUSP is a rigorous formalism for the analytic treatment, from first principles and with no free parameters (but for the modeling of subhalo tidal heating), of halo clustering in hierarchical cosmologies. Taking advantage of the fact that the properties of halos do not depend on their assembly history, it uses the one-to-one correspondence between halos at t and peaks with known properties in the smoothed density field at t i to infer the abundance and properties of the final objects by monitoring the monolithic ellipsoidal collapse and virialization (through shell-crossing) of their seeds.
The achievements of CUSP are multiple: it overcomes numerous fundamental difficulties (Ds) for the analytic follow up of DM clustering; it allows one to extend the results of simulations to any hierarchical cold or warm DM cosmology, any halo mass-definition, and any arbitrary halo radius, mass, and redshift; and it unravels the origin of all macroscopic halo properties and their characteristic features. In particular, it answers the following intriguing questions (Qs) risen by simulations: • Q i) Why are the scaled spherically averaged profiles of halos so nearly universal?
• Q ii) Why is the halo MF well predicted assuming monolithic collapse?
• Q iii) Why does the halo multiplicity function privilege the FoF(0.2) halo-finding algorithm?
• Q iv) What is the inner asymptotic slope of the halo density profile?
• Q v) Why is the virial radius of a halo equal to that of a homogeneous spherical system having conserved the energy during collapse and virialization?
• Q vi) What is the relation between the halo triaxial shape and the velocity anisotropy?
• Q vii) Why is the pseudo phase-space density profile a power-law of index −1.875?
• Q viii) Where does the universal anisotropydensity relation come from?
• Q ix) What are the typical halo inner and outer asymptotic ellipticity and prolateness?
• Q x) Why is the subhalo cumulative MF close to a power-law of index −1?
• Q xi) Why is the number density profile of accreted subhalos indepencdent of subhalo masses and nearly proportional to the halo density profile?
• Q xii) Why is the truncated-to-original subhalo mass ratio independent of the subhalo mass?
• Q xiii) How much diffuse DM is there in halos, and what is its spatial distribution?
• Q xiv) Why are the properties of halos with very different assembly histories so similar?
• Q xv) Where does the halo assembly bias come from?
Furthermore, CUSP has led to the following findings: there is a one-to-one halo (non-nested) peak connection, extensible to subhalos and nested peaks of any level; during accretion halos grow inside-out so that virialization preserves the radial mapping of the system; in major mergers halos are scrambled so that violent relaxation causes the memory loss of the system; the inner properties of halos (except for those related to stripping) do not allow one to distinguish whether and when halos have suffered major mergers; the total abundance of accreted and stripped subhalos coincides except for their normalization depending on the abundance of accreted and total dDM; the stripping-related properties of ordinary halos having suffered major mergers can be inferred from those derived assuming PA (Salvador-Solé, Manrique, & Botella 2021c); and the filtering of the initial density field allowing one to monitor halo clustering must be carried out by means of a Gaussian window for consistency with the entropy increase produced in major mergers.
As mentioned, CUSP shows that the properties of halos do not depend on their assembly history. What is then the origin of the conflict between the results of simulations supporting this result (Huss et al. 1999;Hansen et al. 2006;Wang & White 2009;Barber et al. 2012), and those indicating that halo properties depend on the frequency of major mergers, the so-called 'assembly bias' (Gottlöber et al. 2001(Gottlöber et al. , 2002Sheth & Tormen 2004;Fakhouri & Ma 2009Hahn et al. 2009;Goh et al. 2019;Chen et al. 2020;Hellwing et al. 2021;Ramakrishnan, Paranjape, & Sheth 2021;Wang et al. 2020). According to CUSP, what really causes halo properties to depend on the local density as found in the latter simulations is not that halos lying in different local densities suffer major mergers at a higher rate (Fakhouri & Ma 2009), but the different background densities of their seeds. Indeed, the local density of protohalos affects not only the frequency of major mergers, which as shown has no repercussion in the inner halo properties, but also and mainly their accretion rate (Fakhouri & Ma 2010), which sets their inner properties through the inside-out growth of purely accreting halos, and of halos suffering major mergers as well.
In fact, if instead of monitoring the mean δ(S) trajectory for unconstrained peaks we had monitored that trajectory for peaks lying in backgrounds of different density contrasts δ (e.g. Manrique & Salvador-Solé 1995;Juan et al. 2014b) we would have been led to the mean profiles of halos in the corresponding final local densities.
Clearly, these mean profiles would be somewhat different from those found for all halos regardless of their location. We thus see that the local density of halos affects their typical inner properties, even though major mergers leave no imprints on them (Q xv).
Even though CUSP solves the main issues regarding DM halos it can still be improved. We are not thinking so much in minor effects such as the mass loss by halos in encounters and its eventual reaccretion (van den Bosch 2002; Wang et al. 2011) or dynamical friction which alters the properties of substructure regarding massive subhalos, but in more fundamental effects such as the gravitational torque between neighbouring (proto)halos or the gravitational drag of baryons. The former is at the base of the small angular momentum and tidally-supported elongation of halos (Doroshkevich 1970;White 1984;Salvador-Solé & Solanes 1993), and the latter yields notable deviations in halo structure from those of pure DM objects derived here (Gnedin & Ostriker 1999;Gnedin & Zhao 2002;Governato et al. 2012;. Funding for this work was provided by the Spanish MINECO under projects CEX2019-000918-M of IC-CUB (Unidad de Excelencia 'María de Maeztu') and PID2019-109361GB-100 (co-funded with FEDER funds) and the Catalan DEC grant 2017SGR643. We thank Prof. Ravi K. Sheth for fruitful comments.