Dark Matter Halos: The Dynamical Basis of Effective Empirical Models

Abstract

We investigate the dynamical basis of the classic empirical models (specifically, Sérsic-Einasto and generalized NFW) that are widely used to describe the distributions of collisionless matter in galaxies. We submit that such a basis is provided by our 𝛼-profiles, shown to constitute solutions of the Jeans dynamical equilibrium with physical boundary conditions. We show how to set the parameters of the empirical in terms of the dynamical models; we find the empirical models, and specifically Sérsic-Einasto, to constitute a simple and close approximation to the dynamical models. Finally, we discuss how these provide a useful baseline for assessing the impact of the small-scale dynamics that may modulate the density slope in the central galaxy regions.

1. Introduction

The classic Sérsic [1] models met a wide and lasting success as empirical representations of the projected (2-dimensional) light distributions in spheroidal galaxies (for a review, see Kormendy et al. [2]). Einasto [3] developed and used a similar shape to describe in simple terms 3-dimensional stellar mass profiles.

On the other hand, recent extensive 𝑁-body simulations (e.g., [4–8]) indicate that the Sérsic and Einasto functional forms also provide good patterns to represent the spherically-averaged mass distributions in dark matter (DM) halos ranging from galaxies to galaxy clusters. These hold at levels comparable to, or even better than the popular NFW formula [9].

Still, no agreed understanding is available to explain the value in both the real and the virtual worlds of the Sérsic and Einasto representations (see [2, 10]). Can we identify the underlying astrophysical basis?

2. Empirical Models

Before addressing the issue, we note that these models belong to two main families: generalized NFW (see [11–13]; hereafter gNFW) and Sérsic-Einasto (see [10, 14, 15]; hereafter SE).

2.1. Density Runs

The density runs of the SE family may be represented in the form𝜌(𝑟)̂=𝑟̂^−𝜏𝑒^−𝑢(𝑟̂^𝜂₋₁₎𝜏,𝑢=2−𝜂. (1) Here, quantities are normalized to their value at 𝑟₋₂, the reference radius where the logarithmic slope 𝛾≡−dlog𝜌/dlog𝑟 takes on the value 2; typically, in nearby elliptical galaxies, 𝑟₋₂ corresponds to sizes of order 10 kpc, a few times the half-light radius 𝑅_𝑒.

The parameters 𝜏 and 𝜂 describe the inner slope and the middle curvature of the density run, respectively. The original Einasto profile belongs to this family and is obtained when 𝜏=0. Note, however, that by deprojecting from the plane of the sky, a Sérsic 2-dimensional run 𝑒^{−𝑠1/𝑛} with index 𝑛≈3-4 (suited for normal ellipticals, see [2]) produces a cuspy inner run as in (1) with 𝜏⋍1−1.19/2𝑛+0.22/4𝑛²≈0.8 significantly different from 0 and less than 1, as shown by Prugniel and Simien [15].

On the DM side, recent simulations (see [4, 6, 8]) only provide an upper bound 𝜏<0.9 for the inner slope. When the original Einasto profile (with 𝜏=0) is adopted, the best-fit to simulated DM halos obtains for 𝜂≈0.2; we will come back to this value later on.

In turn, the density runs of the gNFW family may be written in the following form(In the literature these runs are sometimes referred to as 𝛼𝛽𝛾-models, and equivalently defined via the parameters 𝛾=𝜏, 𝛼=1/𝜂, and 𝛽=𝜏+𝜂𝜉.):𝜌(𝑟)̂=𝑟̂^−𝜏(1+𝑤1+𝑤𝑟^𝜂̂)^𝜉,𝑤=−2−𝜏;2−𝜏−𝜂𝜉 (2) the parameters 𝜏, 𝜂, and 𝜉 describe the central slope, the middle curvature, and the outer decline of the density run, respectively. Note that familiar empirical profiles are recovered for specific values of the triple (𝜏,𝜂,𝜉); for example, Plummer's [16] corresponds to (0,2,2.5), Jaffe's [17] to (2,1,2), Herquist's [11] to (1,1,3), and NFW to (1,1,2).

2.2. Toward a Single Family

The main apparent difference between SE and gNFW is constituted by the former's exponential decline versus the latter's powerlaw falloff 𝜌∝𝑟^{−(𝜏+𝜂𝜉)} for large 𝑟.

On the other hand, (2) is to be considered for large values of 𝜉 anyway, since a steep density run in the halo outskirts is indicated by observations of light distribution in spheroidal galaxies (other than cDs, see [2]) and of DM distributions from weak lensing in galaxies and galaxy clusters (e.g., [18–20]).

The circumstance is easily translated into the formal statement that the gNFW family converges to the SE for large 𝜉. This is seen on recasting 𝜌𝑟^𝜏̂ from (2) in exponential form, to read𝑒^{𝜉ln[(1+𝑤)/(1+𝑤}𝑟̂^𝜂_)]⋍𝑒^{𝜉𝑤(1−𝑟𝜂)}⋍𝑒^−𝑢(𝑟̂^𝜂₋₁₎; (3) for approximating the middle and last terms, we have used the circumstance that 𝜉≫1 implies 𝑤≫1 and so 𝜉𝑤⋍(2−𝜏)/𝜂≡𝑢 applies. Thus the two families in (1) and (2) actually become one in this limit.

Thus in the following we focus mainly on the SE family, and proceed to discuss its dynamical basis in terms of the Jeans equation.

3. The Dynamical Model

The dynamical model of DM halos hinges upon the radial Jeans equation that expresses the self-gravitating, equilibrium of collisionless matter (see [21]). The Jeans equation reads1𝜌d(𝜌𝜎²_𝑟)=d𝑟−𝐺𝑀(<𝑟)𝑟²−2𝛽𝜎²_𝑟𝑟, (4) in terms of the density 𝜌(𝑟), the related cumulative mass 𝑀(<𝑟)≡4𝜋∫^𝑟₀d𝑥𝑥²𝜌(𝑥), and the radial velocity dispersion 𝜎²_𝑟(𝑟). The last term on the r.h.s. describes the effects of anisotropic random velocities via the standard Binney [22] parameter 𝛽≡1−𝜎²_𝜃/𝜎²_𝑟.

Note that the Jeans equation is designed to describe a (quasi-)static equilibrium, away from extreme major merger events like in the case of the Bullet Cluster (see [23]). But even in relaxed conditions, solving Jeans requires an “equation of state,” that is, a functional relation expressing the DM pressure 𝜌𝜎²_𝑟 in terms of density (and possibly radius) only.

3.1. Equation of State

In seeking for such a relation, one can make contact with the classic theory of the nonlinear collapse for DM halos in an expanding Universe; here self-similar arguments play the role of a pivotal pattern (see [24–26]). This still applies to modern views of the halo development (e.g., [27–31]), that comprise two stages: an early collapse that builds up the halo main body via a few major merger events and sets its phase-space structure by dynamical relaxation of DM particle orbits; this tails off into a secular development of the outskirts by smooth accretion and minor mergers.

The essence of the macroscopic equilibrium is conveyed by the self-similar scaling 𝜎²_𝑟∝𝑀/𝑟 adding to the geometric relation 𝜌∝𝑀/𝑟³. The macroscopic import of the halo phase-space structure is conveyed by combining these two quantities into the “phase-space density” 𝜌/𝜎³_𝑟 or equivalently into the functional 𝐾(𝑟)≡𝜎²_𝑟/𝜌^2/3 often referred to as DM “entropy” (see [24, 26]). For the latter quantity, one easily derives the scaling 𝐾(𝑟)∝𝑟𝑀^1/3 implying𝐾(𝑟)∝𝑟^𝛼, (5) whence one expects a slope 𝛼 slightly exceeding 1.

To focus the values of 𝛼, in Lapi and Cavaliere [27, 28] we have developed a full semianalytic treatment of the halo growth in the standard accelerating Universe (see [32]). We found constant values of 𝛼, that fall within the narrow range 1.25–1.3; on average, such values grow weakly with the mass of the halo body, from galaxies to rich clusters.

The halo development process has been probed, and the two-stage view confirmed by intensive 𝑁-body simulations (e.g., [33–39]). These also confirm that (i) a (quasi-)static macroscopic equilibrium is attained at the end of the fast collapse and is retained during the subsequent stage of secular, smooth mass addition; (ii) a persistent feature of such an equilibrium is constituted by powerlaw correlations holding in the form 𝜎²_𝑟/𝜌^2/3∝𝑟^𝛼, although it is still widely debated whether the radial or the total velocity dispersion best applies (see also the discussions by [6, 36]).

In building up our dynamical models, we focus on the quantity 𝐾≡𝜎²_𝑟/𝜌^2/3∝𝑟^𝛼 that involves the radial dispersion 𝜎²_𝑟 (see also [40]). Operationally, this provides a direct expression for the radial pressure term 𝜌𝜎²_𝑟=𝐾𝜌^5/3∝𝑟^𝛼𝜌^5/3 in the Jeans (4); anisotropies are accounted for by the last term on the r.h.s., as discussed in Section 3.3 below.

3.2. The DM 𝛼-Profiles

In terms of 𝐾(𝑟)∝𝑟^𝛼, the Jeans equation may be recast into the compact form3𝛾=5𝑣(𝛼+^2𝑐𝜎²_𝑟+2𝛽), (6) with 𝛾≡−dlog𝜌/dlog𝑟 representing the logarithmic density slope and 𝑣²_𝑐≡𝐺𝑀(<𝑟)/𝑟 the circular velocity. Remarkably, by double differentiation, this integrodifferential equation for 𝜌(𝑟) reduces to a handy second order differential equation for 𝛾 (see [40, 41]).

Tackling first the isotropic case 𝛽=0, we recall that the solution space of (6) spans the range 𝛼≤35/27=1.296: the specific solution for the upper bound, and the behaviors of the others have been analytically investigated by Taylor and Navarro [26], Austin et al. [41], and Dehnen and McLaughlin [40]. In Lapi and Cavaliere [28], we explicitly derived the solutions in the full range 𝛼=1.25−1.296, that are marked by a monotonically decreasing run and satisfy physical boundary conditions: a finite central pressure or energy density (equivalent to a round minimum of the gravitational potential) and a steep outer run implying a finite and rapidly converging (hence a definite) overall mass. We dubbed 𝛼-profiles these physical solutions.

We shall use the following basic features of the latter. In the halo body at the point 𝑟₀, the 𝛼-profile is tangent to the pure powerlaw solution 𝜌∝𝑟^−𝛾0 of the Jeans equation; there 𝑣²_𝑐/𝜎²_𝑟∝𝑟^{2−𝛾0/3−𝛼} applies, to imply from  (6)𝛾₀=6−3𝛼. (7) This is consistent for 𝛼=1.25 with the self-similar slope. Note that the point 𝑟₀ lies in the neighborhood of the radius 𝑟₋₂ (see Section 2.1); specifically 𝑟₀≈1.74-1.51𝑟₋₂ holds for 𝛼≈1.25-1.3.

On the other hand, a monotonic density run implies the term 𝑣²_𝑐/𝜎²_𝑟∝𝑟^{2−𝛾/3−𝛼} to vanish at the center; this results in an inner powerlaw 𝜌∝𝑟^−𝛾𝑎 with the slope𝛾_𝑎=35𝛼. (8) This differs from zero as long as the entropy run grows from the center following 𝐾∝𝑟^𝛼 with 𝛼>0.

Finally, a finite mass implies 𝑣²_𝑐/𝜎²_𝑟∝𝑟^{−1+2𝛾/3−𝛼} to hold in the outskirts, so as to yield a typical outer decline 𝜌∝𝑟^−𝛾𝑏 with slope𝛾_𝑏=32(1+𝛼). (9) This exceeds the value 3 and so constitutes the hallmark of a rapidly saturating mass; the circumstance makes less compelling here the role of a virial boundary.

Thus, compared to NFW the inner slope of the dynamical model is considerably flatter and the outer slope steeper; compared to the original Einasto profile, the main difference occurs in the inner regions where the dynamical model is (moderately) steeper.

3.3. Anisotropy

It is clear from (6) that anisotropies will steepen the density run for positive 𝛽, and flatten it for negative 𝛽. The latter condition is expected to prevail in the inner region, where tangential components develop from the angular momentum barrier [30, 42]. Moving outwards, radial motions are expected to prevail, so raising 𝛽 up to values around 0.5 at 𝑟≈𝑟₀; outwards of this, 𝛽 is expected to saturate or even decrease, as one enters a region increasingly populated by DM particles on eccentric orbits with vanishing radial dispersions at their apocenters (see [24]).

This view is supported by numerical simulations (see [6, 40, 41, 43]), which in detail suggest the average anisotropy-density slope relation𝛽(𝑟)≈𝛽(0)+𝛽^′[𝛾(𝑟)−𝛾_𝑎], (10) to hold with parameters 𝛽(0)≈−0.1, 𝛽^′≈0.2, and the constraint 𝛽(𝑟)≲0.5. Note that at all radii the inequality 𝛾(𝑟)≥2𝛽(𝑟) is satisfied; this has been conjectured to constitute a necessary condition for a self-consistent spherical model with positive distribution function (see [44]).

In Lapi and Cavaliere [27], we extended the dynamical model to such anisotropic conditions in the full range 𝛼≈1.25-1.3, inspired by the analysis of Dehnen and McLaughlin [40] for the upper bound of 𝛼. We note that the latter is now slightly modified to 35/27−4𝛽(0)/27≈1.31; likewise, the point 𝑟₀ where 𝛾=𝛾₀=6−3𝛼 applies moves slightly inwards, so that 𝑟₀≈1.58−1.38𝑟₋₂ now holds.

The main outcome, however, is that the density profile is somewhat flattened at the center relative to the isotropic case; the inner slope now reads𝛾_𝑎=356𝛼+5𝛽(0). (11) In particular, even a limited central anisotropy (corresponding to values 𝛽(0)≈−0.1) causes an appreciable flattening down to 𝛾_𝑎≈0.63-0.66 for 𝛼≈1.25-1.3.

On the other hand, we stress that such small phenomenological anisotropies near the center imply the radial 𝜎²_𝑟 and the total dispersions 𝜎²=𝜎²_𝑟[1−2𝛽/3] to be very close.

4. From Dynamical to Empirical Models

Here we discuss how the parameters of the empirical profiles (see Section 2) can be set based on our dynamical model (see Section 3); in such conditions, it will turn out that such profiles constitute close approximations to the model over a wide radial range.

4.1. Parameters from Dynamics

First we consider the original Einasto profile (𝜏=0 in (1)), since this has been widely used in the context of DM halo simulations. Here, 𝜏 is fixed to 0, and the only free parameter is the curvature 𝜂. This we set by requiring the logarithmic density slope𝛾(𝑟)̂=2𝑟^𝜂̂ (12) to equal 𝛾₀ at the point 𝑟₀̂. So we find the expression𝜂=log(𝛾₀/2)log𝑟₀̂, (13) that takes on values 𝜂≈0.15-0.2; see Tables 1 and 2; remarkably, these turn out to agree with those derived from fitting the outcomes of state-of-the-art 𝑁-body simulations in terms of the same Einasto density run, as performed by Navarro et al. [6].

Table 1 

Parameters of empirical models in the isotropic case.

Table 2 

Parameters of empirical models in the anisotropic case.

On the other hand, the flat central slope of the Einasto profile is at variance with the value 𝛾_𝑎=3𝛼/5 given by our dynamical models based on Jeans; to wit, consistency between pure Einasto and Jeans would require at the center a flat entropy distribution and a vanishing pressure.

Actually, the simulations quoted in Section 3.1 within their finite mass resolution provide only an upper limit 𝜏<0.9 to the central slope. This grants scope to the full SE family of (1).

The latter features two parameters, the inner slope 𝜏 and the middle curvature 𝜂. These we set by requiring the logarithmic density slope𝛾(𝑟)̂=𝜏+(2−𝜏)𝑟^𝜂̂ (14) to equal 𝛾_𝑎 for 𝑟̂→0, and 𝛾₀ at 𝑟₀̂; so we find𝜏=𝛾_𝑎,𝜂=log(𝛾₀−𝛾_𝑎)/(2−𝛾_𝑎)log𝑟₀̂. (15) Thus we predict the central slope to take on values 𝜏=0.6-0.8 and the corresponding curvature parameter to take on values 𝜂=0.2-0.3; see Tables 1 and 2. It will be worth fitting the outcomes of 𝑁-body simulations based on these extended SE profiles with 𝜏>0.

Finally, we report the corresponding results for the empirical gNFW family. This features three parameters: inner slope 𝜏, middle curvature 𝜂, and strength of the outer decline 𝜉; these we set by requiring the logarithmic density slopê𝑟̂𝛾(𝑟)=𝜏+𝜂𝜉[(2−𝜏)/(2−𝜏−𝜂𝜉)]^𝜂[(2−𝜏)/(2−𝜏−𝜂𝜉)]𝑟^𝜂̂−1 (16) to equal 𝛾_𝑎 for 𝑟̂→0, 𝛾₀ at 𝑟₀̂, and 𝛾_𝑏 for 𝑟→∞. So we find𝜏=𝛾_𝑎,𝜂=log[(𝛾₀−𝛾_𝑎)(2−𝛾_𝑏)/(𝛾₀−𝛾_𝑏)(2−𝛾_𝑎)]log𝑟₀̂,𝛾𝜉=^𝑏−𝛾_𝑎𝜂. (17)

The parameters so determined are listed in Tables 1 and 2 for both the isotropic and the anisotropic conditions.

4.2. Results and Comments

With the parameters focused as discussed in the previous subsection, Figure 1 illustrates how the empirical compare with our dynamical models. We plot the density run of the latter (specifically for the 𝛼-profile with 𝛼≈1.25 suitable to galactic halos), compared to those of the Einasto, SE, and gNFW models. The left and right panels refer to isotropic and anisotropic conditions, respectively; the popular NFW profile is also shown for reference. To make comparisons easier, we plot in the lower panels the corresponding logarithmic density slopes.

(a)

(b)

Figure 1 

Density profiles in isotropic (a) and anisotropic conditions (b) for the dynamical model and its approximation in terms of empirical models with the parameters derived in Section 4 and summarized in Tables 1 and 2. The lower panels highlight the corresponding logarithmic density slopes.

It turns out that the closest approximation to the dynamical model is provided by SE, which shares with it not only the central slope by construction, but also the body and the outer behaviors. The original Einasto profile provides an acceptable approximation in the middle and outer ranges, but not at the center, because of its flatness. On the other hand, the gNFW family provides an acceptable approximation in the inner and middle ranges, but not in the outskirts where its slope is too flat. Finally, the NFW profile provides an acceptable approximation in the middle range.

Similar conclusions concern the profiles of circular velocities 𝑣²_𝑐(𝑟)≡𝐺𝑀(<𝑟)/𝑟, that are analytically dealt with in the Appendix and illustrated in Figure 2.

(a)

(b)

Figure 2 

Same as Figure 1 for the profiles of circular velocity, and for the corresponding logarithmic slopes.

We stress that the handy SE representation is convenient in analyzing data in several contexts, including the DM particle annihilation signal expected from the Galactic Center (see [45]), rotation curves of dwarf and normal spiral galaxies (see [46]); individual and statistical properties of elliptical and spiral galaxies (see [47]), strong and weak gravitational lensing (see [27])—currently observed in clusters (e.g., [48]) and soon in massive elliptical galaxies (see discussion by [49]) and X-ray emission from the intracluster plasma (see [50–52]).

5. Discussion

We first stress that the dynamical model (as well as its approximations in terms of empirical models) is in keeping with the basic features of standard DM, that is, its cold and collisionless nature. In fact, it implies 𝜎²_𝑟(𝑟)→0 for large 𝑟≫𝑟₋₂, a behavior expected in the outskirts for cold matter dominating the potential well.

At the inner end, with decreasing 𝑟, we expect 𝜎²_𝑟(𝑟) to increase toward a maximum, corresponding to effective conversion of inflow kinetic into random energy. In fact, toward the center Jeans requires dlog𝜎²_𝑟/dlog𝑟=𝛾−𝐺𝑀(<𝑟)/𝑟²→𝛾_𝑎 to hold as the gravitational force vanishes there, to the effect that 𝜎²_𝑟(𝑟)∝𝑟^𝛾𝑎→0.

Concerning the collisionless nature of the DM, the boundary conditions at the center imply a finite, non-zero pressure (and energy density), while a long collisional mean free path allows the pressure gradient d𝑝/d𝑟 to diverge. Conversely, with a short mean free path 𝜆 the pressure gradient cannot diverge on scales 𝑟≳𝜆, where a finite 𝜎² and a flatter 𝛾 apply. In fact, weakly collisional conditions have been proposed to explain the cored light profiles observed in many spheroidal galaxies (see [53]).

On approaching the center of a galactic halo, one expects the basic dynamical model from large-scale Jeans equilibrium to be altered to an increasing degree by small-scale dynamics and/or energetics related to baryons. These processes are specifically related to following issues: transfer of energy/angular momentum from baryons to DM during galaxy formation, scouring baryons by the energy feedback from central active galactic nuclei and any “adiabatic” contraction of the baryons. Such issues will be briefly discussed in turn, with a warning that they enter increasingly debated grounds.

5.1. Energy/Angular Momentum Transfers

Flattening of the inner density profile may be caused by transfer of energy and/or angular momentum from the baryons to the DM during the galaxy formation process (see [54–56]).

In detail, upon transfer of tangential random motions from the baryons to an initially isotropic DM structure, the density in the inner region is expected to behave as [56]𝜌∝𝑟^{−2[𝛾𝑎+2(2−𝛾𝑎)𝛽]/[2+(2−𝛾𝑎)𝛽]}. (18) Thus for 𝛽<0 the profile is flattened relative to the original 𝛾_𝑎, down to the point of developing a core for 𝛽≲−𝛾_𝑎/2(2−𝛾_𝑎)≈−0.3.

However, a reliable assessment of the amount of angular momentum transferred from the baryons to the DM is still wanting, and would require aimed numerical simulations of better resolution than presently achieved.

5.2. Other Processes on Inner Scales

Less agreed processes may affect galactic scales 𝑟≲10² pc, for example, at the formation of a spheroid, central starbursts and a supermassive black hole may easily discharge enough energy (∼10⁶² erg for a black hole mass 𝑀_•∼10⁹𝑀_⊙) with sufficient coupling (≳1%) to blow most of the gaseous baryonic mass 𝑚∝𝑟^3−𝛾𝑎/(3−𝛾_𝑎) out of the gravitational potential well. This will cause an expansion of the DM and of the stellar distributions (see [57]), that flattens the central slope.

In addition, binary black hole dynamics following a substantial merger may eject on longer timescales formed stars from radii 𝑟≈10(𝑀_•/10⁸𝑀_⊙)^{1/(3−𝛾𝑎)} pc containing an overall mass of a few times the black hole's, and so may cause a light deficit in some galaxy cores (see [2, 58, 59]). A full discussion of the issue concerning cored versus cusped ellipticals is beyond the scope of the present paper.

5.3. Adiabatic Contraction?

On the other hand, some steepening of the inner density profile may be induced by any “adiabatic” contraction of the diffuse star-forming baryons into a disc-like structure, as proposed by Blumenthal et al. [60] and Mo et al. [61] but currently under scrutiny, see Abadi et al. [62].

On the basis of the standard treatments, it is easily shown that in the inner region an initial powerlaw 𝜌(𝑟_𝑖)∝𝑟_𝑖^−𝛾 is modified into𝜌∝𝑟^{−3/(4−𝛾𝑎)}; (19) this yields typical slopes around 0.9, steeper than the original 𝛾_𝑎≤0.78 but still significantly flatter than 1.

However, recent numerical simulations (see discussion by [62]) suggest that the treatment of adiabatic contraction leading to (19) is likely to be extreme; actually, in the inner region the contraction is ineffective and the density slope hardly modified. Again, highly resolved 𝑁-body experiments are needed to clarify the issue.

6. Conclusions

We have discussed the dynamical basis of the Sérsic-Einasto empirical models, in terms of well-behaved solutions of the Jeans equation with physical boundary conditions comprising: a finite central energy density, a closely self-similar body, and a finite (definite) overall mass.

We find the SE profile to be particularly suitable to accommodate the general run of the dynamical solution. Specifically, we have discussed how to tune the parameters of SE in terms of the dynamical model; in such conditions, we find the former to constitute a simple and close approximation to the latter.

The resulting SE profile shares with the dynamical model the following features: an outer steep decline, hence a definite overall mass, a closely self-similar body with slope 𝛾₀=6−3𝛼, and an inner slope around 𝛾_𝑎=3𝛼/5, hence flatter than −1. The latter slope provides an useful baseline for discussing alterations of the inner behavior caused by additional baryonic processes.

In conclusion, we submit that the dynamical models discussed here, namely, the 𝛼-profiles, provide the astrophysical basis for understanding the empirical success of the SE profiles in fitting the real and the virtual observables, from galaxies to galaxy clusters.

Acknowledgments

The paper is supported by ASI and INAF. It is a pleasure to acknowledge inspiring correspondence with A. Graham, and stimulating discussions with L. Danese, P. Salucci, and Y. Rephaeli. A. Lapi thanks SISSA and INAF-OATS for warm hospitality.