A Mass Dependent Density Profile from Dwarfs to Clusters

In this paper, we extend the work of Freundlich et al. 2020 who showed how to obtain a Dekel-Zhao density profile with mass dependent shape parameters in the case of galaxies. In the case of Freundlich et al. 2020, the baryonic dependence was obtained using the NIHAO set of simulations. In our case, we used simulations based on a model of ours. Following Freundlich et al. 2020, we obtained the dependence from baryon physics of the two shape parameters, obtaining in this way a mass dependent Dekel-Zhao profile describing the dark matter profiles from galaxies to clusters of galaxies. The extension to the Dekel-Zhao mass dependent profile to clusters of galaxies is the main result of the paper. In the paper, we show how the Dekel-Zhao mass dependent profile gives a good description of the density profiles of galaxies, already shown by Freundlich et al. 2020, but also to a set of clusters of galaxies.


Introduction
As shown by a series of observations and gravitational effects at the cosmological level [1], together with effects at astrophysical scales [2,3], we know that the Universe cannot only be constituted of baryonic matter but should contain a non-baryonic mass/energy component with the property of clustering and which is dubbed dark matter (DM). Similarly, the galaxies appear constituted by an amount of baryons that make them visible, and a large halo of DM in which the baryonic component is embedded. Many studies have been performed to describe the DM halo density profiles. Some have showed that, from dwarf galaxy haloes to galaxy clusters, the density profile can be described by the Navarro-Frenk-White (NFW) profile [4,5], characterized by a double power law at small masses (ρ ∝ r −1 ), and large masses (ρ ∝ r −3 ). Further studies found a different kind of profile whose slope decreases going toward the center of the halo, e.g., [6][7][8][9]. These kinds of profiles are well described by the so-called Einasto profile of which we will speak later. The cuspy behavior of the NFW in the central region of the halo disagrees with observations of low-surface-brightness (LSB), dwarf satellite galaxies, and even to a small extent with galaxy clusters, which show flatter cores, e.g., [10][11][12][13][14][15][16][17][18]. In other words, the smallest of the inner slope of the density profile value predicted by dissipationless N-body simulations is larger than the values obtained by observations [19][20][21][22][23][24], in SPH (Smooth Particle Hydrodynamics) simulations [25,26], or in semi-analytical models [27][28][29][30][31].
This issue presents problems for the ΛCDM model at small scales and is called the Cusp/Core problem [10,11] (other problems often mentioned are (a) the missing satellite problem, namely the discrepancy between the number of subhaloes that N-body simulations predict, e.g., [32] and observations; (b) and the Too-Big-To-Fail (TBTF) problem, characterized by too many and too dense simulated sub-haloes with respect to observations [33,34]).
As mentioned, apart from the case of dwarf galaxies and LSBs, the cusp/core problem is also present at the scales of clusters of galaxies, mainly in their central part (10 kpc [16,17]).
As found by kinematics and lensing constraints in galaxies located in the center of relaxed clusters (cD galaxies, Brightest Central Galaxy, hereafter BCG), while the total mass profile, namely the sum of DM and baryons, is in agreement with the NFW predictions, the cluster's DM profile is flatter than an NFW profile [16,17,35,36].
Several ideas have been proposed to solve that discrepancy: broadly, they are cosmological or astrophysical.
The main idea at the root of astrophysics solutions resides in the existence of a "heating mechanism" that produces an expansion of the DM component of a galaxy. This, in turn, produces a reduction of the inner density. In the "supernovae feedback flattening" mechanism, supernovae explosions and/or AGN (Active Galactic Nucleus) outflows produce potential fluctuations that can heat the DM, producing a flattening of the cusp [25,26,[52][53][54][55][56].
The NFW, Einasto, and related density profiles are independent of the mass scale. In [66], it was shown that the inner slope of the density profile depends on the halo mass. This result was later confirmed by several simulations both in the case of galaxies and clusters e.g., [67][68][69][70][71][72][73][74].
Some papers [67,71,72,75] showed that the inner slope in SPH simulations displays regions with values similar to or larger than the NFW slope, and others that are much flatter, reaching a minimum for masses down to 10 8.5 M , in the case of [67]. This minimum is due to the fact that mass outflows of explosive events overcome halo gravity, giving rise to an expansion of the halo. For larger masses, the slope steepens and reaches values larger than the NFW slope when the stellar mass exceeds 10 10 M . At those masses, baryon accumulation gives rise to adiabatic contraction. The latter effect can be counteracted in the presence of AGN feedback [74,76].
In order to fit the density profiles, several parameterizations have been used. The Einasto profile [6,[77][78][79][80] presents two free shape parameters. Leaving free in the fit these two parameters, one can get good fits to DM cusps. Moreover, analytic expressions for the mass, the surface density, the gravitational potential, and quantities important for lensing can be provided. Unfortunately, the Einasto profile does not give a good fit to the innermost part of the density profile for galaxies [81]. Another profile often used is the generalized NFW (gNFW) model that performs in a similar way to the Einasto profile. Other density profiles used have adopted the averaged form.
(1) as in Zhao [82], where x = r/r c , with r c , a characteristic radius and ρ c , a typical density Note that the non-averaged DZ density profile coincides with gNFW profiles for the choice g = 3, b = 2). This is also the functional form of the gNFW density profile. Considering this form for the actual density profile, if b = n and g = 3 + k/n, where n and k can be any natural numbers, and one can obtain analytically the gravitational potential, the mass, and the velocity dispersion. This family of density profiles is able to provide one which yields very good fits to DM simulations even when baryons are present and is able to follow the profile even when cusps or cores are present. This profile, as shown by Dekel et al. [81], corresponds to n = 2 and k = 1, i.e., b = 2 and g = 3.5 in Equation (1). It remains with two free parameters, namely, a and c, is usually referred to as the Dekel-Zhao (DZ) profile, and captures cores better than other profiles [83]. In order to improve the fit to the density profiles of galaxies and clusters, a new approach has been introduced. The parameters of each fitting profile are expressed in terms of stellar mass, M * , and virial mass, M vir , or their ratio.  [83]. Similarly to Di Cintio et al. [86], they used the NIHAO project to obtain the functional form of the two parameters of the DZ density profile. The validity of Freundlich et al. [83], Di Cintio et al. [86], and all the mass dependent profiles obtained in the literature is limited to galaxies because the simulations at their base do not take account of AGN feedback. Therefore, their mass dependent profile cannot be built up to the clusters of galaxy mass range. In this paper, our goal is to extend the DZ profile to clusters of galaxies.
The NIHAO project is based on SPH simulations, whose core formation base mechanism focuses on "supernovae feedback". We previously mentioned the alternative core formation mechanism based on the exchange of energy and angular momentum between infalling gas clumps and DM. We used it in several papers, and describe the model in the next section. Generating structures of different size with it, we seek in this paper to find a relation between the parameters of the DZ density profile, from galaxies to clusters. The paper is organized as follows: In Section 2. we describe the model used. In Section 3, we summarize the steps in which the model works. In Section 4, we describe the DZ profile in general and in Section 5 the DZ mass dependent profile. In Section 6, we apply the previous profile to clusters of galaxies. Section 7 is devoted to conclusions.
The scheme is implemented according to the following stages: (i) Expansion of the diffuse gas and DM proto-structure in the linear phase, reaching a maximum radius before re-collapse of DM, forming a potential well for baryons to fall; (ii) Formation of stars from baryons radiative clumping in the halo center; (iii) Four parallel processes then follow; (a) baryon AC increases the DM central cusp, e.g., for 10 9 M galaxies, at z 5, see [27], (b) baryon-DM dynamical friction (DF) collapses clumps to the galactic center, (c) the halo central density reduces [57,58] from the transfer to DM, and stars [54,70,111] of DF energy and angular momentum (AM), contrary to the effect of AC; (d) DF and AC balance, opening the possibility for cusps to heat up, some up to core formation, as in spirals and dwarf spheroidals, while others, like giant galaxies, retains their steeper profile and cusp because of their deeper potential wells; (iv) Tidal torques (ordered AM), and random AM join their similar effects to DF.
Finally, the disruption of the smallest gas clumps, due to their partial conversion to stars, and the supernovae explosions repeated gas expulsion decrease stellar density, resulting in an additional slight core enlargement; see [65].

Density Profile Generation
We model the emergence of the density profile in the framework of the spherical model. From an initial, linear Hubble expansion, the expanding density perturbations eventually reaches a maximum turn-around, marking the onset of its recollapse [112,113]. A particle-based Lagrangian approach is used to compute the final density profile, noting the particle's initial and turn-around radii x i and x m (x i ), its turnaround density ρ ta (x m ) and its collapse factor In this approach, one computes the turn-around radius from the density parameter Ω i and the DM and baryon shell's average overdensity δ i Baryons start entirely in gas form, denoted by the "universal baryon fraction" [114], set to 0.167 in [115] over the total mass and set to f b = 0.17 ± 0.01, before stars form as discussed below (Sections 2.2.2 and 2.3).
The profile is then modified by the effects of "specific ordered angular momentum", h, computed from Tidal torque theory (TTT) that describes how larger scales' tidal torques induce smaller scales' angular momentum [94,[116][117][118][119], as well as from "random angular momentum", j, which can be computed from particle orbits, specified by their orbital eccentricity, encoded in the ratio e = r min r max [120], using their pericentric and apocentric radii, resp. r min and r max . The dynamical state of the system induces a correction to the eccentricity [92], computed from the halo maximum and spherically averaged turnaround radii, r ta = x m (x i ) and r max < 0.1r ta , as These effects, as well as DFs, as described by the DF force equation of motion, see Equation A14 in [27], and AC's steepening, following [98], result in the final density profile.

Inclusion of the Baryonic Discs and Clumps Effects
In spiral galaxies, the baryon gas halo evolves into a rotationally supported, stable disk, following the equation of motion to obtain a realistic disc size and mass that also solves the angular momentum catastrophe (AMC) problem, Section 3.

Clump Size Calculation
When disks grow denser, instability appears due to Jean's criterion in spite of shear effect stabilization. This instability leads to clump formation, whose condition was described by [121], defining a limit criterion composed with σ ( 20-80 km/s, for most galaxies hosting a clump), the 1D disk's velocity dispersion, Σ, its surface density, connected to its adiabatic sound speed c s , Ω, its angular velocity and κ, its epicyclic frequency The dispersion relation for perturbations, dω 2 /dk = 0, yields solutions with the mode that grows fastest [122] k inst = πGΣ c 2 s , for Q < 1 (Equation 6 in [65]) that leads to obtain the galaxies clump radii [123] R 7GΣ/Ω 2 1kpc .
The total mass of maximal-velocity-dispersion, marginally unstable discs (Q 1) reaches 3 times that of the cold disc and can convert 10 % of their disk mass M d into clumps [124]. The main properties of clumps found by [125] reflect the general case. For instance, at z 2, 5 × 10 11 M haloes harbor objects of order 10 10 M that remain marginally unstable for 1 Gyr. Smaller haloes have their profiles more efficiently flattened by clumps to DM transfer of AM and energy, as found by [27,[59][60][61][62][63][64][65]126].
Very gas-rich discs are expected to give birth, through self-gravitation instability from the accreting dense gas radiative cooling, to those clumpy structures e.g., [125,127,[141][142][143]. Clump lifetime is central to their impact on halo central density: they can flatten a cusp into a core if they survive stellar feedback disruption long enough to sink to the galactic center. The assessment of a stellar clump's bound state can be obtained through e f , its mass fraction loss by stellar feedback, and ε = 1 − e f , its stellar mass fraction. Then, a threshold at ε ≥ 0.5 marks the limit for most of the stellar clump's mass to remain bound, as confirmed by simulations and analytical models [144]. The efficiency of stellar radiation feedback is evaluated from the expulsion fraction A large sample of densities, sizes, environments, and scales were used in Ref. [145] to obtain e f f 0.01. As typical clumps, having masses M 10 9 M , exhibit e f = 0.15 and ε = 0.85, their mass loss after reaching the galactic halo center should remain small. The expulsion fraction method and this conclusion are unfortunately valid in smaller galaxies, for smaller, more compact clumps, for which no clump can be disrupted before reaching the center.
An alternative method produces clump disruption as a result of a comparison between clump migration time to the structure center and their lifetime. The latter has been studied extensively. In their hydrodynamical simulations exhibiting rotational supported clumps in Jean's equilibrium, Ceverino et al. deduced long lifetimes ( 2 × 10 8 Myr) for them, in agreement with several studies: Ref. [123] found such lifetimes in local systems with Kennicutt-Schmidt law star formation. Such long lifetimes can be explained by a long enough clump galactic center migration time to allow them to retain gas, and form bound star groups, as confirmed by simulations from [146]. Galactic center reaching, long-lived clumps were also produced in other simulations, properly accounting for stellar feedback, e.g., radiative and non-thermal feedback, SNF, radiation pressure, etc., in [129,130,132], also presented by [128], for any reasonable feedback efficiency. Clump ages estimated through metal enrichment, expansion, and gas expulsion time scales, 200 Myr, >100 Myr and 170-1600 Myr, respectively, in [136] also favors strongly long-lived clumps. Finally, clump stability also emerges from similar clump observations, in radius, mass, and in [147][148][149] between low and high redshifts.
DF and TTT balance off to produce migration time; see Equations 1 and 18 of [65,136], obtaining 200 Myr for a 10 9 M clump. The Sedov-Taylor solution, Equations 8 and 9 in [136], produced similar migration and expansion timescales.

Feedback and Star Formation Procedure
The model's stellar feedback proceeds from the prescriptions of Sections 2.2.2 and 2.2.3 in [101,102] for gas cooling, reionisation, star formation, SNF and AGN feedback.
Gas cooling is modeled with a cooling flow, e.g., see Section 2.2.2 in [102,150].
Star formation occurs when gas converts into stars, after settling in a disk. Over a given time interval ∆t that can be set to t dyn , the disc dynamical time, the amount of gas mass converted into stars can be computed as with ψ, the star formation rate, obtained from the mass of gas measured at a density above the threshold n > 9.3/cm 3 , fixed as in [67] as follows, see [101], for more details: SNF explosions inject energy into the halo hot gas, following [152]. The computation of this injected energy is prescribed from a Chabrier IMF [153], consisting of This SN released energy into a reheated disk gas and then compared itself with the reheating energy ∆E hot which that same amount of gas should acquire if its injection in the halo should keep its specific energy constant, that is, if the new gas would remain at equilibrium with the halo hot gas. The amount of disk gas the SN and stellar radiation have reheated, ∆M reheat , since it is all produced from radiation of stellar origin, is proportional to the stellar mass Since the halo hot gas specific energy corresponds to the Virial equilibrium specific kinetic energy V 2 vir 2 , keeping this energy constant under the addition of that reheated gas leads to defining the equilibrium reheating energy as The comparison with the actual energy of the gas injected from the disk into the halo by SNs gives the threshold (∆E SN > ∆E hot ), beyond which gas is expelled, the available energy to expel the reheated gas, and thus the amount of gas ejected from that extra energy Contrary to SNF based models such as [67], our mechanism for cusp flattening initiates before the star formation epoch. Since it uses a gravitational energy source, it is thus less limited in available time and energy. Only after DF shapes the core can Stellar and SN feedback occur, which then disrupts gas clouds in the core, similarly to [65].
AGN feedback points to the effects and the formation of a central Super-Massive-Black-Hole (SMBH). Our approach adopts the SMBH mass accretion, and subsequent AGN feedback models of Booth and Schaye [154], modified by the Martizzi et al. [103,155] prescriptions. When the thresholds 2.4 × 10 6 M /kpc 3 , and 100 km/s, for stellar density and reduced gas density (ρ gas /10), and 3D velocity dispersion, are exceeded, the formation of a seed 10 5 M SMBH occurs and it starts accreting. It has been shown [156] that, above M 6 × 10 11 M , significant AGN quenching occurs.

Confirmations of the Semi-Analytic Model's Robustness
The semi-analytic model's robustness and accuracy was established in several confrontations with observations and previous converging models:
The model predicted, in advance of competing groups in the field, the halo cusp inner slope mass dependence, Figure 2a solid line in [66], in terms of rotation V c , given as 2.8 × 10 −2 M 0.316 vir [157], well in advance of the similar conclusion that one can extrapolate from Figure 6 in [67]; γ.
The inner slope dependence on the ratio between baryonic and total mass of the halo was also predicted by the model; see [30] well before its more publicized claim [67]; δ.
The correct galaxy density profiles were also obtained by the model [27,87] previous to the [25,26] SPH simulations, while that for clusters was predicted in [30], before the results of [69]. Note that these results from the model were obtained with its different dominant mechanism from those of [25,26,69]; ε.
The Tully-Fisher and Faber-Jackson, M − M halo , relationships were compared with simulations from the model (Figures 4 and 5 in [107,108]).

Outline of the Semi-Analytic Model's Main Steps
The semi-analytic approach of our model is much less computing intensive than hydrodynamical and/or N-body simulations, such as NIHAO. Recall that the NIHAO simulation is based on the GASOLINE2 hydrodynamical simulation that includes effects such as [74] photoionisation, Compton cooling, metal cooling, heating from the ultraviolet background, chemical enrichment, star formation, and feedback from massive stars and from SN.
Our model therefore offers a simpler scheme to construct galaxy samples and a faster parameter space exploration. Semi-analytic and N-body/hydro simulations' comparisons have shown good agreement in the studied cases, see [158] and references therein. Our cosmological parameters follow, Section 2 in [74], while the system's baryons start in gas form, with "universal baryon fraction" [114] set as f b = 0.17 ± 0.01, set to 0.167 in [115]. Initial conditions, set from the power spectrum, and evolution follow the description from Appendix B in [27]. The onset, when the system's nonlinear regime is reached, and subsequent evolution of tidal interaction with neighbours is detailed in Appendix C [27]. The generation of random angular momentum during the collapse phase is calculated in Appendix D [27], while the baryonic dissipative collapse is treated in Appendix E [27], when spiral structure or a spheroid shape are generated, as described in Appendix A5 [107]. The depiction of the model's clumps characteristics and formation can be found in Section 2.2 in [107], its star formation is specified in Section 2.3 in [107], while its AGN feedback and BH formation are portrayed in Section 2.3, the final part in [107].

The Dekel-Zhao Profile
By means of simulated haloes using the NIHAO suite of simulations [159], Dekel et al. [81] showed that the functional form of Equation (1) with b = 2 and g = 3.5 gives very good fits to haloes characterized by cores or cusps. As already reported, this parameterization is dubbed DZ. It can be characterized by better fits to simulated profiles with respect to the NFW and the Einasto profile, and better capture of the cores of density profiles. Moreover, the DZ parameterization allows for analytic expressions of the velocity dispersion, the density, the mass, and the gravitational potential. We now recall several important relations.
The inner logarithmic slope, at scale r 1 , is given by [83] s 1 = a + 3.5c 1/2 (r 1 /R vir ) 1/2 1 + c 1/2 (r 1 /R vir ) 1/2 , (15) and the concentration parameter smoothed at the same scale, by The condition that the density is positive yields the condition a ≤ 3, while, from the condition that the inner slope is positive, we obtain the extra condition a + 3.5c 1/2 (r 1 /R vir ) 1/2 ≥ 0. The parameters (a, c) are related to (s 1 , c 2 ). Both a and c can be expressed in terms of s 1 and c 2 as and As in [83], (a, c) are used to express analytic expressions, and (s 1 , c 2 ) are used in numerical tests. The analytic expression for the gravitational potential, the velocity dispersion, and the lensing properties are given in Equation (18), Equation (21), and Section 2.3 of [83], respectively.
Freundlich et al. [83] compared the DZ, the gNFW, and the Einasto profiles using the results of the NIHAO suite of SPH simulations. The authors used some profiles from the NIHAO suite, plotted in their Figure 3, and fitted the NIHAO profiles using the gNFW, DZ, and Einasto profiles. The parameters for each profile were considered free. Freundlich et al. [83] showed that the DZ performs better than the other two parameterizations, and in particular, significantly better than the Einasto profile, and marginally better than the gNFW profile.

The Dekel-Zhao Mass Dependent Profile
As already discussed, Ref. [86] fitted a profile, similar to Equation (1), to SPH simulations, finding the functional form of the three parameters in terms of M /M vir , and Ref. [84] proceeded similarly with the Einasto profile. The result provided a density profile with parameters depending on mass. Freundlich et al. [83] did the same with the DZ profile. In summary, they fitted the logarithm of the density profile of the simulated haloes given by NIHAO, according to the DZ parameterization through a least-square minimization in the range 0.01 ≤ R vir ≤ 1. Given the mass M vir of a profile, one can easily find R vir , and then one remains with two free parameters a, and c. One can proceed similarly with M and M /M vir . The profile radii r were spaced logarithmically, with N 100 radii in the indicated range (0.01 R vir -R vir ). As we already discussed, the parameters s 1 and c 2 are related to a and c. The quality of the fit was evaluated using the rms of the residuals between the simulated log ρ, and the model With the procedure described above, one can obtain the best fit s 1 , and c 2 , in terms of M vir , M , and M /M vir . In Figure 8 of [83], the dependence of s 1 and c 2 is plotted in terms of M vir , M , and M /M vir , derived from the DZ density profile fits.
To capture the behavior of s1 and c 2 in terms of M vir , M , and M /M vir , Freundlich et al. [83] assumed two functions where x 0 , s , s , and ν are adjustable parameters whose values are reproduced in Table 1 of [83]. Concerning c 2 , they used where again x 0 , c , and ν are adjustable parameters. With this procedure, the DZ becomes a mass-dependent profile, displaying some advantages over competing profiles. For example, compared with that of Di Cintio et al. [86], it requires less free parameters and has analytic expressions for dispersion velocity, gravitational potential, and lensing properties.
Until now, we have described how [83] obtained a mass-dependent profile. In our case, we followed a similar procedure, with the difference that, while the simulated haloes' logarithm of the density profile was fitted by Freundlich et al. [83] from the NIHAO simulations, according to the DZ parameterization, we used the results given by our model. The procedures are identical, except for the forms of the functions that capture the behavior of s1, and c 2 in terms of M vir , M , and M /M vir . From our model, we used the functions s 1 (M ) or c 2 (M ) = c 0 + c 1 y + c 2 y 2 + c 3 y 3 + c 4 y 4 10 4 ≤ M ≤ 10 12 y = log 10(1 + (M /x 0 ) n ), while s 1 (M vir ) = c 0 + c 1 y + c 2 y 2 + c 3 y 3 + c 4 y 4 + +c 5 y 5 10 9 ≤ M vir ≤ 10 15 y = log 10(1 + (M vir /x 0 ) n ) and c 2 (M vir ) = c 0 + c 1 y + c 2 y 2 + c 3 y 3 + +c 4 y 4 10 10.3 ≤ M vir ≤ 10 15 y = log 10(1 + (M vir /x 0 ) n ).
The parameters of the fit can be found in Table 1.  Figure 1 shows the dependence of s 1 and c 2 on M and M vir . The shaded gray region represents the rms of the residuals. As described in Section 3.3.1 in [83], it is obtained by means of an iterative process excluding points beyond 3 σ. s 1 and c 2 are obtained from the fits to the density profile of M vir and M . The slope s 1 is obtained using Equations (22) and (23). The black solid line represents the best-fit curve, while the dashed lines represent Equations (45) and (46) of [83]. Similarly, parameter c 2 is obtained using Equations (22) and (24). The best-fit curve is again represented by the black solid line. The relations of s 1 or c 2 in terms of M /M vir are not plotted, since it is directly related to the previous two relations (dependence on M vir and M ). For a fixed value of s 1 , one can obtain M , and M vir , and consider their ratio. Moreover, the relation between s 1 , c 2 , and M /M vir , because of AGN feedback, is characterized by a two-branch relation, as that shown in Figure 3 in [74]. In order to obtain the DZ DM profile relative to each halo, one can follow the prescription given in the following, also proposed in Section 4.1 in [83].
The parameters s 1 and c 2 , in terms of M vir , M , are given by Equations (22)- (23) with the parameters given in Table 1.
The parameters a and c of the DZ function can be obtained from the values of s 1 and c 2 using Equations (17) and (18). (5) The scaling parameter ρ c is given by the following relations: µ = c a−3 (1 + c 1/2 ) 2(3−a) , ρ vir = 3M vir /4πR 3 vir = ∆ρ crit , and ρ c = (1 − a/3)c 3 µρ vir , where r c = R vir /c. (6) Finally, the DZ mass-dependent density profile is given by Equation (14). The corresponding circular velocity profile can be obtained from Equations 4 and 6 of [83] with ρ c = c 3 µρ vir , b = 2, and g = 3.  (22) and (23). The best-fit curve is represented by the black solid line, while the dashed lines represent Equations (45) and (46) of [83]. The parameters of the fit are given in Table 1. The parameter c 2 is obtained using Equations (22) and (24).The best-fit curve is represented by the black solid line. The rms σ in gray is obtained, as discussed in Section 5, by means of an iterative process excluding points beyond 3 σ.
An application of the method is shown in Figure 2. The plot shows the density profiles for given values of the masses (M vir , M , and M /M vir ) taken from the NIHAO suite. The red lines are the profiles from NIHAO, the short dashed lines those obtained by [83] with the DZ mass-dependent profile, the dotted ones the result from Di Cintio et al. [86], and the long-dashed lines are the density profiles that we obtained. Similarly to [83], s 1 is obtained from Equations (22)- (23), while c 2 , and r c are allowed to vary, and ρ c is constrained from the halo mass M vir . In the case of low masses, our density profiles are flatter at low radii. This is due to the fact that, in the NIHAO simulations, the role of the supernovae feedback is larger than in our model. As we discussed in Section 2, in our case, the flattening of the profiles is mainly due to the interaction of gas clumps with DM through dynamical friction.

Clusters of Galaxies' Density Profiles
In Figure 2, we showed the DZ density profile in the case of galaxy types from dwarfs to the Milky Way. Freundlich et al. [83] dedicated their papers to show how the DZ mass dependent profile can be used to describe density profiles in better detail than the gNFW and Einasto, and furthermore compared it with [86]. The main goal of the present paper is to extend [83] to the mass regime of galaxy clusters. To this aim, in Figure 1, we obtained the s 1 and c 2 parameters in terms of M , M vir , up to masses of the order of 10 15 M . In the following, we show how the DZ mass dependent profile can fit the density profile of clusters of galaxies. We choose the [16,17] clusters. As shown in Figure 3 of [17], the mass distribution was studied combining information from weak and strong lensing, and stellar kinematics, allowing the reconstruction of the mass distribution in baryons in the BCG. As shown in Figure 3, the total mass of the clusters studied (MS2137, A963, A383, A611, A2537, A2667, and A2390) is fitted by an NFW profile, while the DM distribution has profiles flatter than the NFW profile, as seen in Figure 5 of [17]. The baryons usually dominate inside a radius of 10 kpc. Concerning the baryon content, Table 3 of [16] gives the luminosity L V and Table 4 the ratio between stellar mass and luminosity, Υ = M * /L V . The product of Υ by L V yields the stellar mass. Table 8, of the same paper, describes the characteristics of the total mass density profile, the virial mass, and radius, given ∆ = 200, M 200 and r 200 . The results obtained by means of X-ray observations are also given. We summarize in Table 2 the characteristic of the clusters, starting with the inner slope calculated around 1 kpc, which corresponds to 0.1 % of the virial radius. We then report L V , Υ, the stellar mass, the virial mass, the ratio between stellar and virial mass, and finally the values of s 1 and c 2 obtained for each cluster with the method described in the previous sections. In Figure 3, we compared the DM density profile with DZ mass dependent profile calculated as described in the previous sections. The DM density profile obtained by [16,17] is represented in blue, while the mass dependent DZ prediction is plotted in black. The width of the blue bands represents the uncertainty, including the 1 σ uncertainties for isotropic models, see [16,17], and a systematic component obtained as described in Section 4.3 of [17]. As shown, there is a good agreement between the DM density profile and the mass dependent DZ profile.

III
Dark matter (our) Figure 3. Comparison of the DM density profile obtained by [16,17] and the density profile predicted by the mass dependent DZ density profile. The width of the blue bands represents the 1 σ uncertainty plus a systematic component (Section 4.3 of [17]). The black bands represent the DZ mass dependent profile, and the 1 σ uncertainty.

Conclusions
In the present paper, we extended the work of Freundlich et al. [83] related to the DZ density profile. By means of the NIHAO suite of simulations [159], Dekel et al. [81] showed that the functional form of Equation (1) with b = 2 and g = 3.5 gives very good fits to haloes characterized by cores or cusps. This parameterization is dubbed DZ. This function depends on two parameters: Freundlich et al. [83] fitted the DZ profile to SPH simulations expressing the two parameters in terms of stellar mass, M * , virial mass M v , and their ratio. To capture the behavior of the two parameters, s 1 and c 2 , in terms of the quoted masses, two peculiar functions for s 1 , and c 2 were assumed. In this way, the DZ profile acquires a dependence on halo mass. While Freundlich et al. [83] obtained a mass dependent DZ profile for the case of galaxies, we extended his work to clusters of galaxies. Freundlich et al. [83] used the NIHAO simulations to obtain the mass dependence of the two parameters of the profile, while we used our own analytical model. In this model, baryonic clumps moving inside structures exchange angular momentum and energy with DM through dynamical friction. Following Freundlich et al. [83], we obtained the dependence from baryon physics of the two shape parameters. In this way, we obtained a mass dependent DZ profile describing the DM profiles from galaxies to clusters of galaxies. After obtaining the mass dependent DZ profile, we compared our model's simulated profiles with the results of [83,86]. The results from our model are in good agreement with simulations and the corresponding profiles. As already reported, the main goal of the present paper is to extend [83] up to the mass of clusters. To reach this goal, we obtained the s 1 and c 2 parameters in terms of M , M vir , up to masses of the order of 10 15 M . In order to see how the DZ mass dependent profile fit the density profile of clusters of galaxies, we choose the [16,17] clusters. The mass distribution in these clusters was obtained combining several techniques, such as weak and strong lensing, and stellar kinematics. We studied the clusters MS2137, A963, A383, A611, A2537, A2667, and A2390. The DM distribution displays profiles flatter than the NFW profile, as shown in Figure 5 of [17]. Given the parameters of the clusters obtained by [16,17], we obtained values of s 1 and c 2 for each cluster with the method described in the previous sections. We compared the DM density profile obtained by [16,17] with those obtained with our model. As shown, there is a good agreement between the DM density profile and our mass dependent DZ profile.
Author Contributions: All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.