Self-gravitating Equilibria of Non-minimally Coupled Dark Matter Halos

N-body, dark matter (DM)-only simulations in the standard ΛCDM cosmology suggest an almost universal shape of the density distributions ρ(r) within virialized DM halos of different masses and redshifts. This fact has been established since the seminal work by Navarro et al. (1996), who proposed an empirical fitting formula (the so-called NFW profile) for the DM density run $\rho (r)\propto \tfrac{1}{r{\left(r+{r}_{s}\right)}^{2}}$ that still constitutes the standard lore to this day (but see Navarro et al. 2010 for refinements); inward of the scale radius r_s a density cusp, i.e., a steep central divergence, is found.

However, such a behavior is at variance with the observational evidences inferred from well-measured rotation curves in many dwarf, DM-dominated galaxies, which point toward a constant, finite inner density ρ₀ within a core radius r₀ (see, e.g., McGaugh et al. 2001; Gentile et al. 2004; de Blok et al. 2008; Walker & Penarrubia 2011; Weinberg et al. 2015; Genzel et al. 2020; for a review and further references, see Bullock & Boylan-Kolchin 2017). The observed density distribution is usually described with the phenomenological Burkert (1995) profile $\rho (r)\propto \tfrac{1}{(r+{r}_{0})({r}^{2}+{r}_{0}^{2})}$.

A nontrivial point to stress is that, besides the flat shape of the inner density run, the measurements indicate a nearly universal relationship between the core density and radius, in terms of a constant value for the "core surface density" ρ₀ × r₀ among different galaxies (see Salucci & Burkert 2000; Donato et al. 2009; Gentile et al. 2009; Burkert 2015; Kormendy & Freeman 2016). In addition, another puzzling aspect of the DM phenomenology is the existence of tight scaling laws with baryonic quantities. These include the baryonic Tully–Fisher relation (Tully & Fisher 1977; McGaugh 2012), the core radius versus disk scale-length relation (see Donato et al. 2004), the radial acceleration relation (see McGaugh et al. 2016), and the universal DM–baryon constant (see Chan 2019).

Admittedly, through the last 20 years the cusp–core controversy has provoked scientific debate. One class of solutions invokes physical processes that can cause violent fluctuations in the inner gravitational potential and/or transfer of energy and angular momentum from the baryons to DM, thus possibly erasing the central density cusp; some of the most explored possibilities include dynamical friction (see El-Zant et al. 2001, 2016; Tonini et al. 2006; Romano-Díaz et al. 2008), and feedback effects from stars and active galactic nuclei (see Governato et al. 2012; Teyssier et al. 2013; Pontzen & Governato 2014; Peirani et al. 2017; Freundlich et al. 2020a). In particular, hydrodynamic simulations including baryonic physics have shown that a variety of halo responses are created, depending on the stellar/halo masses (e.g., Freundlich et al. 2020b).

An alternative solution, perhaps more fascinating, is to abandon the cold DM hypothesis and look at nonstandard particle candidates (see review by Salucci 2019). A few examples include: self-interacting DM particles with cross section ∼1 cm² g⁻¹, which can create a core as they are heated up via elastic two-body collisions and evacuated from the inner region (e.g., Spergel & Steinhardt 2000; Vogelsberger et al. 2014); a self-interacting or non-interacting (alias fuzzy DM) Bose–Einstein condensate of ultralight particles (likely axions) with masses ∼10⁻²² eV, for which the core stems from the equilibrium between quantum pressure and gravity (e.g., Hu et al. 2000; Bohmer & Harko 2007; Schive et al. 2014a, 2014b; Harko 2014; Hui et al. 2017; Bernal et al. 2018); warm DM, made of fermionic particles with masses of order a few keV (e.g., sterile neutrinos), which can originate cores where gravity is counterbalanced by quantum degeneracy pressure from the Pauli exclusion principle (e.g., Dodelson & Widrow 1994; Shi & Fuller 1999; Kusenko 2009; Destri et al. 2013; Adhikari et al. 2017). A more radical perspective circumventing the core–cusp problem envisages that no DM is present and tries to explain the galactic dynamics via a modification of gravity, as in the modified Newtonian dynamics (MOND) phenomenological approach (see Milgrom 1983, 2009; Bruneton & Esposito-Farèse 2007; Bekenstein 2004, 2009; for a review and further references, see Famaey & McGaugh 2012).

Here we take yet another viewpoint, retaining standard cold DM but envisaging that its dynamics may be subject to a non-minimal coupling with gravity (see Bruneton et al. 2009; Bertolami & Paramos 2010; Bettoni et al. 2011; Bettoni & Liberati 2015; Ivanov & Liberati 2020). The rationale of such attempts is trying to keep the collisionless DM phenomenology on large cosmological scales, and at the same time to introduce MOND-like behavior in galaxies by attributing to DM some non-negligible coupling with gravity. The word "non-minimal" simply means that such a coupling is embodied in the action via an additional interaction term with respect to standard scalar–tensor gravity theories. We stress that the non-minimal coupling is not necessarily a fundamental feature of the DM particles but might develop dynamically when the averaging/coherence length L associated with the fluid description of the collective behavior of matter is comparable to the local curvature scale. We specifically consider and theoretically justify a form of such a coupling that in the Newtonian limit amounts to a modification of the Poisson equation by a term L²∇² ρ proportional to the Laplacian of the DM density ρ itself, thus expressing an effective coupling of the DM fluid with the local gravitational curvature. Note that extensions of gravity including non-minimal couplings have also been investigated on cosmological scales, since under certain conditions they can mimic and explain the properties of the dark energy component (e.g., Bettoni et al. 2012; Bertolami & Paramos 2014).

In the present work, we show that non-minimally coupled DM halos in self-gravitating equilibria feature a density distribution that closely follows the Burkert profile out to several core scale radii. Moreover, we show that these can fit very well the measured rotation curves of dwarf, DM-dominated galaxies, and are consistent with the observed universal core surface density. In more detail, the plan of the paper is the following: in Section 2 we provide the theoretical motivation for the adopted non-minimally coupled DM framework; in Section 3 we introduce the basic formalism to describe self-gravitating equilibria of non-minimally coupled DM halos; in Section 3.1 we describe the effective equation of state for the DM fluid; in Section 3.2 we derive the fundamental equation ruling the DM density profile and study its solution space; in Section 4 we compare our non-minimally coupled solutions with the density profiles classically adopted in the literature to fit simulations and/or observations; in Section 5 we provide a first glimpse of how our non-minimally coupled DM mass distributions can fit the measured rotation curves of dwarf, DM-dominated galaxies; in Section 6 we show that non-minimally coupled DM halos are indeed consistent with the observed universal behavior of the core surface density; finally, in Section 7 we summarize our findings and highlight future prospects.

Throughout this work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration 2020) with rounded parameter values: matter density Ω_M = 0.3, dark energy density Ω_Λ = 0.7, baryon density Ω_b = 0.05, and Hubble constant H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.7; unless otherwise specified, G ≈ 6.67 × 10⁻⁸ cm³ g⁻¹ s⁻² indicates the standard gravitational (Newton) constant.

The motivation for introducing a non-minimal coupling between the DM matter field and gravity is twofold. On the theoretical side, it is allowed by the Einstein equivalence principle (Di Casola et al. 2015), and it might be required for the renormalizability of quantum field theories in curved spacetimes (e.g., Sonego & Faraoni 1993; Bruneton et al. 2009). On the observational side, it may help in explaining the striking relationships between DM and baryons recalled in Section 1 (e.g., baryonic Tully–Fisher, core radius versus disk scale-length, acceleration relations, etc.) that are not trivially understood via galaxy formation processes. In fact, such relationships may indicate a very special DM–baryon interaction. However, on the one hand, evidence of any direct coupling is missing. On the other hand, modified-gravity models such as MOND seem, at least in galactic environments, able to describe the baryon dynamics as that of a freely falling particle on some modified gravitational background.

Following these hints one can be led to conjecture, as in Bruneton et al. (2009), that the physical metric experienced by baryons may not coincide with the gravitational one, but it may be also determined by the properties of the DM field. Indeed, it was shown by Bekenstein (1993) that the most generic transformation between a physical and gravitational metric preserving causality and the weak equivalence principle can be of the general form

$$\begin{eqnarray}&&{\widetilde{g}}_{\mu \nu }={e}^{2\varphi }\left[{ \mathcal A }({ \mathcal X }){g}_{\mu \nu }+{ \mathcal B }({ \mathcal X }){{\rm{\nabla }}}_{\mu }\varphi {{\rm{\nabla }}}_{\nu }\varphi \right],\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal A }$ and ${ \mathcal B }$ are functions to be specified, φ is an extra scalar field and ${ \mathcal X }=-\tfrac{1}{2}\,{g}_{\mu \nu }{{\rm{\nabla }}}^{\mu }\varphi {{\rm{\nabla }}}^{\nu }\varphi ;$ such a relation between metrics is called a disformal one. Following the aforementioned idea, we can now ask which kind of interaction might be reduced to an effective coupling via a disformal metric of the above form with the scalar field playing the role of DM.

To this purpose, one can start from a general action of the form

$$\begin{eqnarray}\begin{array}{rcl}S & = & {S}_{\mathrm{EH}}[{g}_{\mu \nu }]+{S}_{\mathrm{bar}}[{g}_{\mu \nu },\psi ]\\ & & +\,{S}_{\mathrm{DM}}[{g}_{\mu \nu },\varphi ]+{S}_{\mathrm{int}}[{g}_{\mu \nu },\psi ,\varphi ],\end{array}\end{eqnarray} \tag{ 2 }$$

where the terms on the right-hand side identify, respectively, the Einstein–Hilbert (standard general relativity) action ${S}_{\mathrm{EH}}=\tfrac{{c}^{4}}{16\pi G}\,\int {d}^{4}x\,\sqrt{-g}\,R$ in terms of the Ricci scalar R, the baryonic and DM actions S_bar and S_DM, and the interaction one S_int. Here the scalar fields ψ and φ are thought of as collective variables encoding baryons and DM, respectively. It is then easy to show that DM can produce an effective metric for the baryons of the form specified by Equation (1) if S_bar[g_{μ ν} , ψ] + S_int[g_{μ ν} , ψ, φ] ≈ S_bar[g_{μ ν} + h_{μ ν} , ψ] with h_{μ ν} ∝ ∇_μ φ∇_ν φ. Up to order ${ \mathcal O }({h}^{2})$, this is created by an interaction term of the form ${S}_{\mathrm{int}}[{g}_{\mu \nu },\psi ,\varphi ]\propto \int {d}^{4}x\,\sqrt{-g}\,{T}_{\mathrm{bar}}^{\mu \nu }\,{{\rm{\nabla }}}_{\mu }\varphi {{\rm{\nabla }}}_{\nu }\varphi$, where ${T}_{\mathrm{bar}}^{\mu \nu }$ is the baryonic matter stress–energy tensor.

Note that, if now we express the full action in terms of the physical metric ${\widetilde{g}}_{\mu \nu }\equiv {g}_{\mu \nu }+{h}_{\mu \nu }$ (which is equivalent to choosing a frame in which baryons follow the geodesics of this metric, i.e., the Jordan frame), then it turns out that the DM field gets non-minimally coupled to gravity as

$$\begin{eqnarray}\begin{array}{rcl}S & = & {S}_{\mathrm{EH}}[{\widetilde{g}}_{\mu \nu }]+{S}_{\mathrm{bar}}[{\widetilde{g}}_{\mu \nu },\psi ]+{S}_{\mathrm{DM}}[{\widetilde{g}}_{\mu \nu },\varphi ]\\ & & +\,\epsilon {L}^{2}\,\displaystyle \int {d}^{4}x\,\sqrt{-\tilde{g}}\,{\widetilde{G}}^{\mu \nu }\,{{\rm{\nabla }}}_{\mu }\varphi {{\rm{\nabla }}}_{\nu }\varphi ,\end{array}\end{eqnarray} \tag{ 3 }$$

where L is a coupling length that must be present for dimensional consistency, = ±1 is a constant that represents the polarity of the coupling (undetermined a priori), ${\widetilde{G}}^{\mu \nu }$ is the Einstein tensor expressed in terms of the physical metric ${\widetilde{g}}_{\mu \nu }$, and the DM field has been implicitly redefined through a conformal factor.

Three remarks are in order. First, note that the only other non-minimal coupling term with the same physical dimensions (still leading to second-order field equations) would be proportional to the Ricci scalar as ${ \mathcal X }\,\widetilde{R}$, which, however, is equivalent to that appearing in Equation (3) modulo a surface term (see Bettoni & Liberati 2013). Second, the coupling term ${\widetilde{G}}^{\mu \nu }\,{{\rm{\nabla }}}_{\mu }\varphi {{\rm{\nabla }}}_{\nu }\varphi$ appearing in the above action is proportional to a term of the Horndeski Lagrangian, which constitutes the most general scalar–tensor theory giving rise to second-order field equations; however, here we are not proposing a fundamental theory of modified gravity, but just entailing the possibility that DM in galactic halos dynamically develops a non-minimal coupling with the metric characterized by an effective length scale L. Third, it can be shown that by choosing a special form of the tensor h_{μ ν} , one can reproduce the MONDian regime in galaxies (see Bruneton et al. 2009, their Equations (2.6)–(2.7)). In the present paper, we rely on the much simpler form h_{μ ν} ∝ ∇_μ φ∇_ν φ as above, which will not provide a MONDian limit, but will produce kinematics consistent with observations (at least in dwarf, DM-dominated galaxies); from this point of view, the coupling length L cannot be straightforwardly interpreted in terms of the acceleration parameter a₀ appearing in the MONDian dynamics.

In order to take the Newtonian limit of the above theory it is convenient to adopt fluid variables. In the case of a complex scalar field ⁵ this is easily achieved by adopting the standard Madelung representation $\varphi \sim \sqrt{\rho }{e}^{i\theta }$. However, a fluid limit is possible also for a non-minimally coupled real scalar field (see Bettoni et al. 2012). In the end it can be shown that the non-minimal coupling considered above leads to a modified Poisson equation of the form (see, e.g., Bettoni et al. 2014)

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}^{2}{\rm{\Phi }}=4\pi G\,[(\rho +{\rho }_{\mathrm{bar}})-\epsilon \,{L}^{2}\,{{\rm{\nabla }}}^{2}\rho ],\end{eqnarray} \tag{ 4 }$$

where Φ is the Newtonian potential, and ρ_bar and ρ are the baryon and DM mass densities. This modified Poisson equation implies that the source for gravity is not just the total matter density itself, but also the DM inhomogeneities or spatial variations. ⁶ In the present paper we will focus mainly on dwarf, DM-dominated galaxies and thus we will neglect the baryonic component ρ_bar hereafter.

While the discussed non-minimal coupling could be associated with some modified theory of gravity, in the extant literature at least two main mechanisms have been contemplated to produce dynamically a non-minimal coupling characterized by a length scale L: either it can emerge from some collective behavior of the DM particles associated with a coherence length (for example, via a Bose–Einstein condensation mechanism; see Bettoni et al. 2011 for an extended discussion) or it appears through an averaging procedure associated with the fluid description of matter. Although in general we adhere to the latter viewpoint, discussing the physical emergence of L is beyond the scope of this paper; indeed in what follows we shall investigate the implications of the modified Poisson Equation (4) for the self-gravitating equilibria of DM halos, without an a priori prejudice about the origin and scale of L.

The self-gravitating equilibria of DM halos can be specified in the fluid approximation (Teyssier et al. 1997; Subramanian 2000; Lapi & Cavaliere 2011; Nadler et al. 2017) via the continuity, Euler (also called Jeans in this context), and Poisson equations

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\partial }_{t}\,\rho +{\boldsymbol{\nabla }}\cdot (\rho \,{\boldsymbol{v}})=0,\\ {\partial }_{t}\,{\boldsymbol{v}}+({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{v}}+\tfrac{1}{\rho }\,{\boldsymbol{\nabla }}\,p=-{\boldsymbol{\nabla }}\,{\rm{\Phi }},\\ {{\boldsymbol{\nabla }}}^{2}{\rm{\Phi }}=4\pi G(\rho -\epsilon \,{L}^{2}\,{{\rm{\nabla }}}^{2}\rho ).\end{array}\right.\end{eqnarray} \tag{ 5 }$$

Here v is the bulk velocity, $p=\rho \,{\sigma }_{r}^{2}$ is the pressure dynamically generated by the random motions (and specified in terms of a radial velocity dispersion ${\sigma }_{r}^{2}$ or more generally an anisotropic stress tensor ${\sigma }_{{ij}}^{2}$) of the DM particles in approximate virial equilibrium within the gravitational potential well ⁷ . The last term on the right-hand side of the Poisson equation represents the non-minimal coupling with length scale L and polarity as discussed in Section 2; we will see that for the purpose of originating physically acceptable DM density distributions, a negative polarity = −1 is required, while L will turn out to be closely related to the DM core radius. Note that non-minimal coupling terms do not appear in the Euler equation since they are found to be sub-leading in the non-relativistic limit (i.e., expansion in 1/c²; see Bettoni et al. 2014).

3.1. Equation of State

To close the system, the pressure must be related to the density via an equation of state (EOS). In the present context, we focus on the EOS created by the gravitational assembly of DM halos via accretion and mergers from the cosmic web. This EOS stems from the progressive stratification of the DM pseudo-entropy K ≡ p/ρ^5/3, or equivalently of the coarse-grained phase-space density $\rho /{\sigma }_{r}^{3}\propto {K}^{-3/2}$, in terms of a simple power-law profile K(r) ∝ r^α . Although the physical origin of this scale-free behavior is not fully understood (see Nadler et al. 2017; Arora & Williams 2020), N-body simulations (Peirani et al. 2006; Navarro et al. 2010; Ludlow et al. 2011; Gao et al. 2012; Butsky et al. 2016; Nolting et al. 2016) have shown this to approximately hold over more than three order of magnitude in radius within virialized halos, with power-law index α ≈ 1.25–1.3. Such values are indeed expected on the basis of simple self-similar arguments (see Bertschinger 1985; Lapi & Cavaliere 2009a, 2011; Nadler et al. 2017) and are also broadly consistent with observations (see Lapi & Cavaliere 2009b; Chae 2014; Munari et al. 2014). On this basis, but to keep some degree of generality, ⁸ we adopt the EOS parameterization

$$\begin{eqnarray}&&p(\rho ,r)=\lambda \,{\rho }^{{\rm{\Gamma }}}\,{r}^{\alpha }={\rho }_{\alpha }\,{\sigma }_{\alpha }^{2}{\left(\displaystyle \frac{\rho }{{\rho }_{\alpha }}\right)}^{{\rm{\Gamma }}}{\left(\displaystyle \frac{r}{{r}_{\alpha }}\right)}^{\alpha },\end{eqnarray} \tag{ 6 }$$

where r_α is a reference radius and we have defined ρ_α ≡ ρ(r_α ) and ${\sigma }_{\alpha }^{2}\equiv {\sigma }_{r}^{2}({r}_{\alpha });$ in the following we adopt as fiducial values ${\rm{\Gamma }}\approx \tfrac{5}{3}$ and α ≈ 1.3.

Under static ( v = 0), spherically symmetric, and isotropic conditions (see the Appendix for a generalization) the relevant Equations (5) become

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{1}{\rho }\,\tfrac{{dp}}{{dr}}=-\tfrac{d{\rm{\Phi }}}{{dr}},\\ \tfrac{1}{{r}^{2}}\,\tfrac{d}{{dr}}\left({r}^{2}\,\tfrac{d{\rm{\Phi }}}{{dr}}\right)=4\pi G\left[\rho -\epsilon \,{L}^{2}\,\tfrac{1}{{r}^{2}}\,\tfrac{d}{{dr}}\left({r}^{2}\,\tfrac{d\rho }{{dr}}\right)\right],\end{array}\right.\end{eqnarray} \tag{ 7 }$$

supplemented with the trivial mass conservation constraint ${ \mathcal M }(t)=$ const, where

$$\begin{eqnarray}&&{ \mathcal M }\equiv 4\pi {\int }_{0}^{\infty }{dr}\,{r}^{2}\,\rho (r),\end{eqnarray} \tag{ 8 }$$

is the total mass ${ \mathcal M }$ of the DM halo.

3.2. The Fundamental Equation and its Solutions

It is convenient to introduce normalized variables $\bar{r}=r/{r}_{\alpha }$ $\bar{\rho }\equiv \rho /{\rho }_{\alpha }$ and define the quantities $\kappa \equiv 4\pi G\,{\rho }_{\alpha }\,{r}_{\alpha }^{2}/{\sigma }_{\alpha }^{2}$ and $\eta \equiv \epsilon \,{L}^{2}/{r}_{\alpha }^{2}$. To understand the physical meaning of κ, one can choose ${r}_{\alpha }={r}_{\max }$ to be the point at which the circular velocity ${v}_{c}^{2}(r)\equiv G\,M(\lt r)/r$ peaks at a value ${v}_{c}^{2}({r}_{\max })\,=4\pi \rho ({r}_{\max }){r}_{\max }^{2}$, so that $\kappa ={v}_{c}^{2}({r}_{\max })/{\sigma }^{2}({r}_{\max })$ is seen to compare the estimate ${\sigma }^{2}({r}_{\max })$ for the random kinetic energy with that, ${v}_{c}^{2}({r}_{\max })$, for the gravitational potential.

Eliminating dΦ/dr from Equations (7) and using the EOS Equation (6) yields the following fundamental equation for the density:

$$\begin{eqnarray}&&\bar{\rho }^{\prime\prime} +\displaystyle \frac{({\rm{\Gamma }}-2)\tfrac{\bar{\rho }{{\prime} }^{2}}{\bar{\rho }}+\tfrac{\alpha (2{\rm{\Gamma }}-1)+2\,{\rm{\Gamma }}}{{\rm{\Gamma }}}\,\tfrac{\bar{\rho }^{\prime} }{\bar{r}}+\tfrac{\alpha (\alpha +1)}{{\rm{\Gamma }}}\,\tfrac{\bar{\rho }}{{\bar{r}}^{2}}-2\,\eta \,\kappa \,\tfrac{{\bar{\rho }}^{2-{\rm{\Gamma }}}\,\bar{\rho }^{\prime} }{{\rm{\Gamma }}\,{\bar{r}}^{\alpha +1}}+\kappa \,\tfrac{{\bar{\rho }}^{3-{\rm{\Gamma }}}}{{\rm{\Gamma }}\,{\bar{r}}^{\alpha }}}{1-\eta \,\kappa \,\tfrac{{\bar{\rho }}^{2-{\rm{\Gamma }}}}{{\rm{\Gamma }}\,{\bar{r}}^{\alpha }}}=0,\end{eqnarray} \tag{ 9 }$$

while the mass conservation constraint now reads ${ \mathcal M }={\rho }_{\alpha }\,{r}_{\alpha }^{3}\,{f}_{{ \mathcal M }}$ in terms of the shape factor ${f}_{{ \mathcal M }}=4\pi \,{\int }_{0}^{\infty }d\bar{r}\,{\bar{r}}^{2}\,\bar{\rho }(\bar{r})$. The solution space of such an equation is amazingly rich, and for ${\rm{\Gamma }}=\tfrac{5}{3}$ and η = 0 it has been quite extensively studied in the literature to describe the radial structure of standard ΛCDM halos (see Williams et al. 2004; Hansen 2004; Austin et al. 2005; Dehnen & McLaughlin 2005; Lapi & Cavaliere 2009a). In the following we provide the generalization with a generic Γ and then with the addition of the non-minimal coupling.

To understand the general features of the solutions, it is convenient to look for power-law behaviors $\bar{\rho }\simeq {\bar{r}}^{-\gamma }$. Substituting in the fundamental equation yields

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Gamma }}({\rm{\Gamma }}-1)\left(\gamma -\displaystyle \frac{\alpha }{{\rm{\Gamma }}}\right)\left(\gamma -\displaystyle \frac{\alpha +1}{{\rm{\Gamma }}-1}\right)\\ \quad -\,\eta \,\kappa \,\displaystyle \frac{\gamma (\gamma -1)}{{\bar{r}}^{\gamma (2-{\rm{\Gamma }})+\alpha }}=-\displaystyle \frac{\kappa }{{\bar{r}}^{\gamma (2-{\rm{\Gamma }})+\alpha -2}}.\end{array}\end{eqnarray} \tag{ 10 }$$

Focusing first on the minimally coupled case with η = 0, it is evident that trivial power-law solutions with slope $\gamma ={\gamma }_{\alpha }\,\equiv \tfrac{2-\alpha }{2-{\rm{\Gamma }}}$ are admitted, implying $\kappa ={\kappa }_{{\rm{PL}}}\equiv 2\,\tfrac{\left({\rm{\Gamma }}-\alpha )(\alpha +4-3{\rm{\Gamma }}\right)}{{\left(2-{\rm{\Gamma }}\right)}^{2}}$. These values are actually not physically acceptable at small and large radii since the gravitational force and mass would diverge, but they provide the behavior of any solution at intermediate radii. In fact, to find general solutions it is numerically convenient to choose r_α as the reference radius where the logarithmic slope of the density is –γ_α and then integrate inward and outward. In terms of the normalized variables, this corresponds to setting the boundary conditions as $\bar{\rho }(1)=1$ and $\bar{\rho }^{\prime} (1)=-{\gamma }_{\alpha }$. When integrating inward of r_α , it is found that the density profile features a physically acceptable behavior only for a specific value of κ = κ_Γ,α , somewhat different from the κ_PL defined above. In particular, for such a value κ_Γ,α the profile asymptotes for $\bar{r}\ll 1$ to $\bar{\rho }\propto {\bar{r}}^{-{\gamma }_{0}}$ with ${\gamma }_{0}\equiv \alpha /{\rm{\Gamma }}$, consistently with the power-counting in Equation (10). For κ > κ_Γ,α the profile has wiggles (change of sign in the second derivative) and steepens toward the center to imply a diverging gravitational force, while for κ < κ_Γ,α it develops a central hole. The optimal value κ_Γ,α can be found by a shooting technique (i.e., automatically solving inward the differential equation with different slopes –γ_α at r_α until the desired, physical inner asymptotic behavior is found), as represented in Figure 1.

Zoom In Zoom Out Reset image size

Figure 2 illustrates the resulting full density profiles for three different values of α. The outer behavior of the profile turns out to be physical only for $\alpha \leqslant {\alpha }_{\mathrm{crit}}\equiv \tfrac{{\rm{\Gamma }}(5{\rm{\Gamma }}-6)}{3{\rm{\Gamma }}-2}$. For α > α_crit the outer slope is too flat to imply a diverging mass. For α < α_crit the solution $\bar{\rho }\propto {\bar{r}}^{-{\gamma }_{\infty }}$ attains a slope ${\gamma }_{\infty }\equiv \tfrac{\alpha +1}{{\rm{\Gamma }}-1}$ at a finite large radius before an outer cutoff. For α = α_crit the cutoff is pushed to infinity and the slope ${\gamma }_{\infty ,\mathrm{crit}}\equiv \tfrac{5{\rm{\Gamma }}+2}{3{\rm{\Gamma }}-2}$ is attained only asymptotically for $\bar{r}\to \infty$, consistently with the power-counting from Equation (10); correspondingly, the intermediate and inner slopes read ${\gamma }_{\alpha ,\mathrm{crit}}\equiv \tfrac{5{\rm{\Gamma }}-2}{3{\rm{\Gamma }}-2}$ and ${\gamma }_{0,\mathrm{crit}}\equiv \tfrac{5{\rm{\Gamma }}-6}{3{\rm{\Gamma }}-2}$, respectively. We illustrate these different behaviors in Figure 2 for the fiducial value ${\rm{\Gamma }}=\tfrac{5}{3}$. In such a case, the power-law solutions have γ_α = 6 − 3α and κ_PL = 6(5 − 3α)(α − 1), while the physical solutions feature an inner slope ${\gamma }_{0}=\tfrac{3\,\alpha }{5}$ and an outer slope ${\gamma }_{\infty }=\tfrac{3(1+\alpha )}{2};$ the values of κ_α,Γ are reported as a function of α in the inset of Figure 2. For $\alpha ={\alpha }_{\mathrm{crit}}=\tfrac{35}{27}$ one gets ${\gamma }_{0,\mathrm{crit}}\equiv \tfrac{7}{9}$, ${\gamma }_{\alpha ,\mathrm{crit}}=\tfrac{19}{9}$, ${\gamma }_{\infty ,\mathrm{crit}}=\tfrac{31}{9}$, and ${\kappa }_{\alpha ,{\rm{\Gamma }}}=\tfrac{200}{81}\approx 2.5$.

Zoom In Zoom Out Reset image size

When including the non-minimal coupling $\eta \equiv \epsilon \,\tfrac{{L}^{2}}{{r}_{\alpha }^{2}}$, the solution space changes appreciably. First of all, (hence η) must be negative, otherwise the inner profile diverges at a finite radius. Then for any negative value of η, there is again an optimal value of κ = κ_Γ,α,η such that the inner profile is physical, with limiting central slope γ₀ = 0, i.e., a core. This is again in accordance with the power-counting in Equation (10) since now the second term on the left-hand side dominates the behavior for small $\bar{r}\ll 1$. Other solutions with κ smaller or larger than κ_Γ,α,η are not physically acceptable since they have non-monotonic behaviors with the density first flattening and then steepening toward a central slope γ₀ = 1, or they develop a central hole. As for the outer behavior, the profile has a cutoff at a finite radius ${ \mathcal R }$, setting the effective halo boundary, which is smaller for more negative values of η. We illustrate the physically acceptable profiles for different η and the related values of the constant κ_Γ,α,η in Figure 3. The corresponding distributions of mass, circular velocity, and velocity dispersion for a few values of the non-minimal coupling η are also illustrated in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now compare the shape of our physical solutions to some classic literature density profiles, characterized by different analytic expressions and numbers of parameters, that are commonly adopted to fit simulations and/or observations (see also Figure 15 of Freundlich et al. 2020b for a comprehensive account). To this purpose, it is convenient to use a radial coordinate $\hat{r}\equiv r/{r}_{-2}$ normalized to the radius r₋₂ where the logarithmic density slope $\tfrac{d\mathrm{log}\rho }{d\mathrm{log}r}=-2$, and to rescale the density profile $\hat{\rho }\equiv \rho /\rho ({r}_{-2})$ accordingly. Specifically, we will consider the following density profiles.

• 1.  
  α β γ profiles.The α β γ profiles (Zhao 1996; see also Widrow 2000) feature the shape

  $$\begin{eqnarray}&&\hat{\rho }(\hat{r})={\hat{r}}^{-\tau }{\left(\displaystyle \frac{1+w}{1+w{\hat{r}}^{\omega }}\right)}^{\xi },\end{eqnarray} \tag{ 11 }$$

  where the three parameters τ, ω, and ξ describe, respectively, the central slope, the middle curvature, and the outer decline of the density run, while $w\equiv -\tfrac{2-\tau }{(2-\tau -\omega \,\xi )}$. Familiar empirical profiles are recovered for specific values of the triplet (τ, ω, ξ): e.g., Plummer's profile corresponds to (0, 2, 5/2), Jaffe's to (2, 1, 2), and Hernquist's to (1, 1, 3). The standard NFW profile (Navarro et al. 1996), which is classically used to fit N-body simulations in the ΛCDM model, is obtained for the parameter triple (1, 1, 2); for a generalization with different inner slope, often referred to as gNFW, the parameters (τ, 1, 3 − τ) apply (Mamon et al. 2019; for a more complex cored version see also Read et al. 2016). It is worth mentioning that recently Freundlich et al. (2020b) have considered an α β γ profile (referred to also as the Zhao–Dekel model) with parameters $(\tau ,\tfrac{1}{2},7-2\tau )$, which can provide good fits to the density profiles from N-body, DM-only, and hydro simulations including baryonic effects.
• 2.  
  Sérsic–Einasto profile.The Sérsic–Einasto profile (see An & Zhao 2013) is defined as

  $$\begin{eqnarray}&&\hat{\rho }(\hat{r})={\hat{r}}^{-\tau }\,{e}^{-u({\hat{r}}^{\omega }-1)},\end{eqnarray} \tag{ 12 }$$

  where τ is the inner density slope, ω is a shape parameter, and $u\equiv \tfrac{2-\tau }{\omega }$. The classic cored Einasto shape (Sersic 1963; Einasto 1965; Prugniel & Simien 1997; Graham et al. 2006; also Lazar et al. 2020 for a more complex analytical expression) is recovered for τ = 0. N-body simulations in the standard ΛCDM model are usually well described by the parameter values ω ≈ 0.15–0.2 and τ ≈ 0 (although values τ ≲ 0.8 are not ruled out given the resolution of current simulations).
• 3.  
  Soliton profile.The soliton profile features the shape (see Schive et al. 2014a, 2014b)

  $$\begin{eqnarray}&&\hat{\rho }(\hat{r})={\left(\displaystyle \frac{1+\omega /{\hat{r}}_{c}^{2}}{1+\omega {\hat{r}}^{2}/{\hat{r}}_{c}^{2}}\right)}^{8},\end{eqnarray} \tag{ 13 }$$

  where the normalized core radius reads ${\hat{r}}_{c}=\sqrt{7\,\omega }$. In numerical simulations of non-interacting Bose–Einstein condensate DM (alias fuzzy DM), halos are described by a combination of this solitonic profile with ω ≈ 0.091 in the inner region and an NFW profile in the outskirts.
• 4.  
  BEC profile.The profile followed by an interacting Bose–Einstein condensate (BEC; see Bohmer & Harko 2007; Harko 2014) in the Thomas–Fermi limit reads

  $$\begin{eqnarray}&&\hat{\rho }(\hat{r})=\displaystyle \frac{1}{\hat{r}}\,\displaystyle \frac{\sin (\pi \,\hat{r}/\hat{R})}{\sin (\pi /\hat{R})},\end{eqnarray} \tag{ 14 }$$

  where the normalized halo boundary is defined by the equality $\pi /\hat{R}+\tan (\pi /\hat{R})=0$. Incidentally, note that such a profile actually corresponds to the solution of our Equation (10) for Γ = 2 and α = η = 0.
• 5.  
  Burkert profile.The Burkert profile can be written as (see Burkert 1995; Salucci & Burkert 2000)

  $$\begin{eqnarray}&&\hat{\rho }(\hat{r})=\displaystyle \frac{(1+{\hat{r}}_{c})(1+{\hat{r}}_{c}^{2})}{(\hat{r}+{\hat{r}}_{c})({\hat{r}}^{2}+{\hat{r}}_{c}^{2})},\end{eqnarray} \tag{ 15 }$$

  where ${\hat{r}}_{c}$ is the normalized core radius, defined by the nonlinear algebraic equation $2\,{\hat{r}}_{c}^{3}+{\hat{r}}_{c}^{2}-1=0$. This is a cored profile generally exploited to fit observations of (dwarf) spiral galaxies.

In Figure 5 we compare the shape of our solutions (for definiteness ${\rm{\Gamma }}=\tfrac{5}{3}$ and $\alpha ={\alpha }_{\mathrm{crit}}=\tfrac{35}{27}$ are adopted) to a few of the above profiles. First it is evident that our density run for η = 0 describes quite well the NFW and Einasto profiles commonly used to fit N-body simulations in the standard ΛCDM cosmology. For −0.1 < η < −0.01 our solutions develop a core and the shape out to a few/several r₋₂ is remarkably close to the Burkert profile, commonly exploited to fit observations of dwarf galaxies. At larger distances there is progressive deviation from the Burkert profile, since the latter has been designed to have a limiting slope close to the NFW one, while our profiles get truncated at a finite radius ${ \mathcal R }$. However, for an appreciable range of η values this is not a concern since the truncation occurs at radii much beyond r₋₂ that are scantily (if at all) probed by observations (see also Section 5). We stress that it is a remarkable property of the non-minimally coupled solutions to reproduce the Burkert shape not only in the inner region, but also over an extended radial range outward of r₋₂. For comparison, other models such as those based on Bose–Einstein condensate DM predict a cored density profile, but the deviation from the Burkert shape is appreciable, in terms of a prominent cutoff or steep decline, just outside the core. Note that our solutions for η = −0.05 can be reasonably described by an α β γ model (see Equation (11)) with parameters $(\tau ,\omega ,\xi )=\left(0,\tfrac{3}{2},\tfrac{5}{2}\right)$ or by a Sérsic–Einasto model (see Equation (12)) with parameters $(\tau ,\omega )=\left(0,\tfrac{3}{4}\right)$ up to several r₋₂ before the final cutoff.

Zoom In Zoom Out Reset image size

We now compare the non-minimally coupled mass distributions with observed galaxy rotation curves. To avoid dealing with baryons (which are not included in our treatment), we focus on dwarf, strongly DM-dominated galaxies. Specifically, we exploit the coadded rotation curve built by Lapi et al. (2018) based on the high-quality measurements of about 20 dwarf galaxies with I-band magnitude M_I ≳ −18.5 from the original sample by Persic & Salucci (1996). The sample features an average disk scale-length R_e ≲ 1.5 kpc, disk mass M_⋆ ≲ 10⁹ M_⊙, and halo mass ${ \mathcal M }\lesssim {10}^{11}\,{M}_{\odot }$. The coadded rotation curve is well measured out to a galactocentric distance of about 7 kpc. We then use a Levenberg–Marquardt least-squares minimization routine to fit the measured rotation curve with the NFW and Burkert profiles, and with our physical solutions for different values of the non-minimal coupling parameter η; the outcomes are illustrated in Figure 6.

Zoom In Zoom Out Reset image size

As is well known, the NFW fit (equivalent to our solution with η = 0) struggles to fit the measured dwarf galaxy rotation curves, yielding a reduced χ² ≈ 4. On the other hand, the Burkert profile performs much better, providing a good fit with a reduced χ² ≈ 0.58. Our non-minimally coupled solution for η ≈ −0.05 provides a fit of quality comparable to the Burkert one, yielding practically the same reduced χ² ≈ 0.61. The current data are compatible within 3σ with any value of η in the range from −0.1 to −0.01. Accurate determinations of the coadded rotation curve out to 10 kpc or beyond would be necessary to determine the non-minimal coupling parameter η to a good level of precision. Interestingly, one can also see from Figure 6 that there is a tendency to favor values of η slightly less negative than −0.05, which is actually what is to be expected in dwarf galaxies with halo masses ${ \mathcal M }\lesssim {10}^{11}\,{M}_{\odot }$ on the basis of universal scaling arguments (see Section 6 and in particular Equation (19)).

From the fit with η ≈ −0.05 (or equivalently from the Burkert one) we derive a best-fit value of the core radius around r₀ ≈ 2.9 ± 0.1 kpc, which is about twice the average disk scale-length of the systems in the considered sample. This is remarkably consistent with the empirical, yet still puzzling, relationship between core radius and disk scale-length as determined by Donato et al. (2004).

A remark is in order concerning the good agreement of our simple non-minimally coupled DM model with the observed rotation curve. On the one hand, Bruneton et al. (2009) have shown that choosing an appropriate yet rather complex shape (see their Equations (2.6)–(2.7)) of the coupling h_{μ ν} characterizing the interaction term in the action of Equation (3) can lead to a MONDian phenomenology on galactic scales, which is well known to reproduce rotation curves. However, our model is based on a much simpler coupling h_{μ ν} ∝ ∇_μ φ∇_ν φ, and hence does not lead exactly to MONDian dynamics in its Newtonian limit; from this perspective, the rather good performance of our simple non-minimally coupled model in reproducing the measured rotation curves is valuable and far from trivial.

As mentioned in Section 1, it has been well established observationally (see Salucci & Burkert 2000; Burkert 2015) that, at least for dwarf galaxies with halo masses ${ \mathcal M }\lesssim {10}^{11}\,{M}_{\odot }$, the product of the core density and core radius (i.e., a sort of core surface density) is an approximately universal constant with values ${\rho }_{0}\,{r}_{0}\approx {75}_{-45}^{+55}\,{M}_{\odot }$ pc⁻² among different galaxies, see Figure 7. This somewhat unexpected property poses a serious challenge to any theoretical model of core formation (e.g., Deng et al. 2018; Burkert 2020). We now aim to show that non-minimally coupled DM halos are instead consistent with such a remarkable scaling law.

Zoom In Zoom Out Reset image size

First, adopting for definiteness ${\rm{\Gamma }}=\tfrac{5}{3}$ and $\alpha ={\alpha }_{\mathrm{crit}}=\tfrac{35}{27}$, we compute the physical solutions of Equation (9) subject to the boundary conditions ρ(r_α ) = ρ_α and $\rho ^{\prime} ({r}_{\alpha })=-{\gamma }_{\alpha }\,\tfrac{{\rho }_{\alpha }}{{r}_{\alpha }}$, for several values of the non-minimal coupling parameter η. From the normalized profiles thus obtained we fit as a function of $| \eta | =\tfrac{{L}^{2}}{{r}_{\alpha }^{2}}$ in the range −0.1 ≲ η ≲ −0.01 the following relations involving the normalized boundary radius ${ \mathcal R }$, the core radius r₀, the core density ρ₀, and the mass shape factor ${f}_{{ \mathcal M }}$ defined below Equation (9):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{ \mathcal R }}{{r}_{\alpha }} & \simeq & 0.8\,| \eta {| }^{-0.7},\,\,\,\,\,\,{f}_{{ \mathcal M }}\simeq 9.1\,| \eta {| }^{-0.4},\\ \displaystyle \frac{{\rho }_{0}}{{\rho }_{\alpha }} & \simeq & 0.9\,| \eta {| }^{-0.8},\,\,\,\,\,\,\displaystyle \frac{{r}_{0}}{{r}_{\alpha }}\simeq 1.6\,| \eta {| }^{0.5}.\end{array}\end{eqnarray} \tag{ 16 }$$

We then combine the above scaling with the expression for the total mass ${ \mathcal M }={f}_{{ \mathcal M }}\,{\rho }_{\alpha }\,{r}_{\alpha }^{3}$ and with the definition of the virial radius ${ \mathcal R }={(3{ \mathcal M }/4\pi {{\rm{\Delta }}}_{\mathrm{vir}}{\rho }_{{\rm{c}}}{E}_{z})}^{1/3}\approx 120\,{E}_{z}^{-1/3}({ \mathcal M }/{10}^{11}$ M_⊙)^1/3; here ρ_c ≈ 2.8 × 10¹¹ h² M_⊙ Mpc⁻³ is the critical density, E_z ≡ Ω_M (1 + z)³ + Ω_Λ takes into account the formation redshift z of the halo, and Δ_vir is the nonlinear threshold for virialization, with values around 100 at z ≈ 0 and increasing toward 180 for z ≳ 1. We eventually derive

$$\begin{eqnarray}\begin{array}{rcl}{r}_{0} & \simeq & 1.6\,L,\,\,\,\,\,\,\,{r}_{\alpha }\simeq 1.1\,{L}^{0.6}\,{{ \mathcal R }}^{0.4}\\ {\rho }_{0} & \simeq & 0.3\,{{\rm{\Delta }}}_{\mathrm{vir}}\,{\rho }_{{\rm{c}}}\,{E}_{z}{\left(\displaystyle \frac{{ \mathcal R }}{L}\right)}^{2.1}={\rho }_{\alpha }{\left(\displaystyle \frac{{ \mathcal R }}{L}\right)}^{0.6}\,;\end{array}\end{eqnarray} \tag{ 17 }$$

remarkably, the core radius r₀ turns out to be proportional, with a coefficient of order 1, to the non-minimal coupling length scale L.

To proceed further we need a relation between the core radius r₀ (or L) and the halo mass ${ \mathcal M };$ this is thought not to be of fundamental nature but rather to stem from two other relationships involving the baryonic mass: (i) the relation between the stellar (disk) mass and the halo mass (e.g., Moster et al. 2013), which is known to originate from baryonic processes related to galaxy formation; (ii) the relation between the core radius r₀ and the disk scale-length (in turn related to the stellar mass; e.g., Donato et al. 2004), which rather is still not completely understood. In the present pilot study we are not including the baryonic component and we cannot infer the ${r}_{0}-{ \mathcal M }$ relation from first principles; thus we will adopt the outcome ${r}_{0}\approx 4.5\,{({ \mathcal M }/{10}^{11}{M}_{\odot })}^{0.6}$ kpc from the dynamical modeling study by Salucci et al. (2007), and see what this implies for the core surface density. Specifically, from Equations (17) we obtain

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}\equiv {\rho }_{0}\times {r}_{0}\approx 50\left(\displaystyle \frac{{{\rm{\Delta }}}_{\mathrm{vir}}}{100}\right){E}_{z}^{0.3}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}\,,\end{eqnarray} \tag{ 18 }$$

which is independent of the halo mass and/or core radius, and only weakly dependent on formation redshift. In Figure 7 we report the above for three values of the formation redshift, z ≈ 0, 0.5, 1.5, finding it remarkably consistent with the average observed relation and its scatter. We stress the extreme relevance of this finding: the universality of the core surface density has proven to be extremely challenging for alternative DM models, even for those that are barely consistent with the ${r}_{0}-{ \mathcal M }$ relation we have assumed. For example, as pointed out by Burkert (2020), fuzzy DM can reproduce the ${r}_{0}-{ \mathcal M }$ relation, albeit with some (uncertain) hypothesis on core formation redshift. However, such a model is considerably out of track as to the core surface density scaling, since it robustly predicts ${\rho }_{0}\propto {r}_{0}^{-4}$. The same issue concerns many other DM models inspired by particle physics, as extensively discussed, e.g., by Deng et al. (2018).

As an aside, from Equations (7) and the adopted ${r}_{0}-{ \mathcal M }$ relation we can also derive three other interesting scaling laws. First, the dependence of η on halo mass reads

$$\begin{eqnarray}&&| \eta | \approx 0.04{\left(\displaystyle \frac{{ \mathcal M }}{{10}^{11}{M}_{\odot }}\right)}^{0.35}\,{E}_{z}^{0.3}\,;\end{eqnarray} \tag{ 19 }$$

this confirms that values of the non-minimal coupling in the range η = −0.1 to −0.01 cover the typical mass range of dwarf galaxies. We stress that this dependence of η on halo mass/formation redshift will induce slightly different shapes in the profiles, implying a weak violation of self-similarity. Though challenging, it will be interesting to look for such behaviors in real data (see Section 5). Second, it may be interesting to derive the dependence on the coupling, and hence on mass, of the inner logarithmic slope γ_0.1 measured at a reference radius of r ≈ 0.1r₋₂ ∼ a few percent of ${ \mathcal R }$. We get the scaling ∣γ_0.1∣ ≃ 0.035∣η∣^−0.55, which after Equation (19) translates into a mass dependence

$$\begin{eqnarray}&&| {\gamma }_{0.1}| \simeq 0.2{\left(\displaystyle \frac{{ \mathcal M }}{{10}^{11}{M}_{\odot }}\right)}^{-0.2}\,{E}_{z}^{-0.16}.\end{eqnarray} \tag{ 20 }$$

Thus there is a slight tendency for less massive halos to have flatter profile at a fixed radius in the inner region (note that asymptotically at the center all the non-minimally coupled DM density profiles are flat). In a future work, it would be interesting to investigate how such a scaling is altered by the presence of baryons in halos of different masses, and how the outcome will compare with the results from ΛCDM hydrodynamical simulations including feedback effects (e.g., Tollet et al. 2016; Freundlich et al. 2020b), which show a nontrivial mass dependence for the inner shape of the halo. Third, we can compute the halo concentration as ${c}_{\alpha }\equiv { \mathcal R }/{r}_{\alpha }$ (using ${ \mathcal R }/{r}_{-2}$ yields a similar result), which turns out to be

$$\begin{eqnarray}&&{c}_{\alpha }\simeq 10{\left(\displaystyle \frac{{ \mathcal M }}{{10}^{11}{M}_{\odot }}\right)}^{-0.15}\,{E}_{z}^{-0.2},\end{eqnarray} \tag{ 21 }$$

in broad agreement with, and actually slightly smaller than, the outcome of N-body, DM-only simulations in the standard ΛCDM cosmology (e.g., Bullock et al. 2001; Maccio et al. 2007).

We have investigated self-gravitating equilibria of halos constituted by DM non-minimally coupled to gravity. A non-minimal coupling may be present in modified-gravity theories or it might be dynamically generated when the averaging/coherence length L associated with a fluid description of the DM collective behavior is comparable to the local curvature scale. We have provided a theoretical motivation for a form of such a coupling that in the Newtonian limit amounts to a modification of the Poisson equation by a term L²∇² ρ proportional to the Laplacian of the DM density ρ itself (see Section 2). We have further adopted an effective power-law equation of state p ∝ ρ^Γ r^α relating the DM dynamical pressure p to density ρ and radius r, as expected from phase-space density stratification during the gravitational assembly of halos in a cosmological context (see Section 3.1). In the absence of the non-minimal coupling, we have confirmed previous findings that the DM density run ρ(r) features a central density cusp and an overall shape mirroring the outcomes of N-body simulations in the standard ΛCDM cosmology, as described by the classic NFW or Einasto profiles (see Section 3.2).

We have remarkably found that, when the non-minimal coupling is switched on, it causes the DM density profile to develop an inner core and a shape closely following the Burkert profile out to several core scale radii (see Section 4). In addition, we have highlighted that our non-minimally coupled solutions can fit, with an accuracy comparable to the Burkert profile, the coadded rotation curves of DM-dominated dwarf galaxies (see Section 5). Finally, we have shown that non-minimally coupled DM halos are consistent with the observed scaling relation between the core radius r₀ and the core density ρ₀ in terms of a universal core surface density ρ₀ × r₀ among different galaxies, which has proven to be challenging for many other DM models (see Section 6).

A future development of this work will involve the study of the DM density profile in the presence of baryons. In fact, the non-minimal coupling to gravity constitutes a natural and effective way to tightly link the DM and baryon properties. On the one hand, this could help to understand puzzling scaling relationships between the DM and the baryonic component, and to characterize the physical processes underlying the emergence of the non-minimal coupling length scale L. On the other hand, it will allow us to probe the effectiveness of our solutions in fitting the measured rotation curves of normal, not just dwarf, rotation-dominated galaxies. In parallel, we plan to extend the static investigation pursued in the present paper to time-dependent conditions, by implementing the non-minimal coupling inside a full N-body numerical simulation. On more general grounds, it would be worth exploring the effect of non-minimal coupling on large, cosmological scales, especially in connection with the dark energy phenomenology.

In conclusion, we have proposed that a non-minimal coupling between matter and gravity could constitute a crucial ingredient toward an improved description of realistic DM structures in a cosmological framework. We very much hope this novel perspective will contribute to shedding light on some of the remaining mysteries concerning the DM component in cosmic structures.

We thank the referee for helpful comments and suggestions. We acknowledge D. Bettoni, L. Danese, M. Nori, and P. Salucci for stimulating discussions and critical reading. This work has been partially supported by PRIN MIUR 2017 prot. 2017-3ML3WW, "Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes." A.L. has taken advantage of the MIUR grant "Finanziamento annuale individuale attivitá base di ricerca" and of the EU H2020-MSCA-ITN-2019 Project 860744 "BiD4BEST: Big Data applications for black hole Evolution Studies." S.L. acknowledges funding from the grant PRIN MIUR 2017 prot. 2017-MB8AEZ.

In this Appendix we discuss the self-gravitating equilibria of non-minimally coupled DM halos when anisotropic conditions apply. These can be included in our treatment by modifying the second of Equations (7) as

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\,\displaystyle \frac{{dp}}{{dr}}+2\beta \,\displaystyle \frac{{\sigma }_{r}^{2}}{r}=-\displaystyle \frac{d{\rm{\Phi }}}{{dr}},\end{eqnarray} \tag{ A1 }$$

where $\beta \equiv 1-\tfrac{{\sigma }_{\theta }^{2}}{{\sigma }_{r}^{2}}$ is the anisotropy parameter of Binney (1978) in terms of the tangential and radial velocity dispersions σ_θ and σ_r , respectively. N-body simulations suggest β(r) to increase from central values β₀ ≲ 0, meaning near isotropy, to outer values β ≳ 0.5, meaning progressive prevalence of radial motions. This overall trend can be physically understood in terms of efficient dynamical relaxation processes toward the inner regions, which tend to enforce closely isotropic conditions, while in the outskirts the infall energy of accreting matter is more easily converted by phase mixing into radial random motions (see Lapi & Cavaliere 2011 for details). Specifically, simulations suggest the effective linear expression (e.g., Hansen & Moore 2006)

$$\begin{eqnarray}&&\beta (r)\simeq {\beta }_{0}+{\beta }_{1}\,[\gamma (r)-{\gamma }_{0}]={\beta }_{0}-{\beta }_{1}{\gamma }_{0}-{\beta }_{1}\,r\,\displaystyle \frac{\rho ^{\prime} }{\rho },\end{eqnarray} \tag{ A2 }$$

in terms of the logarithmic density slope $\gamma (r)\equiv -\tfrac{d\mathrm{log}\rho }{d\mathrm{log}r}$, with γ₀ ≡ γ(0) being a value yet to be determined, β₀ ≲ 0, and β₁ ≈ 0.2.

Adopting Equation (A2) and following the same derivation as in the main text, the fundamental Equation (9) now reads

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{c}\begin{array}{c}\left[1-\eta \,\kappa \,\displaystyle \frac{{\bar{\rho }}^{2-{\rm{\Gamma }}}}{\left({\rm{\Gamma }}-2{\beta }_{1}){\bar{r}}^{\alpha }\right.}\right]\,{\bar{\rho }}^{{\rm{{\prime} }}{\rm{{\prime} }}}+({\rm{\Gamma }}-2)\displaystyle \frac{\bar{\rho }{{\rm{{\prime} }}}^{2}}{\bar{\rho }}\\ \,+\,\displaystyle \frac{\alpha (2{\rm{\Gamma }}-1)+2\,{\rm{\Gamma }}+2({\rm{\Gamma }}-1){\beta }_{0}-2[\alpha +2+({\rm{\Gamma }}-1){\gamma }_{0}]{\beta }_{1}}{{\rm{\Gamma }}-2{\beta }_{1}}\,\displaystyle \frac{{\bar{\rho }}^{{\rm{{\prime} }}}}{\bar{r}}\\ \,+\,\displaystyle \frac{\left[\left.\alpha +1]\,[\alpha +2({\beta }_{0}-{\beta }_{1}{\gamma }_{0}\right)\right]}{{\rm{\Gamma }}-2{\beta }_{1}}\,\displaystyle \frac{\bar{\rho }}{{\bar{r}}^{2}}\\ \,-\,2\,\eta \,\kappa \,\displaystyle \frac{{\bar{\rho }}^{2-{\rm{\Gamma }}}\,{\bar{\rho }}^{{\rm{{\prime} }}}}{\left({\rm{\Gamma }}-2{\beta }_{1}){\bar{r}}^{\alpha +1}\right.}+\kappa \,\displaystyle \frac{{\bar{\rho }}^{3-{\rm{\Gamma }}}}{\left({\rm{\Gamma }}-2{\beta }_{1}){\bar{r}}^{\alpha }\right.}=0.\end{array}\end{array}\end{array}\end{eqnarray} \tag{ A3 }$$

Looking for power-law behaviors $\bar{\rho }\simeq {\bar{r}}^{-\gamma }$, one obtains

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Gamma }}({\rm{\Gamma }}-1)\left[\gamma -\displaystyle \frac{\alpha +2{\beta }_{0}}{{\rm{\Gamma }}}-2\displaystyle \frac{{\beta }_{1}}{{\rm{\Gamma }}}(\gamma -{\gamma }_{0})\right]\\ \ \ \times \,\left[\gamma -\displaystyle \frac{\alpha +1}{{\rm{\Gamma }}-1}\right]-\eta \,\kappa \,\displaystyle \frac{\gamma (\gamma -1)}{{\bar{r}}^{\gamma (2-{\rm{\Gamma }})+\alpha }}=-\displaystyle \frac{\kappa }{{\bar{r}}^{\gamma (2-{\rm{\Gamma }})+\alpha -2}},\end{array}\end{eqnarray} \tag{ A4 }$$

which, remarkably, allows one to self-consistently determine the central slope γ₀.

In fact, for minimally coupled halos (η = 0), the anisotropic solutions feature a modified inner slope ${\gamma }_{0}=\tfrac{\alpha +2{\beta }_{0}}{{\rm{\Gamma }}}$ with respect to the isotropic case, while retaining the same slopes at intermediate radii ${\gamma }_{\alpha }=\tfrac{2-\alpha }{2-{\rm{\Gamma }}}$ and in the outer region ${\gamma }_{\infty }=\tfrac{\alpha +1}{{\rm{\Gamma }}-1};$ as mentioned above, β₀ ≈ 0 so that the changes are minor (if any, β₀ ≲ 0 so that the inner profile is flattened a bit). In addition, the critical solution is characterized by a value ${\alpha }_{\mathrm{crit}}=\tfrac{{\rm{\Gamma }}(5{\rm{\Gamma }}-6)+2{\beta }_{0}({\rm{\Gamma }}-1)({\rm{\Gamma }}-2)}{3{\rm{\Gamma }}-2};$ in particular, ${\alpha }_{\mathrm{crit}}=\tfrac{35-4{\beta }_{0}}{27}$ holds for ${\rm{\Gamma }}=\tfrac{5}{3}$. The corresponding inner, intermediate, and outer slopes read ${\gamma }_{0,\mathrm{crit}}=\tfrac{5{\rm{\Gamma }}-6+2{\beta }_{0}{\rm{\Gamma }}}{3{\rm{\Gamma }}-2}$, ${\gamma }_{\alpha ,\mathrm{crit}}=\tfrac{5{\rm{\Gamma }}-2+2{\beta }_{0}({\rm{\Gamma }}-1)}{3{\rm{\Gamma }}-2}$, and ${\gamma }_{\infty ,\mathrm{crit}}=\tfrac{5{\rm{\Gamma }}+2+2{\beta }_{0}({\rm{\Gamma }}-2)}{3{\rm{\Gamma }}-2}$, respectively.

For non-minimally coupled halos, an inner core with γ₀ ≈ 0 is enforced anyway by the second term in Equation (A4), so that the variations with respect to the isotropic case are minor and limited to the outermost regions. In Figure A1 we show how the non-minimally coupled density profile with η = −0.05, ${\rm{\Gamma }}=\tfrac{5}{3}$, and $\alpha =\tfrac{35-4{\beta }_{0}}{27}={\alpha }_{\mathrm{crit}}$ is affected by anisotropies. For realistic values −0.1 ≲ β₀ ≲ 0.1 and β₀ ≈ 0.2, the profile is marginally affected in the inner region and at intermediate radii, while it tends to extend toward slightly larger radii (i.e., the cutoff moves outward) due to the progressive prevalence of radial anisotropy in the halo outskirts; such an effect, though minor, is more pronounced for larger (more positive) β₀.

Zoom In Zoom Out Reset image size