Physically Motivated Fit to Mass Surface Density Profiles Observed in Galaxies

Galaxies are self-gravitating structures that, among all possible equilibrium configurations, choose only those consistent with a stellar mass surface density profile approximately resembling a Sérsic function (e.g., Sersic 1968; Caon et al. 1993; Trujillo et al. 2001b; Blanton et al. 2003; Graham & Driver 2005; van der Wel et al. 2012). ⁴ Settling into this particular configuration could be due to either some fundamental physical process (as, e.g., thermodynamic equilibrium sets the velocities of the molecules in the air) or to the initial conditions that gave rise to the system (Binney & Tremaine 2008). The mass distribution in galaxies is commonly explained as the outcome of cosmological initial conditions (Cen 2014; Nipoti 2015; Ludlow & Angulo 2017; Brown et al. 2020). The option of a fundamental process like thermodynamic equilibrium determining the configuration is traditionally discredited because of two conceptually very different reasons.

If the structure were set by thermodynamic equilibrium, following the principles of statistical physics, it should correspond to the most probable configuration of a self-gravitating system, and thus, it should result from maximizing its entropy. The use of the classical Boltzmann–Gibbs entropy to characterize self-gravitating systems leads to a distribution with infinite mass and energy (Binney & Tremaine 2008; Padmanabhan 2008), disfavoring the thermodynamic equilibrium explanation. However, this difficulty of the theory has been overcome as follows. In the standard Boltzmann–Gibbs approach, the long-range forces that govern self-gravitating systems are not taken into account. Systems with long-range interactions admit long-lasting metastable states described by a maximum entropy formalism based on Tsallis (S_q ) nonadditive entropies (Tsallis 1988, 2009; Chavanis & Sire 2005). Observational evidence for the S_q statistics has been found in connection with various astrophysical problems (Livadiotis & McComas 2013; Silva et al. 2013). In particular, the maximization under suitable constraints of the Tsallis entropy of a Newtonian self-gravitating N-body system leads to polytropic distributions (Plastino & Plastino 1993; Lima & de Souza 2005), which can have finite mass and a shape resembling the dark matter (DM) distribution found in numerical simulations of galaxy formation (Navarro et al. 2004; Calvo et al. 2009; An & Zhao 2013; Sánchez Almeida & Trujillo 2021). In addition to having an identifiable physical origin, the polytropic shape has gained practical importance lately because of its association with real self-gravitating astrophysical objects. The mass density profiles in the centers of dwarf galaxies are well reproduced by polytropes without any tuning or degree of freedom (Sánchez Almeida et al. 2020). The same type of profile also explains the stellar surface density profiles observed in globular clusters (I. Trujillo & J. Sanchez Almeida 2021, in preparation). Polytropic DM haloes fit the velocity dispersion observed in many early type galaxies (Saxton & Ferreras 2010).

The second argument against thermodynamic equilibrium explaining galaxy shapes relies on the nature of DM. In the current cosmological model, DM provides most of the gravitational pull needed for the ordinary matter to collapse forming visible galaxies. The simplest form of DM is collisionless and so is unable to reach thermodynamic equilibrium within the Hubble time (Binney & Tremaine 2008; Power et al. 2003; Ludlow et al. 2019), which poses a potential problem for any distribution arising from thermodynamic equilibrium. However, this assumption on the nature of the DM causes some of the so-called small-scale problems of the cold dark matter (CDM) model (e.g., Weinberg et al. 2015; Bullock & Boylan-Kolchin 2017; Del Popolo & Le Delliou 2017), in particular, the cusp–core problem. Simulated CDM haloes have cusps in their central mass distribution (e.g., Navarro et al. 1997; Wang et al. 2020), which disagree with the central plateau or core often observed in galaxies (e.g., Oh et al. 2015). The various physical processes invoked to turn cusps into cores, in essence, produce the thermalization of the overall gravitational potential within a finite timescale. Among others, the proposals have been feedback of the baryons on the DM particles through gravitational forces (Governato et al. 2010; Di Cintio et al. 2014; Freundlich et al. 2020), scattering with massive gas clumps (Elmegreen & Struck 2013; Struck & Elmegreen 2019), forcing by a central bar (Hohl 1971), or assuming an artificially large DM collision cross section that shortens the two-body collision timescales below the Hubble time (Spergel & Steinhardt 2000; Davé et al. 2001; Elbert et al. 2015). When the DM particles of a numerical simulation are allowed to interact through any of these mechanisms, it leads to a gravitational potential conforming with polytropes (Sánchez Almeida & Trujillo 2021), as expected if the maximum Tsallis entropy were dictating the stationary state.

Thus, the existing arguments against thermodynamic equilibrium setting galaxy shapes are questionable. Here, we contribute to the discussion by studying whether polytropes can reproduce the stellar mass distribution in galaxies. Previous works along this direction are restricted to only a few galaxies and yielded inconclusive results (e.g., Chavanis & Sire 2005; Novotný et al. 2021). The question is addressed in various complementary ways. The paper starts off by formally introducing polytropes and their projection in the plane of the sky (Section 2). Then, Section 3 presents our python code to fit surface density profiles using the projection of a polytrope in the plane of the sky (hereinafter denoted as propol). As we explain above, Sérsic profiles are commonly used to model surface density profiles. We use the fitting code to address whether Sérsic profiles and propols are related. We find that the range of observed Sérsic shapes seem to emerge from the range of physically possible polytropes (Section 4). The correspondence holds even when stars do not strictly follow the total mass distribution. In Section 5, the propol fitting code is systematically applied to a large number of galaxies (∼750) with carefully measured mass density profiles, including all morphological types, and spanning five orders of magnitude in stellar mass ($7\lt {\rm{log}}[{M}_{\star }/{{\rm{M}}}_{\odot }]\lt 12$, where M_⋆/M_⊙ stands for the total mass in stars in solar mass units; Trujillo et al. 2020a). The overall good fits provided by propols worsen for massive galaxies. Section 6 is devoted to understanding why through the examination of a few galaxies that are particularly well observed (Spavone et al. 2017). The main results are highlighted and discussed in Section 7. Various technical details are spelled out in Appendices A–D. The fitting code is made publicly available at https://github.com/jorgesanchezalmeida/polytropic-fits.

A polytrope of index m is defined as the spherically symmetric self-gravitating structure resulting from the solution of the Lane–Emden equation for the (normalized) gravitational potential ψ (Chandrasekhar 1967; Binney & Tremaine 2008),

$$\begin{eqnarray}\displaystyle \frac{1}{{s}^{2}}\displaystyle \frac{d}{{ds}}\left({s}^{2}\displaystyle \frac{d\psi }{{ds}}\right)=\left\{\begin{array}{ll}-3{\psi }^{m} & \psi \gt 0,\\ 0 & \psi \leqslant 0.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

The symbol s stands for the scaled radial distance in the 3D space and the mass volume density is recovered from ψ as

$$\begin{eqnarray}&&\rho (r)=\rho (0)\psi {\left(s\right)}^{m},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&r=b\,s,\end{eqnarray} \tag{ 3 }$$

where r stands for the physical radial distance and ρ(0) and b are two arbitrary constants. Equation (1) is solved under the initial conditions ψ(0) = 1 and d ψ(0)/ds = 0. ⁵

A number of properties are relevant for the present work. Polytropes with m < 5 have finite size, i.e., ρ(r) = 0 for r larger than a given radius (which depends on m). There is a preferred range of polytropic indexes,

$$\begin{eqnarray}&&3/2\lt m\leqslant 5,\end{eqnarray} \tag{ 4 }$$

for reasons detailed by Binney & Tremaine (2008) (see also Plastino & Plastino 1993). Polytropes with m > 3/2 are stable according to the Doremus–Felix–Baumann theorem and Antonov second law for the stability of stellar systems. Stellar polytropes with m ≤ 3/2 are unrealistic because, at the cutoff point associated with particles of vanishing energy, the polytropic phase-space distribution has a singularity. Finally, polytropes with m > 5 are unphysical because they have infinite mass and energy. The limit m → ∞ corresponds to the isothermal sphere, which is the solution corresponding to the thermodynamic equilibrium defined by the Boltzmann–Gibbs entropy (Hunter 2001; Binney & Tremaine 2008). Since the isothermal sphere corresponds to m > 5, it has infinite mass and energy. In general, the polytropes have to be evaluated by numerical integration of Equation (1), with two exceptions (e.g., Binney & Tremaine 2008) corresponding to m = 1, with

$$\begin{eqnarray}\psi (s)=\left\{\begin{array}{ll}\displaystyle \frac{\sin (\sqrt{3}\,s)}{\sqrt{3}\,s} & s\lt \pi /\sqrt{3},\\ 0 & \mathrm{elsewhere},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

and m = 5 (Schuster 1884 sphere or Plummer 1911 model),

$$\begin{eqnarray}&&\psi (s)=\displaystyle \frac{1}{\sqrt{1+{s}^{2}}},\end{eqnarray} \tag{ 6 }$$

expressions employed in Section 3 to test our numerical codes.

Note that all polytropes have cores, in the sense that $d\mathrm{ln}\rho /d\mathrm{ln}r\to 0$ when r → 0. This result follows from the initial condition d ψ(0)/ds = 0 and Equations (1)–(3). Moreover, the cores of all polytropes look the same after suitable normalization (Sánchez Almeida & Trujillo 2021),

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\rho (r)}{{\rho }_{\alpha }} & \simeq & {\left[1+\displaystyle \frac{\alpha }{2m}\left(1-\displaystyle \frac{{r}^{2}}{{r}_{\alpha }^{2}}\right)\right]}^{m}\\ & \simeq & 1+\displaystyle \frac{\alpha }{2}\left(1-\displaystyle \frac{{r}^{2}}{{r}_{\alpha }^{2}}\right),\end{array}\end{eqnarray} \tag{ 7 }$$

with r_α defined as

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}r}({r}_{\alpha })=-\alpha ,\end{eqnarray} \tag{ 8 }$$

so that ρ_α = ρ(r_α ). Equation (7) holds for r < r_α and α ≪ 2m.

We will compare the predictions from polytropic mass distributions with observed mass surface densities. Thus, the predicted volume densities have to be projected in the plane of the sky to make the comparison. The surface density Σ(R) corresponding to the volume density ρ(r) is given by its Abel transform,

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=2\,{\int }_{R}^{\infty }\displaystyle \frac{r\,\rho (r){dr}}{\sqrt{{r}^{2}-{R}^{2}}},\end{eqnarray} \tag{ 9 }$$

with R the projected distance from the center (e.g., Binney & Tremaine 2008). Then, Σ(R) can be expressed in terms of the normalized Abel transform

$$\begin{eqnarray}&&f(x,m)=2\,{\int }_{x}^{\infty }\displaystyle \frac{s\,{\psi }^{m}(s){ds}}{\sqrt{{s}^{2}-{x}^{2}}}\end{eqnarray} \tag{ 10 }$$

as

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=a\,f(R/b,m),\end{eqnarray} \tag{ 11 }$$

where x = R/b is the normalized projected distance to the center, and a = ρ(0) × b has units of surface density. The projected polytrope Σ(R) (propol) inherits its shape from f but depends on two parameters encoding the central surface density a and the absolute radial distance b. Equation (11) will be used to fit observed surface density profiles, with a, b, and m the free parameters to be obtained from the fit. The propol shape f can be evaluated analytically when m = 5 (Equation (6)), and it corresponds to

$$\begin{eqnarray}&&f(x,5)=\displaystyle \frac{4/3}{{\left(1+{x}^{2}\right)}^{2}},\end{eqnarray} \tag{ 12 }$$

with its half-mass–radius at x = 1.

Comparing propols with observed stellar surface densities assumes the stellar mass to scale with the total mass responsible for the overall gravitational potential. The stellar mass results from replacing ρ in Equation (9) with the stellar density ρ_⋆. If ρ_⋆ ∝ ρ, then the shape Σ(R) in Equation (11) is the same for ρ and ρ_⋆. However, ρ_⋆ is not proportional to ρ in real galaxies. According to the current cosmological model, most of the gravitation is provided by DM, which is more spread-out than the centrally concentrated stars responsible for the visible light. This change can be easily accommodated within the above formalism. Assuming a moderate change in the stellar to total mass ratio, then (Appendix A)

$$\begin{eqnarray}&&{\rho }_{\star }(r)\simeq {\rho }_{\star }(r^{\prime} ){\left[\displaystyle \frac{\rho (r)}{\rho (r^{\prime} )}\right]}^{1+\gamma },\end{eqnarray} \tag{ 13 }$$

with

$$\begin{eqnarray*}&&\gamma =\displaystyle \frac{d\mathrm{ln}({\rho }_{\star }/\rho )}{d\mathrm{ln}\rho },\end{eqnarray*}$$

and r and $r^{\prime}$ are any two radii. In galaxies, both ρ_⋆/ρ and ρ drop outward (e.g., Kravtsov 2013; Iocco et al. 2015) leading to γ > 0. An educated guess (Appendix A) renders γ ∼ 0.3, which corresponds to a change of one order of magnitude in ρ_⋆/ρ when ρ changes by three orders of magnitude. Thus, considering the difference between stellar mass and total mass is just a matter of replacing f(x, m) in Equation (11) with g(x, m, γ),

$$\begin{eqnarray}&&g(x,m,\gamma )=2\,{\int }_{x}^{\infty }\displaystyle \frac{s\,{\left[\psi (s)\right]}^{(1+\gamma )m}\,{ds}}{\sqrt{{s}^{2}-{x}^{2}}}.\end{eqnarray} \tag{ 14 }$$

In order to illustrate the range of shapes described by the propols, Figure 1 shows f and g for m from 1–9 when γ = 0 (a) and γ = 0.3 (b). As it happens with the polytropes, all propols have a core (a central plateau) and their size is finite when m < 5 (see Figure 1(a)). As expected, the profiles drop with projected radius much faster when γ ≠ 0 (compare panels (a) and (b) in Figure 1). The characteristic length scales corresponding to theses polytropes are shown in Figure 2. We include the effective radius R_e (i.e., the half-mass–radius of the propol) as well as the core radius R_c , defined as the radius where f (or g) drops to 90% of its maximum value. Given m, both R_e and R_c scale linearly with b. Indexes m > 5 are not shown since these propols have infinite mass and their R_e is undefined (actually, it is infinite). Note that R_e is smaller but of the order of b, and that the ratio R_e /R_c spans from 2–4 for m from 1–5. The dashed lines in Figure 2 correspond to varying stellar to total mass ratio with γ = 0.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Fitting mass surface densities with propols requires evaluating f(x, m) (or g[x, m, γ]) many times. Since there is no general analytic expression for the polytropes, neither is there for f (or g). In order to speed up the fitting procedure, the evaluation of f was carried out by interpolating on a precomputed grid f(x_i , m_i ). We precompute several grids suited to the particular problem to be treated. We include grids covering different ranges of m, different sampling intervals, and different values of γ. Their main characteristics are listed in Table 1. The calculation of the grids and the assessment of uncertainty are described in Appendix B. From these tests, we know that the propol shape, f, is evaluated with a relative precision better than 10⁻⁶, which suffices for most practical applications since other ubiquitous sources of error (from interpolation to observational uncertainties) are larger.

We use a least squares approach to carry out the fits. The final least squares algorithm minimizes the merit function,

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{i}{\left[{y}_{i}-\mathrm{log}\,a-\mathrm{log}f({R}_{i}/b,m)\right]}^{2},\end{eqnarray} \tag{ 15 }$$

where y_i is the log surface density observed at R = R_i , and the subscript i varies from 1 to the number of observed radii in the profile. The minimization of the merit function is carried out using the classical Levenberg–Marquardt minimization algorithm (Press et al. 1986) as implemented in python (scipy.least_squares). The free parameters to be determined are m, $\mathrm{log}a$, and b, and their formal errors are estimated from the diagonal of the covariance matrix. While fitting, the propol are computed by 2D interpolation on the selected grid (Table 1). We use a linear 2D interpolation as provided by interp2d in python (scipy.interpolate).

As we point out in Section 1, the mass within galaxies drops with radial distance following a law approximately given by Sérsic functions (Sersic 1968; Caon et al. 1993; Graham & Driver 2005),

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={\rm{\Sigma }}(0)\exp [-{c}_{n}{\left(R/{R}_{e}\right)}^{1/n}],\end{eqnarray} \tag{ 16 }$$

with Σ(R) the mass surface density at a distance R from the center, R_e the radius enclosing half of the mass, and c_n a constant that only depends on the parameter n, called Sérsic index. The Sérsic index n controls the shape of the profile, and has been observed to vary mostly from 0.5 to 6 (Blanton et al. 2003; van der Wel et al. 2012), going from disk-like galaxies to elliptical galaxies, and from low-mass dwarfs to the brightest galaxy in a cluster. As we also explain in Section 1, the range of shapes includes exponential disks (n = 1; de Jong & van der Kruit 1994) and the de Vaucouleurs R^1/4 law characteristic of massive ellipticals (n = 4; de Vaucouleurs 1948).

The question arises as to whether the empirical Sérsic profiles and the theoretical propols are equivalent. All propols have cores, inherited from their parent polytrope (see Figure 1). Sérsic profiles also present cores, but of different nature. Except when n ≲ 1, they are too small to be comparable to the cores of the propols. This mismatch in the cores can be shown analytically: the logarithmic derivative of the surface density given by a Sérsic function (Equation (16)) can be approximated as (e.g., Graham & Driver 2005),

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{\rm{\Sigma }}(R)}{d\mathrm{ln}R}\simeq -(2-0.33/n)\cdot {\left(R/{R}_{e}\right)}^{1/n},\end{eqnarray} \tag{ 17 }$$

which tends to zero when R → 0, but very slowly when n ≳ 1 because of the 1/n exponent. Defining the core radius as the radius where the logarithmic derivative of the profile is small enough (0.3 is often used; e.g., Oh et al. 2015), then the resulting core radius R_0.3, implicitly defined as

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{\rm{\Sigma }}({R}_{0.3})}{d\mathrm{ln}R}\simeq -0.3,\end{eqnarray} \tag{ 18 }$$

turns out to be comparable with the effective radius for n = 0.5 (R_0.3/R_e ≃ 0.5) but becomes tiny for n = 5 (≃10⁻⁴). Therefore, except for n ≲ 1, the propols cannot reproduce the Sérsic profiles for R ≲ R_e .

As it has been evidenced elsewhere (e.g., Trujillo et al. 2004) and we will show in Section 5, rather than a problem for the propols, this mismatch often reveals the disability for the Sérsic profiles to reproduce some of the central cores observed in galaxies. Thus, to study the equivalence between Sérsic profiles and propols, we will fit Sérsic profiles with propols using the code presented in Section 3, but avoiding the central cores when n ≳ 1. Figure 3 shows examples of such fits. The range of fitted radii is indicated in the figure using symbols. Even if the cores are excluded, the fits still comprise a factor of 20 in radii and up to a factor of 10⁵ in surface density. Within this range, the agreement is well within any realistic observational error (see the residuals shown index in the secondary panels). Figure 3 only shows three representative Sérsic profiles. The results for the whole range of observed Sérsic indexes are summarized in Figure 4. The blue symbols and lines represent the fits carried out using grid2 (Table 1) and avoiding the cores as shown in Figure 3. We note that the agreement between Sérsic profiles and propols is excellent, with the rms of the residuals well below any realistic observational error: around 0.02 dex or 5% (Figure 4(d)) over one order of magnitude in radius (Figure 4(b)) and five orders of magnitude in surface density (Figure 4(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The correspondence between the original Sérsic index and the polytropic index of the fitted propol is shown in Figure 4(a). Interestingly, the range of observed Sérsic indexes, 0.5 ≲ n ≲ 6.0 (e.g., Trujillo et al. 2001a; Blanton et al. 2003; Tasca & White 2011; van der Wel et al. 2012), gives rise to a range of polytropic indexes between 2 and 5 (the blue symbols and lines), right in the range where polytropes are physically sensible (Section 2, Equation (4)). In other words, the range of allowed polytropic indexes (Equation (4)) corresponds to Sérsic indexes in the observed range, a result whose implications are hard to judge at this point because, when n ≳ 2, the strict equivalence between Sérsic profiles and propols breaks down when the cores are included.

In order to study the dependence of the above results on details of the fit, we repeat the same exercise varying some of the hypotheses and hyper-parameters that define the fits. Thus, if rather than completely excluding the cores, we partly include them, then the fits worsen and render polytropic indexes m > 5 for n ≳ 3 (see the green lines and symbols in Figure 4(a), and Appendix C for examples of fits and residuals). The increase of rms is significant (from 0.02–0.05; see the green and blue lines in Figure 4(d)). If rather than using pure propols, we assume ρ_⋆/ρ to vary radially (γ ≠ 0 in Equation (13)), then the fits also get worse: see Figure 4(d), the orange line as well as the examples in Appendix C. In this case, grid2g3 in Table 1 was used for fitting. As it happens when the cores are partly included, m becomes >5 when n ≳ 3 (Figure 4(a)); examples are given in Appendix C. When a different more coarse propol grid is used for fitting (e.g., grid3 in Table 1) the results are indistinguishable from those shown as blue lines and symbols in Figure 4, which are based on grid2. Finally, rather than fitting propols to Sérsic profiles, we also tried fitting Sérsic profiles to propols. When the range of radii are similar to the range used when fitting propols to Sérsic profiles, then the fits are equally good. The red line and symbols in Figure 4(a) correspond to this case. Summing up all these results, we conclude that the good correspondence between the outskirts of Sérsic profiles and propols remains valid in quite general terms. Outside the optimal range of radii and γ, the fits worsen and the best-fitting propols start having m > 5, which corresponds to polytropes of unbound mass and energy (Section 2). These trends also appear when fitting mass profiles in real galaxies (Sections 5 and 6) and in globular clusters (I. Trujillo & J. Sanchez Almeida 2021, in preparation).

Differences between Sérsic profiles and propols are observed in the center (R ≲ R_e ) and outskirts (R ≫ R_e ), but this mismatch with Sérsic profiles is also shared with real galaxies. They also deviate from pure Sérsic functions because of the presence of cores in their centers (e.g., Trujillo et al. 2004; Del Popolo & Le Delliou 2017), and because of the existence of truncations and stellar haloes in their outskirts (e.g., Pohlen & Trujillo 2006; Erwin et al. 2008). Thus, we try direct fits of observed mass density profiles in galaxies to see whether propols provide a fair description of the observation. As an exploratory work, we do not intend to reproduce the profiles in detail, rather, we want to assess the overall goodness of the fits and if they are comparable in quality to the traditional fits based on Sérsic profiles.

5.1. Description of the Data Set

We select the data set described in detail by Trujillo et al. (2020a, 2020b) to study how propols fit real galaxies. It comprises around 750 galaxies of the local Universe (redshift < 0.09), including all morphological types (as classified by Nair & Abraham 2010; Maraston et al. 2013), and spanning five orders of magnitude in stellar mass ($7\lt {\rm{log}}[{M}_{\star }/{{\rm{M}}}_{\odot }]\lt 12;$ see Figure 5). The sample excludes objects with signs of strong interactions (e.g., having shells or tails; Rampazzo et al. 2020) to keep only galaxies with relaxed outskirts. The profiles were derived from images of the IAC Stripe 82 Legacy Project (Fliri & Trujillo 2016; Román & Trujillo 2018), which is a re-processing of the SDSS Stripe 82 (Frieman et al. 2008) with improved sky subtraction and reaching a surface brightness limit. ⁶ of 29.1, 28.6, and 28.1 mag arcsec² in the bands g, r and i, respectively. The scattered light from stars in the field was removed, and the surface brightness was corrected for galaxy inclination as described by Trujillo et al. (2020a). The point-spread function (PSF), as measured from stars in the field, has a central Gaussian core with FWHM ≃ 1.4''. Stellar mass surface densities were computed from surface brightness using mass-to-light ratios inferred from the observed colors, as described by Roediger & Courteau (2015). The underlying hypotheses of this calibration include the use of Bruzual & Charlot (2003) models for the stellar populations and a Chabrier (2003) initial mass function. All in all, the procedure provides mass profiles for all galaxies down to a surface density fainter than 1 M_⊙ pc⁻². Examples of profiles are shown in Figures 6–7 (the blue symbols).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. Fitting Surface Density Profiles with Propols

The profiles described in Section 5.1 were fitted with the python code introduced in Section 3. In order to minimize observational biases, we did not include radii smaller than 0.7'' and g-band surface brightness fainter than 29 mag arcsec⁻². We also excluded profiles with less than four points after thresholding. The first constraint avoids radii that may be influenced by the PSF whereas the second discards points at or below the noise level. Apart from those discarded, all radii contribute equally to the goodness of the fit (Equation (15)). We take this approach since the residuals of the fits are mainly produced by unknown systematic errors, arising from the lack of symmetry in the galaxies, the correction of inclination, the assumed mass-to-light ratio and, obviously, the deviations of the true galaxy profiles from the propol shape. A histogram with the stellar masses of the selected galaxies is shown in Figure 5—this sample will be called reference to be distinguished from other cutouts from the original sample where some of the restriction are lifted or modified. The reference sample uses grid2 for fitting (Table 1).

Examples of good fits are shown in Figure 6. Examples of fits that are not so good are included in Figure 7. To decide whether the fits are good or not, we use as benchmark Sérsic profile fits to the same data. Good corresponds to propol fit with rms smaller than the Sérsic fit, and vice versa. Thus, together with the observed profiles (blue solid symbols) Figures 6 and 7 include the propol fit (orange solid lines) and the Sérsic profile fit (green dashed lines). One can find good and poor fits in all the range $7.5\leqslant {\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\leqslant 11.5$, although there is a clear trend of the propol fits to worsen with increasing stellar mass, as we report on below. Note that the propols with m < 5 have finite size (see Figure 1). This size is always larger than the last point used for fitting, but it is often smaller than the radii of some of the excluded points (the open symbols in Figures 6 and 7).

Zoom In Zoom Out Reset image size

Diagnostic diagrams comparing the Sérsic and propol fits are included in Figure 8, with a symbol for each galaxy color coded according to $\mathrm{log}{M}_{\star }$. Figure 8(a) shows Sérsic index versus polytropic index whereas Figure 8(b) shows rms of the Sérsic fit residuals versus rms of the propol fit. The dashed line in Figure 8(a) corresponds to the equivalence between Sérsic index and polytropic index worked out in Section 4 and represented by a blue line with symbols in Figure 4(a). The dashed line in Figure 8(b) corresponds to y = x so points above (below) the line correspond to propol fits better (worse) than the corresponding Sérsic fit. A number of conclusions can be drawn from these figures as well as similar figures for different galaxy selections and grids included in Appendix D to avoid cluttering. These conclusions are:

• 1.  
  The relative goodness of the propol fits with respect to the Sérsic fits depends on the galaxy stellar mass, with the goodness of the propol fits increasing with decreasing stellar mass. Thus, propol fits are systematically better than Sérsic fits for ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\lesssim 9$ and systematically worse for ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\gtrsim 10$. This effect is clear when comparing Figure 8(b) with the same figure when only low-mass galaxies are plotted (${\rm{log}}[{M}_{\star }/{{\rm{M}}}_{\odot }]\lt 10;$ Figure 13(b)).
• 2.  
  Even if the worsening of the fits with increasing stellar mass is clear, there are high-mass galaxies having good propol fits and low-mass galaxies where Sérsic fits do a better job (Figures 6 and 8(b)).
• 3.  
  The relation between the Sérsic index and the Polytropic index only follows the theoretical relation worked out in Section 4 when m ≲ 4 (see the dashed line in Figure 8(a)). The deviations at larger indexes can be understood as the effect of including galaxy cores when fitting. The theoretical line was computed to minimize the rms of the residuals, which forces excluding cores when m ≳ 1 (see Figure 3). To show that cores are producing the deviation, we repeat the fits excluding most of the core when fitting (excluding the central R < 2''). The resulting relation is shown in Figure 16(a), where the points corresponding to the most massive objects have shifted toward the theoretical curve. Simultaneously, the systematic difference between the rms of the propol and Sérsic fits at high mass tends to go away (Figure 16(a)).
• 4.  
  Including the expected difference between the DM and the stellar mass density profiles does not improve the quality of the fits. Figure 17 shows another rendering of the diagnostic plot using grid2g1 rather than grid2, where the parameter that characterizes the difference (γ in Equation (13)) is set to 0.1 rather than 0 (Table 1).
• 5.  
  Fits based on other samplings of the propol grid and extending the range of indexes (e.g., grid3 in Table 1) do not modify the above conclusions in any significant way.
• 6.  
  Polytropes represent spherically symmetric structures, however, deviations from this underlying hypothesis are to be expected in real galaxies. Since early type galaxies are rounder than late type galaxies, they were anticipated to have better propol fits than late types. However, we do not see any obvious systematic worsening of the fits associated with late type objects (compare Figures 14 and 15). Differences can be easily ascribed to the fact that, in our sample, early types tend to be more massive (Figure 5). We do not have a clear physical interpretation for this negative result. It may be due to the fact that both early and late types are massive enough for the fitting errors to hide the putative differences between the two populations.

Zoom In Zoom Out Reset image size

The discussion above is made in terms of goodness relative to Sérsic fits, but we also quantify the fraction of good fits in absolute terms. In this context, we define as good any fit having rms <0.1, which is the limit for the good-fit profiles shown in Figure 6. As Figure 8(b) evidences, this fraction very much depends on M_⋆. Thus, we get 84% of good fits when ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\lt 8$, 81% when ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })$ is in between 8 and 9, 56% when in between 9 and 10, 34% when in between 10 and 11, and 23% when larger than 11. The percentage of good fits when all galaxies are included is 49%.

The propol fits tend to crowd around the index m = 5, with a dispersion that increases as the stellar mass decreases (Figure 9). As we explain in Section 2, the profiles having m > 5 have infinite mass, which may pose a problem of interpretation. The problem can be sorted out assuming that the existence of significant deviations from the propol profiles outside the fitted range, which is not unexpected after all. Breaks in the mass profiles cannot be reproduced by propols, although they are quite common (e.g., Pohlen & Trujillo 2006). Moreover, even when they are not present and the observed profiles vary smoothly, the outer parts of some galaxies may have not reached thermodynamic equilibrium yet (Section 6).

Zoom In Zoom Out Reset image size

Massive early type galaxies are often poorly reproduced by propols, as we uncover in Section 5. Here we try to understand why. We analyze profiles of six massive elliptical galaxies observed within the VST Early type GAlaxy project by Spavone et al. (2017). They were chosen because their profiles represent state-of-the-art observations of very large objects, so that the potential systematic effects arising from low noise or insufficient spatial sampling are minimized. Rather than mass profiles, the authors provide surface brightness profiles. Assuming a constant mass-to-light ratio along the galaxy, the mass and light surface density profiles are identical except for a scaling factor. Thus, we fit light profiles assuming them to present mass profiles, specifically, we fit the surface brightness in the g band, μ_g , provided by Spavone et al. (2017). (Trials with mass profiles are discussed later on, but they do not significantly improve or worsen the goodness of the fit.) The example of one of such fits (NGC 1399) is shown in Figure 10. Fits to the other massive ellipticals (NGC 3923, NGC 4365, NGC 4472, NGC 5044, and NGC 5846) are similar and will not be shown. For reference, the stellar masses of these galaxies are in the range between 4.3 × 10¹¹ and 1.5 × 10¹² M_⊙ (Spavone et al. 2017).

Zoom In Zoom Out Reset image size

The profiles are well defined for radii in between 1'' and 5'. Radii smaller than 1'' are discarded to avoid seeing smearing the profile, while radii larger than 5' are affected by breaks and truncations in the profile. The propol fit to NGC 1399 is shown in Figure 10(a). There are clear systematic deviations between the best-fitting propol and the observed profiles. The Sérsic profile included in the same plot clearly do a better job, although the deviation at the core starts to be substantial. The observed profile shows a central core, but the range of shapes provided by the propols cannot fit the core and the outskirts of the profile simultaneously. When the outskirts are not included in the fit, then propols make an excellent job fitting the core, as shown in Figure 10(b). However, the propol fit to only the outskirts is unrealistically far from the observed profile (Figure 10(c)).

Thus, propols only provide a good fit to the cores of these massive ellipticals. The full profile can be represented as an inner polytropic core plus a Sérsic profile describing the outskirts (Figure 10(d)), very much in line with the piece-wise profiles describing the haloes from self-interacting DM numerical simulations (Robertson et al. 2021; Sánchez Almeida & Trujillo 2021). This fact suggests that the inner part of these massive ellipticals is in thermal equilibrium, which does not reach the outer part. Thus, the physical process leading to thermalization should be more effective in the high-density environment corresponding to the inner part. The required scaling with density is quite natural, e.g., it appears when the equilibrium is reached via two-body gravitational interactions or by other types of collisions between DM particles (e.g., Sánchez Almeida & Trujillo 2021). Another pathway to thermalization may be the merger of supermassive black holes (SMBHs), expected to occur at the center of massive ellipticals (López-Cruz et al. 2014; Mazzalay et al. 2016). The motion of two merging black holes produces scouring of stars. In addition, the recoil kicks the merged SMBH out of the center, forcing a final swing of the SMBH that stirs the global gravitational potential (e.g., Merritt 2006; Nasim et al. 2021). Thus, scouring plus recoil may allow the self-gravitating system to reach thermodynamic equilibrium in a timescale much shorter than the two-body collision timescale, and well within the Hubble time.

As an additional sanity check, we also transform the surface brightness profiles μ_g into mass surface densities Σ_⋆ using mass-to-light ratios precomputed by Roediger & Courteau (2015) and the radial profile in the i band provided Spavone et al. (2017). The mass profiles differ very little from the light profiles described above. In particular, the difficulty for the propols to provide a good overall fit remains.

The solutions of Equations (1)–(3) are called polytropes. Given a total mass and energy, they describe spherically symmetric self-gravitating systems that maximize Tsallis entropy. These structures, which have been studied by astronomers for long, are nowadays gaining renewed interest because they provide valid explanation for several seemingly disconnected properties of galaxies (Section 1). In this paper, we carry on exploring the practical interest of polytropes by studying if they account for the stellar mass distribution observed in real galaxies.

We develop a python code (Section 3) to fit surface density profiles using polytropes projected in the plane of the sky (propols). The code is systematically applied to understand if and why propols have the shapes observed in the mass distribution of galaxies. Sérsic profiles are generally regarded as good approximations to galaxy shapes. Thus, we first use the fitting code to address whether Sérsic profiles and propols are related (Section 4). We find that the cores of the propols and Sérsic profiles are inconsistent when the Sérsic index is larger than around 1. However, when the cores are excluded, propols and Sérsic profiles are indistinguishable within observational errors. The rms of the residuals is around 0.02 dex or 5% (Figure 4(d)) over one order of magnitude in radius (Figure 4(a)) and five orders of magnitude in surface density (Figure 4(c)). We find a one-to-one correspondence between the Sérsic index and the polytropic index of the best-fitting propol (Figure 4(a), the blue line). Polytropes are physically meaningful only within a particular range of indexes (Equation (4) which corresponds to Sérsic indexes between 0.4 and 6.0, in excellent agreement with the Sérsic indexes observed in nature (e.g., Blanton et al. 2003; van der Wel et al. 2012; Trujillo et al. 2001a).

The propol fitting code has been systematically applied to a large number of galaxies (∼750) with carefully measured mass density profiles (Section 5). The galaxy sample includes all morphological types and spans five orders of magnitude in stellar mass ($7\lt {\rm{log}}[{M}_{\star }/{{\rm{M}}}_{\odot }]\lt 12;$ Trujillo et al. 2020a). The goodness of each fit is established by comparing the rms of its residuals with those resulting from a Sérsic fit. The relative goodness thus defined depends on the galaxy stellar mass, with the goodness of the propol fits increasing with decreasing stellar mass. The propol fits are systematically better than Sérsic s for ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\lesssim 9$ and systematically worse for $\mathrm{log}({M}_{\star }/{\text{}}{{\rm{M}}}_{\odot })\gtrsim 10$. Even if the worsening of the fits with increasing stellar mass is clear, there are high-mass galaxies having good propol fits and low-mass galaxies where Sérsic fits do a better job than propols (Figures 6 and 8(b)). The clean one-to-one relation between Sérsic index and Polytropic index mentioned above breaks down when the cores are included in the fits, as it usually the case when fitting real galaxies, which often have cores. The observed polytropic indexes tend to cluster around m = 5, although with large scatter (Figure 9). We also study whether the expected difference between the total mass and the stellar mass density profiles modifies the quality of the fits, to find out that it does not. Using a threshold in the rms of the residuals to identify good fits in absolute terms (<0.1 dex in our case), we characterize the fraction of good fits depending on M_⋆. It goes from 84% for ${\rm{log}}({M}_{\star }/{{\rm{M}}}_{\odot })\lt 8$ to 23% when it is larger than 11. When all galaxies are included, the percentage reaches 49%. Section 6 is devoted to understand why the overall good fits provided by propols worsen for massive galaxies. We find that propols are very good at reproducing the central parts, but they do not handle well cores and outskirts all together. However, the combination of a propol core and a Sérsic outskirt matches the observed profiles extremely well (Figure 10(d)). This fact suggests that only the inner part of a massive galaxy is already in thermal equilibrium, which is sensible considering that the physical process leading to thermalization are more effective in high density environments. The merger of SMBHs at the center of massive ellipticals also provides a pathway to thermalization (Section 6).

There is no systematic differences between the goodness of the fits in early and late type galaxies, even though the hypothesis of being spherically symmetric is better satisfied by early types. We do not have a proper interpretation of this negative result. It may be due to the fact that both early and late types are massive enough for the fitting errors to hide the putative differences between the two populations.

On the whole, propols do an excellent job reproducing the stellar mass distribution of low-mass galaxies, as well as the cores of the massive objects. We have also shown that the outskirts of propols are indistinguishable from Sérsic profiles. In addition, propols should also reproduced the tight mass – size relation observed in galaxies when measuring sizes at constant surface brightness (e.g., Trujillo et al. 2020a), because the tight relation is naturally produced by Sérsic profiles (Sánchez Almeida 2020), and so should result from propols as well. Finally, polytropes reproduce the cores in the total mass distribution of dwarf galaxies without being adjusted (Sánchez Almeida & Trujillo 2021).

Rather than being purely empirical constructs, polytropic shapes are expected when self-gravitating structures are in thermal equilibrium as defined by the Tsallis entropy (Plastino & Plastino 1993, see Section 1). Thus, our results are consistent with the ansatz that the principle of maximum Tsallis entropy dictates the internal structure in dwarfs and in the central region of the most massive galaxies. Should this ansatz be correct, it has profound implications ranging from measuring physical properties of galaxies (should the Navarro–Frenk–White profiles (Navarro et al. 1997) be replaced with polytropes when interpreting gravitational lensing signals?) to the basic physics underlying galaxy formation (how is the thermodynamic equilibrium established? Is it due to some sort of nongravitational interaction between DM particles?), and including galaxy modeling through numerical simulations (why do cold DM numerical simulations create cusps rather than cores? See Sánchez Almeida & Trujillo 2021).

We thank Nuskia Chamba for help with the handling of the profiles analyzed in Section 5, and Marilena Spavone and Enrichetta Iodice for providing the surface brightness profiles analyzed in Section 6. Thanks are also due to an anonymous referee for suggestions leading to improve the presentation of the results and some arguments. J.S.A. acknowledges support from the Spanish Ministry of Science and Innovation, project No. PID2019-107408GB-C43 (ESTALLIDOS), and from Gobierno de Canarias through EU FEDER funding, project No. PID2020010050. I.T. acknowledges support from the Spanish Ministry of Science and Innovation, grant No. PID2019-107427GB-C32, from the European Union Horizon 2020 research and innovation program under Marie Skłodowska-Curie grant agreement No. 721463 to the SUNDIAL ITN network, and from the European Regional Development Fund (FEDER), IAC project No. P/300624. The IAC project is partly financed by the Spanish Ministry of Science, Innovation and Universities (MCIU), through the State Budget, and by the Canary Islands Department of Economy, Knowledge and Employment, through the Regional Budget of the Autonomous Community.

Facilities: SDSS (York et al. 2000), IAC Stripe 82 (Fliri & Trujillo 2016).

Software: astropy (Astropy Collaboration et al. 2013).

The Lane–Emden equation describes the distribution of total gravitational mass, whereas observations refer to stellar mass. The two masses do not necessarily scale each other, in which case Equation (11) would be useless for stellar mass. Fortunately, differences between the total density ρ and the stellar density ρ_⋆ can be easily accommodated within our formalism. Assuming a smooth and moderate change in the stellar to total mass ratio, then,

$$\begin{eqnarray}&&\mathrm{ln}[\displaystyle \frac{{\rho }_{\star }}{\rho }(r)]\simeq \mathrm{ln}[\displaystyle \frac{{\rho }_{\star }}{\rho }(0)]+\gamma \,[\mathrm{ln}\rho (r)-\mathrm{ln}\rho (0)]+...,\end{eqnarray} \tag{ A1 }$$

with

$$\begin{eqnarray*}&&\gamma =\displaystyle \frac{d\mathrm{ln}[\tfrac{{\rho }_{\star }}{\rho }]}{d\mathrm{ln}\rho }(0),\end{eqnarray*}$$

so that one can approximate ρ_⋆ as

$$\begin{eqnarray}&&{\rho }_{\star }(r)\simeq {\rho }_{\star }(0){\left[\displaystyle \frac{\rho (r)}{\rho (0)}\right]}^{1+\gamma },\end{eqnarray} \tag{ A2 }$$

or, in general,

$$\begin{eqnarray}&&{\rho }_{\star }(r)\simeq {\rho }_{\star }(r^{\prime} ){\left[\displaystyle \frac{\rho (r)}{\rho (r^{\prime} )}\right]}^{1+\gamma }.\end{eqnarray} \tag{ A3 }$$

Therefore, considering the variation of ρ_⋆/ρ with r is equivalent to deriving ρ from the Lane–Emden Equation (1), computing ρ_⋆ through Equation (A2), which is then projected on the plane of the sky replacing ρ with ρ_⋆ in Equation (9). We want to stress that the above equations express a purely phenomenological relation between ρ(r) and ρ_⋆(r). It does not imply any interpretation of the nature of this relation of the kind of, e.g., the radial acceleration relation (McGaugh et al. 2016). The densities ρ(r) and ρ_⋆(r) are both monotonic functions of r, therefore, their relative variation with r can always be expressed as the Taylor series in Equation (A1). If the ratio of densities changes moderately, then the first two terms of the expansion suffice, leading to Equation (A2).

This approach requires a value for γ, which can be estimated as follows: γ > 0 since the stars are more centrally concentrated than the DM and so both ρ_⋆/ρ and ρ decrease outward in a galaxy. For a significant (one order of magnitude) change from [ρ_⋆(0)/ρ(0)] to [ρ_⋆(r)/ρ(r)], Equation (A2) requires

$$\begin{eqnarray}&&{\left[\rho (r)/\rho (0)\right]}^{\gamma }\leqslant \displaystyle \frac{1}{10}.\end{eqnarray} \tag{ A4 }$$

The density ρ changes by many orders of magnitude in any galaxy. Thus, for a typical ρ(r)/ρ(0) ∼ 10⁻³, Equation (A4) is met for γ = 0.3. If the radial drop is larger, then a smaller γ suffices.

There are two critical points in the evaluation of the grid that defines the shapes of the propols (i.e., f). The first one has to do with the integration of the Lane–Emden Equation (1). We split it into a system of two first order differential equations for ψ and d ψ/ds, which were integrated from s = 0 with the initial conditions ψ(0) = 1 and d ψ/ds(0) = 0 using Lsoda (Hindmarsh 2019) as implemented in python (scipy.odeint). The second one has to do with carrying out the Abel transform of the polytropes (Equation (10)). It was evaluated by direct integration in the 3D radial distance using Simpson's rule (simps provided by scipy.integrate). Each x was integrated independently, with a mesh from 0.01 to 200 equi-spaced in $\mathrm{log}x$ steps of 0.02. The integration parameters were tuned so as to achieve a relative error in the inferred f smaller than 10⁻⁶. The final errors in ψ and f were evaluated resorting to the cases where the Lane–Emden equation has analytic solutions (Equations (5), (6), and (12)). The difference between the resulting numerical and exact solutions are shown in Figure 11. Relative errors in ψ are of the order of 10⁻¹⁰, whereas they are smaller than 10⁻⁶ for f. The actual figure corresponds to grid2 (Table 1).

Zoom In Zoom Out Reset image size

This appendix includes examples of propol fits to Sérsic profiles. The range of hyper-parameters characterizing the fits are varied with respect to the nominal values shown in Figure 3. The left column in Figure 12 corresponds to fits where the cores are partly included, and have to be associated with the green lines and symbols in Figure 4. The right column in Figure 12 shows fits where the stellar to total mass density ratio is not constant with radius. They correspond to the orange lines and symbols in Figure 4.

Zoom In Zoom Out Reset image size

This appendix presents a number of plots that show how the goodness of the propol fits to the observed mass profiles changes when the hyper-parameters of the fit and the sample selection are modified. These plots are referred to and used in Section 5.2. All figures are identical to Figure 8 except that high-mass galaxies are excluded in Figure 13, only late morphological types are included in Figure 14, only early morphological types are included in Figure 15, central cores are not used for fitting in Figure 16, and differences between total mass and stellar mass profiles are considered in Figure 17.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size