Self-gravitating Equilibria of Non-minimally Coupled Dark Matter Halos
DDraft version February 9, 2021
Typeset using L A TEX default style in AASTeX63
Self-gravitating Equilibria of Non-minimally Coupled Dark Matter Halos
Giovanni Gandolfi,
1, 2
Andrea Lapi,
1, 2, 3, 4 and Stefano Liberati
1, 2, 3 SISSA, Via Bonomea 265, 34136 Trieste, Italy IFPU - Institute for fundamental physics of the Universe, Via Beirut 2, 34014 Trieste, Italy INFN-Sezione di Trieste, via Valerio 2, 34127 Trieste, Italy INAF-Osservatorio Astronomico di Trieste, via Tiepolo 11, 34131 Trieste, Italy
ABSTRACTWe investigate self-gravitating equilibria of halos constituted by dark matter (DM) non-minimallycoupled to gravity. In particular, we consider a theoretically motivated non-minimal coupling whichmay arise when the averaging/coherence length L associated to the fluid description of the DM collec-tive behavior is comparable to the local curvature scale. In the Newtonian limit, such a non-minimalcoupling amounts to a modification of the Poisson equation by a term L ∇ ρ proportional to the Lapla-cian of the DM density ρ itself. We further adopt a general power-law equation of state p ∝ ρ Γ r α relating the DM dynamical pressure p to density ρ and radius r , as expected by phase-space densitystratification during the gravitational assembly of halos in a cosmological context. We confirm previ-ous findings that, in absence of the non-minimal coupling, the resulting density ρ ( r ) features a steepcentral 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. Most importantly, we findthat the non-minimal coupling causes the density distribution to develop an inner core and a shapeclosely following, out to several core scale radii, the Burkert profile. In fact, we highlight that theresulting mass distributions can fit, with an accuracy comparable to the Burkert’s one, the co-addedrotation curves of dwarf, DM-dominated galaxies. Finally, we show that non-minimally coupled DMhalos are consistent with the observed scaling relation between the core radius r and core density ρ ,in terms of an universal core surface density ρ × r among different galaxies. Keywords:
Cosmology (343) - Dark matter (353) - Non-standard theories of gravity (1118) INTRODUCTION N − body, dark matter (DM)-only simulations in the standard ΛCDM cosmology suggest an almost universal shapeof the density distributions ρ ( r ) within virialized DM halos of different masses and redshifts. This fact has beenestablished since the seminal work by Navarro et al. (1996), who proposed an empirical fitting formula (the so-calledNFW profile) for the DM density run ρ ( r ) ∝ r ( r + r s ) that still nowadays constitute the standard lore (but see Navarroet 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 curvesin many dwarf, DM-dominated galaxies, that 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 etal. 2015; Genzel et al. 2020; for a review and further references, see Bullock & Boylan-Kolchin 2017). The observeddensity distribution is usually described with the phenomenological Burkert (1995) profile ρ ( r ) ∝ r + r ) ( r + r ) .A non-trivial point to stress is that, besides the flat shape of the inner density run, the measurements indicate anearly universal relationship between the core density and radius, in terms of a constant value for the “core surfacedensity” ρ × r among different galaxies (see Salucci & Burkert 2000; Donato et al. 2009; Gentile et al. 2009; Burkert2015; Kormendy & Freeman 2016). In addition, another puzzling aspect of the DM phenomenology is the existenceof tight scaling laws with baryonic quantities. These include the baryonic Tully-Fisher relation (Tully & Fisher 1977;McGaugh 2012), the core radius vs. disc 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). a r X i v : . [ a s t r o - ph . C O ] F e b Gandolfi, Lapi, Liberati
Admittedly, through the last 20 years the cusp-core controversy has flamed the scientific debate. One class of solutionsinvokes physical processes that can cause violent fluctuations in the inner gravitational potential and/or transfer ofenergy and angular momentum from the baryons to DM, thus possibly erasing the central density cusp; some of themost explored possibilities include dynamical friction (see El-Zant et al. 2001, 2016; Tonini et al. 2006; Romano-Diazet 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 simulationsincluding baryonic physics have shown that a variety of halo responses are originated depending on the stellar/halomasses (e.g., Freundlich et al. 2020b).An alternative solution, perhaps more fascinating, is to abandon the cold DM hypothesis and look at non-standardparticle candidates (see review by Salucci 2019). A few examples include: self-interacting DM particles with cross-section ∼ g − , which can create a core as they are heated up via elastic two-body collisions and evacuatedfrom 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 ultra-light particles (likely axions) with masses ∼ − eV, for whichthe core stems from the equilibrium between quantum pressure and gravity (e.g., Hu et al. 2000; Bohmer & Harko2007; Schive et al. 2014a,b; Harko 2014; Hui et al. 2017; Bernal et al. 2018); warm DM, made of fermionic particleswith masses of order a few keV (e.g., sterile neutrinos), that can originate cores where gravity is counterbalanced byquantum 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-cuspproblem envisages that no DM is present and tries to explain the galactic dynamics via a modification of gravity, likein the MOND phenomenological approach (see Milgrom 1983, 2009; Bruneton & Esposito-Farese 2007; Bekenstein2004, 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 subjectto 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 collisionlessDM phenomenology on large cosmological scales, and at the same time to introduce MOND-like behavior in galaxiesby attributing to DM some non-negligible coupling with gravity. The words “non-minimal” simply mean that such acoupling is embodied in the action via an additional interaction term with respect to standard scalar-tensor gravitytheories. We stress that the non-minimal coupling is not necessarily a fundamental feature of the DM particles butmight develop dynamically when the averaging/coherence length L associated to the the fluid description of the mattercollective behavior is comparable to the local curvature scale. We specifically consider and theoretically justify a formof 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 withthe local gravitational curvature. Note that extensions of gravity including non-minimal couplings have been alsoinvestigated on cosmological scales, since under certain conditions they can mimic and explain the properties of thedark 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 densitydistribution closely following the Burkert profile out to several core scale radii. Moreover, we show that these canfit very well the measured rotation curves of dwarf, DM-dominated galaxies, and are consistent with the observeduniversal core surface density. In more detail, the plan of the paper is the following: in Sect. 2 we theoreticallymotivate the adopted non-minimally coupled DM framework; in Sect. 3 we introduce the basic formalism to describeself-gravitating equilibria of non-minimally coupled DM halos; in Sect. 3.1 we describe the effective equation of statefor the DM fluid; in Sect. 3.2 we derive the fundamental equation ruling the DM density profile and study its solutionspace; in Sect. 4 we compare our non-minimally coupled solutions with the density profiles classically adopted inthe literature to fit simulations and/or observations; in Sect. 5 we provide a first glimpse on how our non-minimallycoupled DM mass distributions can fit the measured rotation curves of dwarf, DM-dominated galaxies; in Sect. 6 weshow that non-minimally coupled DM halos are indeed consistent with the observed universal behavior of the coresurface density; finally, in Sect. 7 we summarize our findings and highlight future prospects.Throughout this work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration 2020) with roundedparameter values: matter density Ω M = 0 .
3, dark energy density Ω Λ = 0 .
7, baryon density Ω b = 0 .
05, and Hubbleconstant H = 100 h km s − Mpc − with h = 0 .
7; unless otherwise specified, G ≈ . × − cm g − s − indicatesthe standard gravitational (Newton) constant. quilibria of non-minimally coupled DM halos A THEORETICAL FRAMEWORK FOR NON-MINIMALLY COUPLED DMThe motivation for introducing a non-minimal coupling between the DM matter field and gravity is twofold. On thetheoretical side, it is allowed by the Einstein equivalence principle (Di Casola et al. 2015), and it might be requiredfor the renormalizability of quantum field theories in curved spacetimes (e.g., Sonego & Faraoni 1993; Bruneton etal. 2009). On the observational side, it may help in explaining the striking relationships between DM and baryonsrecalled in Sect. 1 (e.g., baryonic Tully-Fisher, core radius vs. disc scale-length, acceleration relations, etc.) thatare not trivially understood via galaxy formation processes. In fact, such relationships may indicate a very specialDM-baryon interaction. However, on the one hand, evidence of any direct coupling is missing. On the other hand,modified-gravity models like MOND seem, at least in galactic environments, able to describe the baryon dynamics asthat of 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 experiencedby baryons may not coincide with the gravitational one, but it may be also determined by the properties of theDM field. Indeed, it was shown by Bekenstein (1993) that the most generic transformation between physical andgravitational metric preserving causality and the weak equivalence principle can be of the general form (cid:101) g µν = e ϕ [ A ( X ) g µν + B ( X ) ∇ µ ϕ ∇ ν ϕ ] , (1)where A and B are functions to be specified, ϕ is an extra scalar field and X = − g µν ∇ µ ϕ ∇ ∂ ν ϕ ; such a relationbetween metrics is called a disformal one. Following the aforementioned idea, we can now ask which kind of interactionmight be reduced to an effective coupling via a disformal metric of the above form with the scalar field playing therole of DM.To this purpose, one can start from a general action of the form S = S EH [ g µν ] + S bar [ g µν , ψ ] + S DM [ g µν , ϕ ] + S int [ g µν , ψ, ϕ ] , (2)where the terms on the right hand side identify respectively the Einstein-Hilbert (standard general relativity) action S EH = c πG (cid:82) d x √− g R in terms of the Ricci scalar R , the baryonic S bar and DM S DM actions, and the interactionone S int . Here the scalar fields ψ and ϕ are thought 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 Eq. (1) if S bar [ g µν , ψ ] + S int [ g µν , ψ, ϕ ] ≈ S bar [ g µν + h µν , ψ ] with h µν ∝ ∇ µ ϕ ∇ ν ϕ . Up to order O ( h ), this is originated by aninteraction term of the form S int [ g µν , ψ, ϕ ] ∝ (cid:82) d x √− g T µν bar ∇ µ ϕ ∇ ν ϕ , where T µν bar is the baryonic matter stress-energytensor.Noticeably, if now we express the full action in terms of the physical metric (cid:101) g µν ≡ g µν + h µν (which is equivalent tochoose a frame in which baryons follow the geodesics of this metric, i.e., the Jordan frame), then it turns out that theDM field gets non-minimally coupled to gravity as S = S EH [ (cid:101) g µν ] + S bar [ (cid:101) g µν , ψ ] + S DM [ (cid:101) g µν , ϕ ] + (cid:15)L (cid:90) d x (cid:112) − ˜ g (cid:101) G µν ∇ µ ϕ ∇ ν ϕ , (3)where L is a coupling length that must be present for dimensional consistency, (cid:15) = ± (cid:101) G µν is the Einstein tensor expressed in terms of the physical metric (cid:101) g µν , 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 physicaldimensions (still leading to second order field equations) would be proportional to the Ricci scalar as X (cid:101) R , which ishowever equivalent to that appearing in Eq. (3) modulo a surface term (see Bettoni & Liberati 2013). Second, thecoupling term (cid:101) G µν ∇ µ ϕ ∇ ν ϕ 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 weare not proposing a fundamental theory of modified gravity, but just entailing the possibility that DM in galactic halosdynamically develops a non-minimal coupling with the metric characterized by an effective length-scale L . Third, itcan 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 Eqs. 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 (atleast in dwarf, DM-dominated galaxies); from this point of view, the coupling length L cannot be straightforwardlyinterpreted in terms of the acceleration parameter a appearing in the MONDian dynamics. Gandolfi, Lapi, Liberati
In order to take the Newtonian limit of the above theory it is convenient to adopt fluid variables. In the case of acomplex scalar field this is easily achieved by adopting the standard Madelung representation ϕ ∼ √ ρe iθ . However,a fluid limit is possible also for a non-minimally coupled real scalar field (see Bettoni et al. 2012). In the end it canbe shown that the non-minimal coupling considered above leads to a modified Poisson equation of the form (see e.g.,Bettoni et al. 2014) ∇ Φ = 4 πG [( ρ + ρ bar ) − (cid:15) L ∇ ρ ] , (4)where Φ is the Newtonian potential, and ρ bar and ρ are the baryon and DM mass densities. This modified Poissonequation implies that the source for gravity is not just the total matter density in itself, but also the DM inhomogeneitiesor spatial variations . In the present paper we will focus mainly on dwarf, DM-dominated galaxies and thus we willneglect the baryonic component ρ bar hereafter.While the discussed non-minimal coupling could be associated to some modified theory of gravity, in the extantliterature at least two main mechanisms have been contemplated to produce dynamically a non-minimal couplingcharacterised by a length-scale L : either it can emerge from some collective behavior of the DM particles associated toa coherence length (for example via a Bose-Einstein condensation mechanism; see Bettoni et al. 2011 for an extendeddiscussion), or it appears through an averaging procedure associated to the fluid description of matter. Although ingeneral 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 Eq. (4) for the self-gravitatingequilibria of DM halos, without an a priori prejudice about the origin and scale of L . SELF-GRAVITATING EQUILIBRIA OF DM HALOSThe self-gravitating equilibria of DM halos can be specified in the fluid approximation (Teyssier et al. 1997; Sub-ramanian et al. 2000; Lapi & Cavaliere 2011; Nadler et al. 2017) via the continuity, Euler (also called Jeans in thiscontext), and Poisson equations ∂ t ρ + ∇ · ( ρ v ) = 0 ,∂ t v + ( v · ∇ ) v + 1 ρ ∇ p = −∇ Φ , ∇ Φ = 4 πG ( ρ − (cid:15) L ∇ ρ ) . (5)Here v is the bulk velocity, p = ρ σ r is the pressure dynamically generated by the random motions (and specified interms of a radial velocity dispersion σ r or more generally of an anisotropic stress tensor σ ij ) of the DM particles inapproximate virial equilibrium within the gravitational potential well . The last term on the right hand side of thePoisson equation represents the non-minimal coupling with length-scale L and polarity (cid:15) as discussed in Sect. 2; wewill see that for the purpose of originating physically acceptable DM density distributions, a negative polarity (cid:15) = − L will turn out to be closely related to the DM core radius. Note that non-minimal coupling termsdo not appear in the Euler equation since they are found to be sub-leading in the non-relativistic limit (i.e., expansionin 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 presentcontext, we focus on the EOS originated by the gravitational assembly of DM halos via accretion and mergers from thecosmic web. This EOS stems from the progressive stratification of the DM pseudo-entropy K ≡ p/ρ / , or equivalentlyof the coarse-grained phase-space density ρ/σ r ∝ K − / , in terms of a simple power-law profile K ( r ) ∝ r α . Althoughthe 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; Nolting et al. While for simplicity we have written our formulas in terms of a real scalar field, the aforementioned derivative non-minimal coupling canbe easily generalized to a complex scalar field as (cid:101) G µν ∇ µ ϕ ∇ ν ϕ † . Remarkably, the same formal modification can be obtained starting directly from a cosmological fluid description (see Bettoni & Liberati2015) albeit in this case there is no fundamental reason to couple the fluid directly to the Einstein tensor (and indeed it is necessary tocouple separately the fluid to the Ricci scalar and/or Ricci tensor to get the same kind of modification, see Bettoni et al. 2014). It is alsoworth reporting that the same form of modified Poisson equation can be derived in Born-Infield gravity (e.g., Beltran Jimenez et al. 2018). Note that p is not to be confused with the relativistic pressure p ≈ quilibria of non-minimally coupled DM halos α ≈ . − .
3. Such values are indeed expected on the basis of simpleself-similar arguments (see Bertschinger 1985; Lapi & Cavaliere al. 2009a, 2011; Nadler et al. 2017) and also broadlyconsistent with observations (see Lapi & Cavaliere 2009b; Chae 2014; Munari et al. 2014). On this basis, but to keepsome degree of generality , we adopt the EOS parametrization p ( ρ, r ) = λ ρ Γ r α = ρ α σ α (cid:32) ρρ α (cid:33) Γ (cid:32) rr α (cid:33) α , (6)where r α is a reference radius and we have defined ρ α ≡ ρ ( r α ) and σ α ≡ σ r ( r α ); in the following we adopt as fiducialvalues Γ ≈ and α ≈ . v = 0), spherically symmetric, and isotropic conditions (see Appendix for a generalization) the relevantEqs. (5) become ρ d p d r = − dΦd r , r dd r (cid:32) r dΦd r (cid:33) = 4 πG (cid:34) ρ − (cid:15) L r dd r (cid:32) r d ρ d r (cid:33)(cid:35) , (7)supplemented with the trivial mass conservation constraint M ( t ) =const, where M ≡ π (cid:90) ∞ d r r ρ ( r ) , (8)is the total mass M of the DM halo.3.2. The fundamental equation and its solutions
It is convenient to introduce normalized variables ¯ r ≡ r/r α , ¯ ρ ≡ ρ/ρ α and define the quantities κ ≡ πG ρ α r α /σ α and η ≡ (cid:15) L /r α . To understand the physical meaning of κ , one can choose r α = r max to be the point at which thecircular velocity v c ( r ) ≡ G M ( < r ) /r peaks at a value v c ( r max ) = 4 πρ ( r max ) r , so that κ = v c ( r max ) /σ ( r max ) isseen to compare the estimate σ ( r max ) for the random kinetic energy with that v c ( r max ) for the gravitational potential.Eliminating dΦ / d r from Eqs. (7) and using the EOS Eq. (6) yields the following fundamental equation for the density¯ ρ (cid:48)(cid:48) + (Γ −
2) ¯ ρ (cid:48) ¯ ρ + α (2Γ −
1) + 2 ΓΓ ¯ ρ (cid:48) ¯ r + α ( α + 1)Γ ¯ ρ ¯ r − η κ ¯ ρ − Γ ¯ ρ (cid:48) Γ ¯ r α +1 + κ ¯ ρ − Γ Γ ¯ r α − η κ ¯ ρ − Γ Γ ¯ r α = 0 , (9)while the mass conservation constraint now reads M = ρ α r α f M in terms of the shape factor f M = 4 π (cid:82) ∞ d¯ r ¯ r ¯ ρ (¯ r ).The solution space of such an equation is amazingly rich, and for Γ = and η = 0 it has been quite extensivelystudied in the literature to describe the radial structure of standard ΛCDM halos (see Williams et al. 2004; Hansen2004; Austin et al. 2005; Dehnen & McLaughlin 2005; Lapi & Cavaliere 2009a). In the following we provide thegeneralization 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 powerlaw behaviors ¯ ρ (cid:39) ¯ r − γ .Substituting in the fundamental equation yieldsΓ (Γ − (cid:16) γ − α Γ (cid:17) (cid:18) γ − α + 1Γ − (cid:19) − η κ γ ( γ − r γ (2 − Γ)+ α = − κ ¯ r γ (2 − Γ)+ α − . (10)Focusing first on the minimally coupled case with η = 0, it is evident that trivial powerlaw solutions with slope γ = γ α ≡ − α − Γ are admitted, implying κ = κ PL ≡ (Γ − α ) ( α +4 − − Γ) . These values are actually not physically acceptable For example, some authors (e.g., Schmidt et al. 2008; Hansen et al. 2010) have claimed that it is the quantity ρ/σ εr ∝ r − ζ with ε (cid:46) α = ζ and Γ = 1 + ε . Gandolfi, Lapi, Liberati
Figure 1.
Example aimed at illustrating the shooting technique adopted to solve the fundamental Eq. (9) and to find thephysically acceptable solutions. Radial coordinate and density profile are normalized to the reference radius r α where thelogarithmic slope γ ≡ − d log ρ d log r of the profile is γ α = − α − Γ . The differential equation is integrated inward of r α with boundarycondition ρ ( r α ) = ρ α and ρ (cid:48) ( r α ) = − γ α ρ α r α . Red lines show shot solutions for different values of the constant κ appearing inEq. (9), while the black line is the only physical profile for κ = κ Γ ,α , see text for details (for κ > κ Γ ,α shot solutions are belowthe physical profile, while for κ < κ Γ ,α are above). The two dashed lines indicate the asymptotic slope γ = α Γ in the innerregion and the slope γ α at r = r α . In this example we specifically considered the minimally-coupled case η = 0, and EOSparameters Γ = and α = α crit = , yielding κ Γ ,α ≈ . γ = and γ α = . at small and large radii since the gravitational force and mass would diverge, but provide the behavior of any solutionat intermediate radii. In fact, to find general solutions it is numerically convenient to choose r α as the reference radiuswhere the logarithmic slope of the density is − γ α and then integrate inward and outward. In terms of the normalizedvariables, this corresponds to set the boundary conditions as ¯ ρ (1) = 1 and ¯ ρ (cid:48) (1) = − γ α . 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 value κ Γ ,α the profile asymptotes for ¯ r (cid:28) ρ ∝ ¯ r − γ with γ ≡ α Γ , consistently with the power counting in Eq. (10). For κ > κ Γ ,α the profile has wiggles (changeof 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., automaticallysolving inward the differential equation with different slopes − γ α at r α until the desired, physical inner asymptoticbehavior is found), as represented in Fig. 1.Fig. 2 illustrates the resulting full density profiles, for three different values of α . The outer behavior of the profileturns out to be physical only for α ≤ α crit ≡ Γ(5Γ − − . For α > α crit the outer slope is too flat, to imply a divergingmass. For α < α crit the solution ¯ ρ ∝ ¯ r − γ ∞ attains a slope γ ∞ ≡ α +1Γ − at a finite large radius before an outercutoff. For α = α crit the cutoff is pushed to infinity and the slope γ ∞ , crit ≡ − is attained only asymptotically for¯ r → ∞ , consistently with the power-counting from Eq. (10); correspondingly, the intermediate and inner slopes read γ α, crit ≡ − − and γ , crit ≡ − − , respectively. We illustrate these different behaviors in Fig. 2 for the fiducial valueΓ = . In such a case, the powerlaw solutions have γ α = 6 − α and κ PL = 6 (5 − α ) ( α − quilibria of non-minimally coupled DM halos Figure 2.
Examples aimed at showing the outer behavior of the density profiles from solving the fundamental Eq. (9). Theillustrated profiles refer to the minimally coupled case η = 0 with EOS parameter Γ = and three different values of α : thered line corresponds to α = = α crit for which the density asymptotes to a slope γ ∞ = α +1Γ − ≈ ; the blue line refers to α = 1 . < α crit for which the slope γ ∞ is attained at a finite radius before a cutoff; the green line refers to α = 1 . > α crit for which the outer slope is unphysical since it would lead to a diverging mass. The inset reports the dependence on α of theconstant κ Γ ,α for the full solutions (colored dots and solid black line), and for the pure power-law solutions (dashed black line),see Sect. 3.2 for details. solutions feature an inner slope γ = α and an outer slope γ ∞ = α )2 ; the values of κ α, Γ are reported as a functionof α in the inset of Fig. 2. For α = α crit = one gets γ , crit ≡ , γ α, crit = , γ ∞ , crit = , and κ α, Γ = ≈ . η ≡ (cid:15) L r α , the solution space changes appreciably. First of all, (cid:15) (hence η )must be negative, otherwise the inner profile diverges at a finite radius. Then for any negative value of η , there is againan 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 Eq. (10) since now the second term on the l.h.s. dominatesthe behavior for small ¯ r (cid:28)
1. Other solutions with κ smaller or larger than κ Γ ,α,η are not physically acceptable sincethey 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 R setting theeffective halo boundary, which is smaller for more negative values of η . We illustrate the physically acceptable profilesfor different η and the related values of the constant κ Γ ,α,η in Fig. 3. The corresponding distributions of mass, circularvelocity and velocity dispersion for a few values of the non-minimal coupling η are also illustrated in Fig. 4. COMPARISON WITH LITERATURE PROFILESWe now compare the shape of our physical solutions to some classic literature density profiles, characterized bydifferent analytic expressions and numbers of parameters, that are commonly adopted to fit simulations and/or obser-vations (see also Freundlich et al. 2020b, their Fig. 15 for a comprehensive account). To this purpose, it is convenientto use a radial coordinate ˆ r ≡ r/r − normalized to the radius r − where the logarithmic density slope d log ρ d log r = − ρ ≡ ρ/ρ ( r − ) accordingly. Specifically, we will consider the following density profiles. • αβγ profiles Gandolfi, Lapi, Liberati
Figure 3.
Density profiles (EOS parameters Γ = and α = α crit = ) for different values of the non-minimal coupling η = 0(grey), − .
001 (orange), − .
005 (brown), − .
01 (magenta), − .
025 (red), − .
05 (green), − .
075 (blue), − . κ Γ ,α,η for the physical solutions as a function of η (colored dots and solid line),and for the η = 0 case (dashed black line). The αβγ profiles (Zhao 1996; see also Widrow 2000) feature the shapeˆ ρ (ˆ r ) = ˆ r − τ (cid:32) w w ˆ r ω (cid:33) ξ , (11)where the three parameters τ , ω , ξ describe respectively the central slope, the middle curvature and the outerdecline of the density run, while w ≡ − − τ (2 − τ − ω ξ ) . Familiar empirical profiles are recovered for specific valuesof the triplet ( τ, ω, ξ ): e.g., Plummer’s profile corresponds to (0 , , / , , N − body simulations in theΛCDM model, is obtained for the parameter triple (1 , , τ, , − τ ) apply (Mamon et al. 2019; for a more complex cored versionsee also Read et al. 2016). It is worth mentioning that recently Freundlich et al. (2020b) have considered a αβγ profile (referred also as Zhao-Dekel model) with parameters ( τ, , − τ ), that can provide good fits to thedensity profiles from both N − body, DM-only and hydro simulations including baryonic effects. • Sersic-Einasto profileThe Sersic-Einasto profile (see An & Zhao 2013) is defined asˆ ρ (ˆ r ) = ˆ r − τ e − u (ˆ r ω − , (12)where τ is the inner density slope, ω is a shape parameter and u ≡ (2 − τ ) ω . The classic cored Einasto shape(Sersic 1963; Einasto 1965; Prugniel et al. 1997; Graham et al. 2006; also Lazar et al. 2020 for a more complexanalytical expression) is recovered for τ = 0. N − body simulations in the standard ΛCDM model are usually quilibria of non-minimally coupled DM halos Figure 4.
Profiles of mass (top left panel), circular velocity (top right panel) and velocity dispersion (bottom left panel)corresponding to a few of the density profiles illustrated in the previous Figure (EOS parameters Γ = and α = α crit = ),with η = 0 (black), − .
001 (orange), − .
01 (magenta), − . well described by the parameter values ω ≈ . − . τ ≈ τ (cid:46) . • Soliton profileThe soliton profile features the shape (see Schive et al. 2014a,b)ˆ ρ (ˆ r ) = (cid:32) ω/ ˆ r c ω ˆ r / ˆ r c (cid:33) , (13)where the normalized core radius reads ˆ r c = √ ω . In numerical simulations of non-interacting Bose-Einsteincondensate DM (alias fuzzy DM), halos are described by a combination of this solitonic profile with ω ≈ . • BEC profileThe profile followed by an interacting Bose-Einstein condensate (BEC; see Bohmer & Harko 2007; Harko 2014)in the Thomas-Fermi limit reads ˆ ρ (ˆ r ) = 1ˆ r sin( π ˆ r/ ˆ R )sin( π/ ˆ R ) , (14)where the normalized halo boundary is defined by the equality π/ ˆ R + tan( π/ ˆ R ) = 0. Incidentally, note that sucha profile actually corresponds to the solution of our Eq. (10) for Γ = 2 and α = η = 0. • Burkert profile0
Gandolfi, Lapi, Liberati
Figure 5.
Comparison of non-minimally coupled density distributions with classic literature profiles (see Sect. 4 for details):NFW (red line), Einasto (magenta line), soliton (blue line), interacting Bose-Einstein condensate (cyan line), and Burkert (greenline). The dot black line refers to η = 0, the solid black line is for η = − .
05 and the grey-shaded area illustrates the regioncovered by η in the range from − . − .
01 (upper envelope). The radial coordinate and the density profileshave been normalized to the radius r − where the logarithmic density slope d log ρ d log r = − The Burkert profile can be written as (see Burkert 1995; Salucci & Burkert 2000)ˆ ρ (ˆ r ) = (1 + ˆ r c ) (1 + ˆ r c )(ˆ r + ˆ r c ) (ˆ r + ˆ r c ) , (15)where ˆ r c is the normalized core radius, defined by the nonlinear algebraic equation 2 ˆ r c + ˆ r c − and α = α crit = are adopted) to a few ofthe above profiles. First it is evident that our density run for η = 0 describes quite well the NFW and Einasto profilescommonly used to fit N − body simulations in the standard ΛCDM cosmology. For − . < η < − .
01 our solutionsdevelop a core and the shape out to few/several r − is remarkably close to the Burkert profile, commonly exploitedto fit observations of dwarf galaxies. At larger distances there is progressive deviation from the Burkert profile, sincethe latter has been designed to have a limiting slope close to the NFW one, while our profiles get truncated at a finiteradius R . However, for an appreciable range of η values this is not a concern since the truncation occurs at radii muchbeyond r − that are scantily (if at all) probed by observations (see also Sect. 5). We stress that it is a remarkableproperty of the non-minimally coupled solutions to reproduce the Burkert shape not only in the inner region, but alsoover an extended radial range outward of r − . For comparison, other models such as those based on Bose-Einsteincondensate DM predict a cored density profile, but the deviation from the Burkert shape are appreciable, in terms ofa prominent cutoff or steep decline, just outside the core. Note that our solutions for η = − .
05 can be reasonablydescribed by a αβγ model (cf. Eq. 11) with parameters ( τ, ω, ξ ) = (0 , , ), or by a Sersic-Einasto model (cf. Eq. 12)with parameters ( τ, ω ) = (0 , ) up to several r − before the final cutoff. quilibria of non-minimally coupled DM halos Figure 6.
Comparison of non-minimally coupled halo mass distributions with observed dwarf galaxy rotation curves. Datapoints (red stars; Lapi et al. 2018) refer to the co-added rotation curves of about 20 dwarf galaxies with I-band magnitude M I (cid:38) − .
5, extracted from the original sample by Persic & Salucci (1996). Solid line illustrates the fit via our physical solutionwith non-minimal coupling parameter η ≈ − .
05, while the shaded grey area show the effect of changing η in the range from − . − .
01. For comparison, dashed line illustrates the fit with the Burkert profile, and dotted line that with the NFWprofile. 5.
COMPARISON WITH MEASURED ROTATION CURVESWe now compare the non-minimally coupled mass distributions with observed galaxy rotation curves. To avoiddealing with baryons (which are not included in our treatment), we focus on dwarf, strongly DM-dominated galaxies.Specifically, we exploit the co-added rotation curve built by Lapi et al. (2018) based on the high-quality measurementsof about 20 dwarf galaxies with I − band magnitude M I (cid:38) − . R e (cid:46) . M (cid:63) (cid:46) M (cid:12) and halo mass M (cid:46) M (cid:12) .The co-added 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 profile, andwith our physical solutions for different values of the non-minimal coupling parameter η ; the outcomes are illustratedin Fig. 6.As it is well known, the NFW fit (equivalent to our solution with η = 0) struggles to fit the measured dwarf galaxyrotation curves, yielding a reduced χ ≈
4. On the other hand, the Burkert profile performs much better, providinga good fit with a reduced χ ≈ .
58. Our non-minimally coupled solution for η ≈ − .
05 provides a fit of qualitycomparable to the Burkert one, yielding practically the same reduced χ ≈ .
61. The current data are compatiblewithin 3 σ with any value of η in the range from − . − .
01. Accurate determinations of the co-added rotation curveout to 10 kpc or beyond would be necessary to determine the non-minimal coupling parameter η to a good level ofprecision. Interestingly, one can also see from Fig. 6 that there is a tendency to favor values of η slightly less negativethan − .
05, which is actually what is to be expected in dwarf galaxies with halo masses M (cid:46) M (cid:12) on the basis ofuniversal scaling arguments (see Sect. 6 and in particular Eq. 19).2 Gandolfi, Lapi, Liberati
Figure 7.
Core surface density ρ × r as a function of the core radius r . Data points from Burkert (2020) are well fitted bythe universal value 75 +55 − M (cid:12) pc − (blue shaded area). Our prediction for non-minimally coupled DM halos is illustrated bythe black lines for three different halo formation redshifts z ≈ . . From the fit with η ≈ − .
05 (or equivalently from the Burkert one) we derive a best fit value of the core radius around r ≈ . ± . h µν characterizing the interaction term in the action of Eq. (3)can lead to a MONDian phenomenology on galactic scales, that 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 dynamicsin its Newtonian limit; from this perspective, the rather good performance of our simple non-minimally coupled modelin reproducing the measured rotation curves is valuable and far from trivial. UNIVERSAL CORE SURFACE DENSITYAs mentioned in Sect. 1, it has been well established observationally (see Salucci & Burkert 2000; Burkert 2015)that, at least for dwarf galaxies with halo masses M (cid:46) M (cid:12) , the product of the core density and core radius(i.e., a sort of core surface density) is an approximately universal constant with values ρ r ≈ +55 − M (cid:12) pc − amongdifferent galaxies, see Fig. 7. This somewhat unexpected property poses a serious challenge to any theoretical modelof core formation (e.g., Deng et al. 2018; Burkert 2020). We now aim to show that non-minimally coupled DM halosare instead consistent with such a remarkable scaling law.First, adopting for definiteness Γ = and α = α crit = , we compute the physical solutions of Eq. (9) subject tothe boundary conditions ρ ( r α ) = ρ α and ρ (cid:48) ( r α ) = − γ α ρ α r α , for several values of the non-minimal coupling parameter η .From the so obtained normalized profiles we fit as a function of | η | = L r α in the range − . (cid:46) η (cid:46) − .
01 the followingrelations involving the normalized boundary radius R , the core radius r , the core density ρ and the mass shape quilibria of non-minimally coupled DM halos f M defined below Eq. (9): R r α (cid:39) . | η | − . , f M (cid:39) . | η | − . , ρ ρ α (cid:39) . | η | − . , r r α (cid:39) . | η | . . (16)We then combine the above scaling with the expression for the total mass M = f M ρ α r α and with the definition ofthe virial radius R = (3 M / π ∆ vir ρ c E z ) / ≈ E − / z ( M / M (cid:12) ) / ; here ρ c ≈ . × h M (cid:12) Mpc − is thecritical density, E z ≡ Ω M (1+ z ) +Ω Λ takes into account the formation redshift z of the halo, and ∆ vir is the nonlinearthreshold for virialization, with values around 100 at z ≈ z (cid:38)
1. We eventually derive r (cid:39) . L , r α (cid:39) . L . R . ρ (cid:39) . vir ρ c E z (cid:18) R L (cid:19) . = ρ α (cid:18) R L (cid:19) . ; (17)remarkably, the core radius r turns out to be proportional, with a coefficient of order 1, to the non-minimal couplinglength-scale L .To proceed further we need a relation between the core radius r (or L ) and the halo mass M ; this is thought notto be of fundamental nature but rather to stem from two other relationships involving the baryonic mass: (i) therelation between the stellar (disc) mass and the halo mass (e.g., Moster et al. 2013), which is known to be originatedby baryonic processes related to galaxy formation; (ii) the relation between the core radius r and the disk scalelength(in turn related to the stellar mass; e.g., Donato et al 2004), which is instead still not completely understood. Inthe present pilot study we are not including the baryonic component and we cannot infer the r − M relation fromfirst principles; thus we will adopt the outcome r ≈ . M / M (cid:12) ) . kpc from the dynamical modeling study bySalucci et al. (2007), and see what this implies for the core surface density. Specifically, from Eqs. (17) we obtainΣ ≡ ρ × r ≈ (cid:18) ∆ vir (cid:19) E . z M (cid:12) pc − (18)independent of the halo mass and/or core radius, and only weakly dependent on formation redshift. In Fig. 7 we reportthe above for three values of the formation redshift z ≈
0, 0 .
5, 1 . r − M relation we have assumed. For example, as pointed out by Burkert (2020) fuzzy DM can reproducethe r − M relation, albeit with some (uncertain) hypothesis on core formation redshift. However, such a model isconsiderably out of track as to the core surface density scaling, since it robustly predicts ρ ∝ r − . The same issueconcerns many other DM models inspired by particle physics, as extensively discussed, e.g., by Deng et al. (2018).As an aside, from Eqs. (7) and the adopted r − M relation we can also derive other three interesting scaling laws.First, the dependence of η on halo mass reads | η | ≈ . (cid:18) M M (cid:12) (cid:19) . E . z ; (19)this confirms that values of the non-minimal coupling in the range η = − . − .
01 cover the typical mass rangeof dwarf galaxies. We stress that this dependence of η on halo mass/formation redshift will induce slightly differentshapes in the profiles, implying a weak violation of self-similarity. Though challenging, it will be interesting to lookfor such behaviors in real data (see Sect. 5). Second, it may be interesting to derive the dependence on the coupling,hence on mass, of the inner logarithmic slope γ . measured at a reference radius of r ≈ . r − ∼ a few percent of R .We get the scaling | γ . | (cid:39) . | η | − . that after Eq. (19) translates into a mass-dependence | γ . | (cid:39) . (cid:18) M M (cid:12) (cid:19) − . E − . z . (20)Thus there is a slight tendency for less massive halos to have flatter profile at a fixed radius in the inner region (notethat asymptotically at the center all the non-minimally coupled DM density profiles are flat). In a future work, itwould 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 effects4 Gandolfi, Lapi, Liberati (e.g., Tollet et al. 2016; Freundlich et al. 2020b), that show a non-trivial mass dependence for the halo inner shape.Third, we can compute the halo concentration as c α ≡ R r α (using R r − yields similar result), which turns out to be c α (cid:39) (cid:18) M M (cid:12) (cid:19) − . E − . z , (21)in broad agreement, and actually slightly smaller than the outcome of N − body, DM-only simulations in the standardΛCDM cosmology (e.g., Bullock et al. 2001; Macci´o et al. 2007). SUMMARYWe have investigated self-gravitating equilibria of halos constituted by dark matter (DM) non-minimally coupled togravity. A non-minimal coupling may be present in modified gravity theories or it might be dynamically generated whenthe averaging/coherence length L associated to a fluid description of the DM collective behavior is comparable to thelocal curvature scale. We have theoretically motivated a form of such a coupling that in the Newtonian limit amountsto a modification of the Poisson equation by a term L ∇ ρ proportional to the Laplacian of the DM density ρ itself(see Sect. 2). We have further adopted an effective power-law equation of state p ∝ ρ Γ r α relating the DM dynamicalpressure p to density ρ and radius r , as expected by phase-space density stratification during the gravitational assemblyof halos in a cosmological context (see Sect. 3.1). In absence of the non-minimal coupling, we have confirmed previousfindings that the DM density run ρ ( r ) features a central density cusp and an overall shape mirroring the outcomesof N − body simulations in the standard ΛCDM cosmology, as described by the classic NFW or Einasto profiles (seeSect. 3.2).We have remarkably found that, when the non-minimal coupling is switched on, it causes the DM density profile todevelop an inner core and a shape closely following, out to several core scale radii, the Burkert profile (see Sect. 4). Inaddition, we have highlighted that our non-minimally coupled solutions can fit, with an accuracy comparable to theBurkert profile, the co-added rotation curve of DM-dominated dwarf galaxies (see Sect. 5). Finally, we have shownthat non-minimally coupled DM halos are consistent with the observed scaling relation between the core radius r andthe core density ρ in terms of an universal core surface density ρ × r among different galaxies, that has proven tobe challenging for many other DM models (see Sect. 6).A future development of this work will involve the study of the DM density profile in presence of baryons. Infact, the non-minimal coupling to gravity constitutes a natural and effective way to tightly link the DM and baryonproperties. On the one hand, this could help to understand puzzling scaling relationships between the DM and thebaryonic component, and to characterise the physical processes underlying the emergence of the non-minimal couplinglength-scale L . On the other hand, this will allow us to probe the effectiveness of our solutions in fitting the measuredrotation curves of normal, and not only dwarf, rotation-dominated galaxies. In parallel, we plan to extend the staticinvestigation pursued in the present paper to time-dependent conditions, by implementing the non-minimal couplinginside a full N − body numerical simulation. On more general grounds, it would be worth to explore the effect ofnon-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 crucialingredient toward an improved description of realistic DM structures in a cosmological framework. We very much hopethis novel perspective will contribute to shed light on some of the remaining mysteries concerning the DM componentin cosmic structures. ACKNOWLEDGMENTSWe 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 2017prot. 2017-3ML3WW, “Opening the ALMA window on the cosmic evolution of gas, stars and supermassive blackholes.” A.L. has taken advantage of the MIUR grant “Finanziamento annuale individuale attivit´a base di ricerca”and of the EU H2020-MSCA-ITN-2019 Project 860744 “BiD4BEST: Big Data applications for Black hole EvolutionSTudies.” S.L. acknowledges funding from the grant PRIN MIUR 2017 prot. 2017-MB8AEZ. quilibria of non-minimally coupled DM halos A. ANISOTROPIC CONDITIONSIn this Appendix we discuss the self-gravitating equilibria of non-minimally coupled DM halos when anisotropicconditions apply. These can be included in our treatment by modifying the second of Eqs. (7) as1 ρ d p d r + 2 β σ r r = − dΦd r , (A1)where β ≡ − σ θ σ r is the Binney (1978) anisotropy parameter in terms of the tangential and radial velocity dispersions σ θ and σ r , respectively. N − body simulations suggest β ( r ) to increase from central values β (cid:46)
0, meaning near isotropy,to outer values β (cid:38) .
5, meaning progressive prevalence of radial motions. This overall trend can be physicallyunderstood in terms of efficient dynamical relaxation processes toward the inner regions, that tend to enforce closelyisotropic conditions, while in the outskirts the infall energy of accreting matter is more easily converted by phasemixing into radial random motions (see Lapi et al. 2011 for details). Specifically, simulations suggest the effectivelinear expression (e.g., Hansen & Moore 2006) β ( r ) (cid:39) β + β [ γ ( r ) − γ ] = β − β γ − β r ρ (cid:48) ρ , (A2)in terms of the logarithmic density slope γ ( r ) ≡ − d log ρ d log r , with γ ≡ γ (0) being a value yet to be determined, β (cid:46) β ≈ . (cid:34) − η κ ¯ ρ − Γ (Γ − β ) ¯ r α (cid:35) ¯ ρ (cid:48)(cid:48) + (Γ −
2) ¯ ρ (cid:48) ¯ ρ ++ α (2Γ −
1) + 2 Γ + 2(Γ − β − α + 2 + (Γ − γ ] β Γ − β ¯ ρ (cid:48) ¯ r ++ [ α + 1] [ α + 2( β − β γ )]Γ − β ¯ ρ ¯ r − η κ ¯ ρ − Γ ¯ ρ (cid:48) (Γ − β ) ¯ r α +1 + κ ¯ ρ − Γ (Γ − β ) ¯ r α = 0 , (A3)Looking for powerlaw behaviors ¯ ρ (cid:39) ¯ r − γ one obtainsΓ (Γ − (cid:20) γ − α + 2 β Γ − β Γ ( γ − γ ) (cid:21) (cid:20) γ − α + 1Γ − (cid:21) − η κ γ ( γ − r γ (2 − Γ)+ α = − κ ¯ r γ (2 − Γ)+ α − , (A4)which, remarkably, allows to self-consistently determine the central slope γ .In fact, for minimally coupled halos ( η = 0), the anisotropic solutions feature a modified inner slope γ = α +2 β Γ with respect to the isotropic case, while retaining the same slopes at intermediate radii γ α = − α − Γ and in the outerregion γ ∞ = α +1Γ − ; as mentioned above, β ≈ β (cid:46) α crit = Γ (5Γ − β (Γ −
1) (Γ − − ;in particular, α crit = − β holds for for Γ = . The corresponding inner, intermediate and outer slopes read γ , crit = − β Γ3Γ − , γ α, crit = − β (Γ − − , and γ ∞ , crit = β (Γ − − , respectively.For non-minimally coupled halos, an inner core with γ ≈ η = − .
05, Γ = and α = − β = α crit is affected byanisotropies. For realistic values − . (cid:46) β (cid:46) . β ≈ .
2, the profile is marginally affected in the inner region6
Gandolfi, Lapi, Liberati
Figure A1.
Effects of realistic anisotropic conditions on the minimally-coupled density profile with η = − .
05; EOS parametersΓ = and α = − β = α crit have been adopted. Anisotropy profiles (see Appendix A for details) are described by theexpression β ( r ) = β + β ( γ − γ ), and are illustrated in the inset. Purple line refers to the reference profile in isotropicconditions with β = β = 0, red line to β = 0 and β = 0 .
2, blue line to β = − . β = 0 .
2, and orange line to β = +0 . β = 0 . and at intermediate radii, while it tends to extend toward slightly larger radii (i.e., the cutoff moves outward) due tothe progressive prevalence of radial anisotropy in the halo outskirts; such an effect, though minor, is more pronouncedfor larger (more positive) β . REFERENCES
Adhikari, R., Agostini, M., Ky, N.A., et al. 2017, JCAP, 01, 025An, J.,& Zhao, H. 2013, MNRAS, 428, 2805Arora, A., & Williams, L.R. 2020, ApJ, 893, 53Austin, C. G., Williams, L. L. R., Barnes, E. I., Babul, A., &Dalcanton, J. J. 2005, ApJ, 634, 756Bekenstein, J.D. 2004, PhRvD, 70, 083509 [erratum: 2005,PhRvD, 71, 069901]Bekenstein, J.D. 2009, NuPhA, 827, 555Bekenstein, J.D. 1993, PhRvD, 48, 3641Beltran Jimenez, J., Heisenberg, L., Olmo. L.G. J. & Rubiera-Garcia, D. 2018, Phys. Rept. 727, 1Bernal, T., Fernandez-Hernandez, L.M., Matos, T., & Rodriguez-Meza, M.A. 2018, MNRAS, 475, 1447Bertolami, O., & Paramos, J. 2010, JCAP, 3, 9Bertolami, O., & Paramos, J. 2014, PhRvD, 89d4012Bertschinger, E. 1985, ApJS, 58, 39Bettoni, D., Liberati, S., & Sindoni, L. 2011, JCAP, 11, 007Bettoni, D., Pettorino, V., Liberati, S., & Baccigalupi, C. 2012, 7,27Bettoni, D., & Liberati, S. 2013, PhRvD, 88, 4020Bettoni,D., Colombo, M., & Liberati, S. 2014, JCAP, 004, 1402Bettoni, D., & Liberati, S. 2015, JCAP, 008, 023Binney, J. 1978, MNRAS, 183, 779Bohmer, C. G., & Harko, T. 2007, JCAP, 6, 25 Bruneton, J.-P., Liberati, S., Sindoni, L., & Famaey, B. 2009,JCAP, 3, 21Bruneton, J.-P., & Esposito-Far`ese, G. 2007, PhRvD, 76l4012Bullock, J.S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343Bullock, J.S., Kolatt, T.S., Sigad, Y., et al. 2001, MNRAS, 321,559Burkert, A. 2020, ApJ, 904, 161Burkert, A. 2015, ApJ, 808, 158Burkert, A. 1995, ApJ, 447, L25Butsky, I., Maccio, A.V., Dutton, A.A., et al. 2016, MNRAS, 462,663Chae, K.-H. 2014, ApJ, 788, L15Chan, M.H. 2019, Sci. Rep., 9, 3570de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136,2648Dehnen, W., & McLaughlin, D. E. 2005, MNRAS, 363, 1057Deng, H., Hertzberg, M.P., Namjoo, M.H., & Masoumi, A. 2018,PhRvD, 98, 023513Destri, C., de Vega, H. J., & Sanchez, N. G. 2013, NewA, 22, 39Di Casola, E., Liberati, S. & Sonego, S. 2015, Am. J. Phys, 83, 39Dodelson, S., & Widrow, L.M. 1994, PhRvL, 72, 17Donato, F., Gentile, G., Salucci, P., et al. 2009, MNRAS, 397,1169Donato, F., Gentile, G., Salucci, P., et al. 2004, MNRAS, 353, L17 quilibria of non-minimally coupled DM halos17