Empirical Evidence of Non-Minimally Coupled Dark Matter in the Dynamics of Local Spiral Galaxies?

We look for empirical evidence of a non-minimal coupling (NMC) between dark matter (DM) and gravity in the dynamics of local spiral galaxies. In particular, we consider a theoretically motivated NMC that may arise dynamically from the collective behavior of the coarse-grained DM field (e.g., via Bose-Einstein condensation) with averaging/coherence length $L$. In the Newtonian limit, this NMC amounts to modify the Poisson equation by a term $L^{2} \nabla^{2} \rho$ proportional to the Laplacian of the DM density itself. We show that such a term, when acting as a perturbation over the standard Navarro-Frenk-White (NFW) profile of cold DM particles, can substantially alter the dynamical properties of galaxies, in terms of their total radial acceleration within the disk and rotation velocity. Specifically, we find that this NMC model can properly fit the stacked rotation curves of local spiral galaxies with different velocities at the optical radius, including dwarfs and low-surface brightness systems, at a level of precision comparable to, and in some instances even better than, the phenomenological Burkert profile. Finally we show that, extrapolating down to smaller masses the scaling of $L$ vs. halo mass found from the above rotation curve analysis, the NMC model can adequately reproduce the radial acceleration relation (or RAR) in shape and normalization down to the dwarf spheroidal galaxy range, a task which constitutes a serious challenge for alternative DM models even inclusive of baryonic effects.


INTRODUCTION
The analysis of spiral galaxy rotation curves (RCs) has empirically highlighted since the late 1970s a discrepancy between the amount of luminous matter and the mass budget required to explain the overall kinematic properties of such systems (see Rubin et al. 1978;Bosma 1978). The common lore traces back the missing mass to an unseen component called dark matter (DM), constituted by cold (i.e. non-relativistic) and weakly interacting massive particles. Despite such a cold DM paradigm has proven to be relatively successful on cosmological scales, it struggles to fully describe the observed phenomenology on galactic scales, especially in DM-dominated dwarfs. In this respect, two crucial issues will be our focus here.
The first is the well-known cusp-core controversy about the inner shape of the DM density profile in galaxies. Analysis of observed galaxy RCs in the standard Newtonian framework seems to favor an inner core, while gravityonly simulations in the standard cold DM framework produce the universal Navarro-Frenk-White (NFW; see Boylan-Kolchin & Ma 2004;Navarro 2006;de Blok 2010;Navarro et al. 2017) profile with a cuspy inner behavior. The second, somewhat related, point concerns the so called radial acceleration relation (RAR), linking the (total) radial acceleration g tot inferred from galaxy RCs with different masses/velocities and that associated to the luminous matter distribution g bar mainly probed by photometric observations. The RAR is a remarkably tight relationship (see e.g. Lelli et al. 2017b;, that subsumes/generalizes many well-known dynamical laws of galaxies (see Lelli et al. 2017b), and has been extensively studied (e.g. Burrage et al. 2017;Lelli et al. 2017a;Chae et al. 2019;Di Paolo et al. 2019;Green & Moffat 2019;Tian et al. 2020;Rodrigues & Marra 2020a).
From the theoretical point of view, a plethora of physical effects have been invoked to clear the above issues. As for the cusp-core problem, it has been advocated that dynamical friction (El-Zant et al. 2001;Tonini et al. 2006;Romano-Diaz et al. 2008;Goerdt et al. 2010;El-Zant et al. 2016) or 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;Freundlich et al. 2020b) during the galaxy formation process can induce violent fluctuations in the inner gravitational potential and/or transfer of energy and angular momentum from the baryons to DM, so erasing the central cusp. Another class of solutions conceives DM as constituted by non-standard particle candidates, thus leading to abandon the cold DM hypothesis in favor of more exotic alternatives (see the review by Salucci 2019 and references therein). As for the RAR, it has been claimed to emerge naturally from the self-similarity of cold DM halos (e.g., Navarro et al. 2017) or to be explained by properly accounting for the effects of baryons (e.g., Di Cintio et al. 2014;Di Cintio & Lelli 2016;Santos-Santos et al. 2016;Ludlow et al. 2017;Desmond 2017;Navarro et al. 2017;Wheeler et al. 2019). Note that it has been hinted that cored profiles could have a better chance of correctly reproducing the RAR (Di Cintio & Lelli 2016, Di Paolo et al. 2019, Li et al. 2020.
On a different ground, the incompleteness of the cold DM model on galactic scales has led to the emergence of numerous theories of modified gravity. Perhaps the most famous framework on galactic scales is the Modified Newtonian Dynamics (MOND), that was originally proposed by Milgrom (1983a) and further investigated in a rich literature (see Bekenstein 2004;Bruneton & Esposito-Farèse 2007;Milgrom 2010;Famaey & McGaugh 2012). As the name suggests, MOND aims to explain the mass discrepancy in galaxies through a modification of Newtonian gravity (or more generally of the Newton second law) that comes into action at accelerations well below a definite universal threshold; in its original formulation, DM is not included and baryons are the only source of the gravitational field. As for the two aforementioned issues, it has been claimed that MOND (or theories reducing to it in the weak field limit) can properly fit galactic RCs (de Blok & McGaugh 1998;Sanders & McGaugh 2002), and provide a satisfying description of the RAR (e.g., . In this paper, we entail yet another viewpoint to modify the standard cold DM framework and make it capable of accurately describe the observed galaxy RCs and at the same time to faithfully reproduce the RAR. Specifically, we put forward the possibility that DM could be non-minimally coupled to gravity, as conjectured in a series of previous works from our team and collaborators (see Bruneton et al. 2009;Bertolami & Paramos 2010;Bettoni et al. 2011;Bettoni et al. 2014;Bettoni & Liberati 2015;Ivanov & Liberati 2020;Gandolfi et al. 2021). Introducing such a coupling can retain the success of the cold DM on large cosmological scales while improving its behavior in galactic systems, recovering there a MOND-like (even if not exactly MONDian, since DM is there) dynamics. The word "nonminimal" in this context means that DM, or more precisely its gradients, are directly coupled to the Einstein tensor. We caveat that such non-minimal coupling (NMC) is not necessarily a fundamental feature of the DM particles, but rather may dynamically develop when the averaging/coherence length L associated with the fluid description of the DM collective behavior is comparable to the local curvature scale. In the Newtonian limit the NMC here considered implies a modification of the Poisson equation by a term L 2 ∇ 2 ρ proportional to the DM density ρ (as in Bettoni et al. 2014). This apparently simple addition can significantly change the internal dynamics of galaxies with respect to a pure cold DM framework, and in fact it has already proven to alleviate some problems in DM-dominated systems (see Gandolfi et al. 2021). Incidentally, note that on large scales non-minimally coupled fluids behave under certain conditions similarly to a repulsive dark energy component, and thus the NMC could have a cosmological relevance too (e.g., Bettoni et al. 2011;Bertolami & Páramos 2013).
We will show that the NMC term, when acting as a perturbation on a galaxy system characterised by the cuspy NFW profile for the DM, can substantially alter its dynamical properties. Such a NMC model can thus provide accurate fits to the stacked RCs of spiral galaxies with different velocities at the optical radius, including dwarfs and low-surface brightness systems. Moreover, we will show that the same NMC model can properly account for the RAR as well. The plan of the paper is as follows: in Sect. 2 we briefly review the theoretical background behind our NMC model; in Sect. 3 we analyse a sample of stacked RCs of spiral galaxies with the NMC model; in Sect. 4 we build empirically-based mock RCs of galaxies with different masses and construct the RAR, showing that it is well reproduced by the NMC model; in Sect. 5 we summarize our findings and outline future perspectives and applications of the NMC framework.

A THEORETICAL BACKGROUND FOR THE NMC
In this section, we recall the basic theoretical background behind the NMC hypothesis, referring the reader to Gandolfi et al. (2021) for further details. A very basic NMC model can be built by adding a coupling term S int between DM and gravity in the total Einstein-Hilbert action (in the Jordan frame) with shape: (1) here ϕ is the (real) DM scalar field, = ±1 is the polarity of the coupling, G µν is the Einstein tensor, and L is the NMC characteristic length-scale. Note that L may be not a new fundamental length-scale of Nature, but rather can emerge dynamically from some collective behavior of the coarse-grained DM field (e.g., Bose-Einstein condensation). Therefore, such a NMC model does not consist in a modified gravity theory, but simply in a formalization of an emergent behavior of cold DM in galactic environments. Note that, from a purely theoretical perspective, such form of the NMC is allowed by the Einstein equivalence principle (e.g., Di Casola et al. 2015). We will keep indicated as a bookkeeping parameter, but based on Gandolfi et al. (2021) we will set it to = −1 (repulsive coupling). Adopting the fluid approximation for the field ϕ (as in Bettoni et al. 2012) and taking the Newtonian limit, it can be shown that the NMC boils down to a simple modification of the usual Poisson equation (Bettoni et al. 2014, Gandolfi et al. 2021) where Φ is the Newtonian potential, and ρ bar and ρ are the baryonic and DM densities. In spherical symmetry, Eq.
(2) implies that the total gravitational acceleration writes where M (< r) is the total mass enclosed in the radius r; the first term is the usual Newtonian acceleration and the second term is the additional contribution from the NMC. Plainly, the related RCs v 2 tot (r) = |g tot (r)| r of spiral galaxies predicted in this framework will differ from the standard Newtonian case.
In Gandolfi et al. (2021) we have highlighted that Eq.
(2) gives rise to some interesting features for strongly DMdominated systems in self-gravitating equilibria: the NMC can help to develop an inner core in the DM density profile, enforcing a shape closely following the Burkert one out to several core scale radii; DM-dominated halos with NMC are consistent with the core-column density relation (see e.g. Salucci & Burkert 2000, Donato et al. 2009, Burkert 2015, Behroozi et al. 2013, Burkert 2020, i.e. with the observed universality of the product between the core radius r 0 and the core density ρ 0 . However, the NMC hypothesis needs still to be tested in galaxies with different velocities at the optical radius, where the contribution of the baryonic component to the dynamics can be substantial, which is precisely our aim in the next sections.

TESTING THE NMC WITH STACKED RCs
In this section, we will apply the NMC to mass-model stacked RCs of local spiral galaxies with different velocities at the optical radius and related properties. We will then compare our results with fits obtained from the standard Newtonian case for two other classic DM halo shapes, namely the standard NFW profile (see Navarro et al. 1996) emerging from gravity-only simulations of cold DM, and the phenomenological Burkert profile (see Burkert 1995).
For the analysis, we rely on the samples of stacked RCs collected by Persic et al. (1996) for normal spirals divided in 11 average velocity bins, by Dehghani et al. (2020) for low surface brightness (LSB) spirals divided in 5 average velocity bins and by Karukes & Salucci (2017) for low-luminosity dwarfs. These stacked RCs are built by co-adding high-quality individual RCs of thousands galaxies with similar velocities at the optical radius and related properties, after properly normalizing the velocity and radial variables to reference scales for each galaxy, which are typically the optical radius r opt and the optical circular velocity v opt ≡ v(r opt ); the interested reader can find details of such a procedure in Lapi et al. (2018). The average properties of our sample of stacked RCs are listed in Table 1.
We mass-model the stacked RCs as the sum of a baryonic (disk) component v 2 d (r) = G M d (< r)/r plus a DM contribution v 2 DM (r) = G M DM (< r)/r, with M d (< r) and M DM (< r) the cumulative disk and DM mass, respectively. The overall velocity model plainly reads v 2 tot (r) = v 2 d (r) + v 2 DM (r). The distribution followed by baryonic matter is modeled as a razor-thin exponential disk (see Freeman 1970) with exponential surface density Σ d (r) = Σ 0 exp (−r/r d ) ; here Σ 0 = M d /2π r 2 d is the central value in terms of the total disk mass M d = M d (< ∞) and of the disk scale-length r d ≈ r opt /3.2. The related contribution to the RC is given by (e.g., Binney & Tremaine 1987) where y ≡ r/(2 r d ), while I 0,1 (·) and K 0,1 (·) are modified Bessel functions. Since the fit is performed in a radial range r r opt we have checked that any contribution from a gaseous disk (typically more important at larger radii) is negligible and largely unconstrained, so we include only the stellar disk in the mass-modeling.

DM models
We exploit three different DM models. Two are based on standard Newtonian gravity, but differ in the form of the DM profile shape: NFW or Burkert. The other model is based on the NFW profile but include a perturbative correction to the dynamics via the NMC term of Sect. 2.
• NFW profile The standard NFW profile features the shape (see Navarro et al. 1996;Lokas & Mamon 2001) where δ c is the (dimensionless) characteristic overdensity of the halo, ρ c = 3 H 2 0 /8π G is the local critical density, and r s is the scale radius. The virial mass M v and the concentration c ≡ r v /r s , defined in terms of the virial radius r v ≈ 260 (M v /10 12 M ) 1/3 , can be used to fully characterize the profile since where s ≡ r/r v . From the above, it is clear that the overall galaxy RC can be specified in terms of three parameters: the halo mass M v , the halo concentration c and the disk mass M d .

• NMC model
We include the effect of the NMC as a perturbative correction to the dynamics based on Eq.
(3), retaining the standard NFW profile for the DM. The perturbative parameter in our analysis is the term L 2 /r 2 s , which, as we will show with our results, is always small for the range of masses probed in our study. Plugging Eq. (5) in Eq.
(3) and after some simple yet tedious algebra we obtain the RC v 2 DM (r) = The overall RC model has four free parameters: the halo concentration c, the halo virial mass M v , the NMC length-scale L, and the disk mass M d .

• Burkert profile
The phenomenological Burkert profile features the shape where r 0 is the core radius and ρ 0 the core density. The RC writes (see Salucci & Burkert 2000) v 2 DM (r) = 4 GM 0 r ln 1 + r r 0 − tan −1 r r 0 where M 0 = 1.6 ρ 0 r 3 0 . When using the Burkert profile, we adhere to the customary approach of describing the total RC in terms of three parameters: the core radius r 0 , the core mass M 0 , and the ratio κ ≡ v 2 d (r opt )/v 2 tot (r opt ) of the disk to the total velocity at the optical radius.

Fitting procedure and results
We fit the stacked RC data with the mass models described above, exploiting the emcee python package for Bayesian Monte Carlo Markov Chain parameter estimation (see Foreman-Mackey et al. 2013). We present here representative outcomes concerning one velocity bin for each of the galaxy type: normal spirals, LSBs and dwarfs; the complete analysis for all the other velocity bins produces similar results and is reported in the Appendix and in Tables 2-3-4. First, we consider Bin 5 from the Persic et al. (1996) sample of spiral galaxies, whose outcome is reported in Fig. 1.
The results on the estimated virial mass are consistent for the three profiles. As to the disk mass, it is consistent between NMC and Burkert models, while for the NFW model only a rather loose upper limit can be derived. All in all, the NMC model curve performs appreciably better in terms of reduced χ 2 red ≈ 0.6 with respect to the Burkert χ 2 red ≈ 22.5 and to the pure NFW model χ 2 red ≈ 11, as can be also appreciated graphically. The estimated value of the NMC length-scale is around 0.2 kpc, roughly corresponding to a sixtieth of r s . We also try to perform the fits of the NFW and NMC models by imposing the concentration parameter of the halo to satisfy the relation with the virial mass by Dutton & Macciò (2014). We find that both fits are not appreciably altered, but the posterior distribution of the fitted parameters in the NMC model are still consistent and somewhat narrowed. Fig. 2 refers to Bin 5 in the sample of LSB galaxies by Dehghani et al. (2020). In this case, the Burkert model yields a reduced χ 2 red ≈ 11, the NFW fit yields χ 2 red ≈ 3 and the NMC model performs better yielding χ 2 red ≈ 1.411. The disk mass in all the fits is poorly constrained, as expected since these LSB galaxies have an extremely extended disk mass distribution relative to the region probed by the RC.
Finally, in Fig. 3 the dwarf galaxy bin is analysed. Since it was originally designed on purpose, it is not surprising that in this case the Burkert profile yields the best description of the RC with a reduced χ 2 red ≈ 0.8. However, the NMC model performs decently with χ 2 red ≈ 4, and substantially better than the NFW profile for which χ 2 red ≈ 14. Note that such galaxies are strongly DM dominated in the region probed by the RC, hence the disk mass is vanishingly small and/or unconstrained by all models.
An overall interesting result is that the NMC model predicts higher values of the length-scale L in DM halos of higher virial masses -see Fig. (4). This trending is well reproduced by the scaling L(M v ) ∝ M 0.8 v , a result consistent with the findings of Gandolfi et al. (2021) for DM-dominated dwarf galaxies.
As can be seen by looking at the overall results listed in the Appendix and recapped in Tables 2-3-4, the NMC model yields RC fits that are always superior to the pure NFW one and in several instances comparable or even better than the Burkert model. Furthermore, we performed an F-test to compare the NFW and the NMC models, see Table  4. Overall, the test suggests that the addition of the parameter L effectively improves the fits for the majority of the bins. Two caveats are in order here. First, the Burkert profile is phenomenological, and has been designed specifically to fit the RC of dwarf galaxies. Contrariwise, our NMC model is derived theoretically from first principles (though in a specific scenario), so the fact that its performances on RC fitting for different kind of galaxies improves substantially over the pure NFW shape is in itself encouraging. Second, we will show in the next Section that the NMC model will perform better than the Burkert profile in reproducing the RAR.

TESTING THE NMC WITH THE RAR
The radial acceleration relation (or RAR) was originally proposed in McGaugh et al. (2016) by exploiting the individual high-quality RCs of the SPARC sample (see Lelli et al. 2016a). As argued in Lelli et al. (2017b), the RAR subsumes and generalizes a plethora of well-known dynamical laws of galaxies, such as the Baryonic Tully-Fisher relation (McGaugh 1999;Wheeler et al. 2019), the dichotomy between high and low surface brightness galaxies (de Blok & McGaugh 1997;Tully & Verheijen 1997), and others (Lelli et al. 2013;Lelli et al. 2016b;van Albada & Sancisi 1986;Sancisi 2004;Faber & Jackson 1976;Serra et al. 2016).
In Lelli et al. (2017b), an overall representation of the RAR was introduced in terms of the function with g † andĝ being fitting parameters. Eq. (10) withĝ = 0 accurately represents the RAR for spiral and irregulars, while the additive term depending onĝ describes the flattening of the RAR in the typical acceleration regime proper of dwarf spheroidal (dSphs) galaxies. All in all, the values g † = (1.1±0.1)×10 −10 m s −2 andĝ = (9.2±0.2)×10 −12 m s −2 are derived from the analysis of the overall SPARC sample. In the recent literature, it was argued that the parameter g † could represent an acceleration scale governing the average internal dynamics of galaxies. Since the existence of such a scale in the standard cosmological model is far from trivial, this phenomenon has been interpreted as a possible sign of modified gravity (e.g., Notice, however, that such an interpretation is still highly debated, with some works supporting it (see, e.g., Ghari et al. 2019) and some others ruling it out (e.g., Marra et al. 2020;Rodrigues & Marra 2020b). Our aim here is to determine whether our NMC model can adequately reproduce the RAR, and whether it can do so with values on the NMC length-scale that are consistent with those derived from the previous analysis of stacked RC data. Toward this purpose, first notice that the RAR is a local scaling law that combines data at different radii in galaxies with different masses, which feature different contribution of stellar disc and bulge, gas and DM. Therefore we approach the problem via a semi-empirical method: we first build up mock RCs of galaxies with different properties, and then we sample them to derive the total and baryonic accelerations and construct the RAR, as detailed below.

Mock RC modelling
Our procedure to build up mock RCs consists of the following steps.

• DM mass
We start by randomly drawing a very large number of total DM halo masses M v within the range 8 < log(M v /M ) < 13.3 according to the local halo mass function (a uniform sampling does not impact appreciably the final outcomes).

• Stellar mass
We then derive the stellar mass associated to each galaxy by using the relation found by Behroozi et al. (2013) through abundance matching technique: with log M 1 = 11.514 being a characteristic halo mass, and parameters log ε = −1.777, α = −1.412, δ = 3.508, γ = 0.316. We allow for a log-normal scatter of 0.25 dex.
• Gas mass We determine the gas mass by exploiting the relation found with the stellar mass by Papastergis et al. (2012) and Peeples et al. (2014): allowing a lognormal scatter of 0.15 dex. Note that we are implicitly assuming that in local galaxies the majority of the intestellar medium consists in atomic hydrogen HI and that both the ionized and the molecular components are minor (see Papastergis et al. 2012, Saintonge et al. 2011). The total gas mass is M gas ≈ 1.33 M HI to account for the contribution of He.
• Stellar and gas radial distributions We assume that the gaseous and the stellar components are both distributed in a razor-thin exponential disk. We determine the stellar disk half-mass radius from the stellar mass via the relation by Shen et al. (2003) log applying for M < 10 9 M , and that by Lange et al. (2015) R e = 0.13 M M 0.14 1.0 + M 14.03 · 10 10 M 0.77 kpc (14) for M ≥ 10 9 M . The stellar disk scale-length is given by R d ≈ R e /1.678. We allow for a lognormal scatter of both these relations around 0.1 dex. The gas distribution scale-length is taken as R gas = 2 R d .
• Bulge mass We determine the bulge mass using the relation with the stellar mass by Gadotti (2009) and Moran et al. (2012): with a lognormal scatter of 0.1 dex. The implied bulge-to-total mass ratio is ∼ 30% for Milky-Way like galaxies. We further assume that the bulge is present only if M ≥ 3 × 10 9 M .
• Bulge radial distribution We assume that the bulge mass is radially distributed according to an Hernquist profile (see Hernquist 1990) where R 1/4 is the radius at which the enclosed bulge mass is a quarter of its total value. The half-mass radius R 1/2 = (1 + √ 2)R 1/4 is gauged on the basis of the scaling relation with the bulge mass by Gadotti (2009) log with a lognormal scatter of 0.1 dex.
• DM radial distribution We radially distribute the DM mass according to various profiles, to test their performance on the RAR. For the NFW and the NMC models we use Eq. (5); the concentration parameter c is determined according to the relation with the halo mass from Dutton & Macciò (2014) log with a lognormal scatter of 0.11 dex.
For the Burkert model we use Eq. (8) by setting the core radius r 0 from two conditions: (i) the mass within the virial radius must match M v ; the core radius and core density must satisfy the universal core-column density relation ρ 0 × r 0 ≈ 75 M pc −2 , with a scatter of 0.2 dex (see e.g. Salucci & Burkert 2000, Donato et al. 2009, Burkert 2015, Behroozi et al. 2013, Burkert 2020).
Finally, we consider the profile emerging from the hydrodynamical simulations by Di Cintio et al. (2014) which take into account DM responses to baryonic effects (including stellar feedback); this is basically a generalized NFW profile with shape parameters linked to the stellar-to-halo mass ratio X ≡ M /M v (see also Stinson et al. 2013 • Building up the mock RC For any galaxy of given virial mass M v , we have now specified all the mass components and the associated radial distribution M i (< R), so that the RC can be easily determined from v 2 i (R) = G M i (< R)/R; the only exception is the NMC model for which the DM velocity has an additional term v 2 DM (r) = G M DM (< r)/r − L 2 r 4π G dρ/dr. The overall mock RC is then the sum of all the different contributions v 2 tot = i v 2 i . In Fig. 5 we illustrate four representative mock RCs for galaxies with different halo masses M v , highlighting the diverse behavior when assuming the NFW, the Burkert, the Di Cintio or the NMC halo profiles. As for the baryonic components, in moving toward smaller halo masses, the inner contribution due to the bulge component becomes less prominent, while the gas contribution progressively increase to become even dominant over the stellar disk. As for the DM models, the halo shapes are rather different, with the Burkert profile yielding overall higher velocities in lower mass galaxies. In order to further test the realism of our mock RCs, we compared them to the stacked empirical RCs considered in Sec. (3). The outcome is shown in Fig. (6) showing the compatibility between the mock curves and empirical, stacked ones.

Building the RAR and results
Once the mock RC for each mock galaxy has been characterized, we compute the gravitational acceleration following with the index j = bar, tot specifying the baryonic contribution or the total value including DM. The RAR is then constructed by binning our mock galaxy sample in g bar and extracting the average values and standard deviation of g tot . For fair comparison with the data, the sampled portion of the RC is restricted to twice the optical radius of each mock galaxy. It is clear that this procedure includes in a given bin of g bar objects with different halo masses and at different radii; e.g., an object can display a low baryonic acceleration either because it has a small halo mass or because its RC is sampled at large radii.
In Fig. 7 we illustrate our results on the RAR, for the DM models listed above (color-coded). For comparison, we report as a black line with shaded area the determination by Lelli et al. (2017b) represented by Eq. (10); binned data for spirals and irregulars are represented by grey squares and individual data for dSph galaxies are highlighted with diamonds (filled symbols are for more secure determinations). There is a substantial agreement of the RAR for all the DM models at high baryonic accelerations. This is because such a regime is mainly dominated by the contribution at small radii in high mass galaxies. There the total gravitational acceleration is anyway dominated by baryons, implying g tot ≈ g bar irrespective of the specific DM profile. However, a marked difference among the RAR for different DM models sets in toward lower baryonic accelerations. Such a regime is dominated by the behavior at small/intermediate radii in intermediate and low mass galaxies. There the total baryon acceleration g bar is dominated by the stellar disk, while the total gravitational acceleration is contributed by both the disk and the halo g tot ≈ g DM + g bar ; thus depending on the DM model most of the contribution to g tot may come either from the disk enforcing a behavior of the RAR similar to the high acceleration regime, or from the DM enforcing an upward deviation of the RAR.
All in all, both the RAR associated to the Di Cintio and the Burkert models tend to appreciably deviate downward, to the point of becoming inconsistent with the measured RAR (especially in dSph) for g bar 10 −11 m s −2 . Contrariwise, the RAR of the NFW model displays the opposite behavior, with the corresponding curve flattening and progressively saturating to values slightly above the observed RAR, though still consistent with the upper outliers; nevertheless, one must keep in mind that NFW model suffers from the poor performances in fitting the individual RCs of many dSphs (e.g., de Blok 2010), and the stacked dwarf galaxy RCs analysed in this paper. Finally, the NMC model can reproduce the average measured RAR when extrapolating down to smaller masses, the dependence L(M v ) ∝ M 0.8 v found from the RC analysis of Sec. (3) and to be consistent in the dwarf irregular regime with Eq. (9) of Gandolfi et al. (2021). We find that the RAR thus produced follows a profile intermediate between the NFW and the cored ones. Remarkably, the NMC one is the only model considered here that can simultaneously reproduce the RAR and decently fits the stacked RC of spirals, LSBs and dwarf galaxies.
For reference, we also illustrate the prediction on the RAR for the MOND framework. This may be derived from the relation µ(x) g tot = g bar where the simple interpolating function µ(x) = x/(1 + x) with x ≡ g tot /a 0 and a 0 ∼ 1.2 × 10 −10 m s −2 is generally adopted (e.g., Famaey & Binney 2005, Zhao & Famaey 2006. The MOND outcome is quite close to the measured RAR at high acceleration, while it lacks the progressive flattening at low g bar . However, some authors pointed out that this simple parameterization of MOND is not accurate because of the so called external field effect (EFE; see Milgrom 1983b;Swaters et al. 2010;Candlish et al. 2018) associated to the violation of the strong equivalence principle in the relativistic MOND theory and implying that a galaxy dynamics depends also on the gravitational pull g ext from external fields (e.g., external galaxies or large-scale surroundings). One can account for the EFE by modifying the interpolating function to read µ(x) = (x/1 + x + e) [1 + (2 + e) e/x(1 + e)], with e = g ext /a 0 being the strength of the effect with respect to the MOND acceleration scale (see Timberlake et al. 2021). This parameter was estimated to be around e ≈ 0.033 by Chae et al. (2020) and Chae et al. (2021) from the analysis of individual galaxy RCs (Desmond et al. 2018). The RAR from MOND theory including the EFE deviate downward at low accelerations, and can possibly account for some of the bottom outliers. However, to reproduce the observed RAR for the bulk of the galaxies would require to have negative values of e, which are not supported by observational estimates and known to be theoretically unfeasible in the MOND framework (see Chae et al. 2020).
In Fig. 8 we plot the RAR expected from the NMC model for different values of the NMC length-scale L/r s . Plainly, for vanishing L/r s the NFW outcome is recovered. For L/r s progressively increasing, the NMC model spans the dispersion of the outliers in the RAR at low baryonic accelerations.

SUMMARY
We have looked for empirical evidence of a non-minimal coupling (NMC) between dark matter (DM) and gravity in the dynamics of local spiral galaxies. In particular, taking up the work by Gandolfi et al. (2021) we have considered a theoretically motivated non-minimal coupling that may arise dynamically from some collective behavior of the coarsegrained DM field (e.g., Bose-Einstein condensation) with averaging/coherence length L. In the Newtonian limit, this non-minimal coupling amounts to modify the Poisson equation by a term L 2 ∇ 2 ρ proportional to the Laplacian of the DM density itself.
We have then worked out how such a term, when acting as a perturbation over the standard NFW profile from gravity-only cold DM simulations, can substantially alter the dynamical properties of galaxies, in terms of their total radial acceleration within the disk and rotation velocity. We have then tested such a model against dynamical data of local spiral galaxies. Our main results are as follows: • We have found that the NMC model can fit the stacked RCs of local spiral galaxies with different average velocities at the optical radius, including dwarfs and low-surface brightness systems, at a level of precision superior to the NFW profile and comparable to (in some instances even better than) the phenomenological Burkert profile.
• We have shown that at the same time our NMC model, when extrapolating down to smaller masses the massdependent scaling of the coupling length-scale L found from the RC analysis, can adequately reproduce in shape and scatter the radial acceleration relation (RAR) down to the dwarf spheroidal galaxy range, a task which constitutes a serious challenge for alternative DM profiles even inclusive of baryonic effects.
A possible future extension of the present work may include to trace the physical origin of the NMC in DM halos, especially in terms of the mechanism determining the NMC length-scale L in different galaxies and originating the mass dependence required to fit the RAR at very low baryonic accelerations. In this vein, full N −body simulations incorporating the NMC hypothesis could be exploited to study time-dependent conditions and the overall formation of cosmic structures in such a framework.
We warmly thank T. Ronconi and P. Salucci for the helpful insights and stimulating discussions. AL is supported by the EU H2020-MSCA-ITN-2019 Project 860744 BiD4BESt: Big Data applications for black hole Evolution STudies, and by the PRIN MIUR 2017 prot. 20173ML3WW: Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes.   . The Radial Acceleration Relation or RAR. The solid black line with shades illustrate the average results and its 2σ and 3σ variance from the analysis of the SPARC database by Lelli et al. (2017b); in particular, grey squares refer to the binned outcome for normal spiral galaxies and diamonds to measurements in individual dwarf spheroidal (filled symbols are more secure determinations). Such dwarf spheroidals have large error bars that are not displayed in this plot for visual clarity, thus the extension of the fit line through this cloud is much less certain than for the LTGs. The colored circles illustrate the prediction from our empirical modeling of RCs when adopting different halo profiles: NFW (green), Di Cintio (blue), Burkert (cyan) and NMC with a mass dependent scaling for the coupling length-scale L (red; see text for details). For reference, the MOND expectations without (dashed orange) and with (dotted orange) the external field effect (to the value e = +0.033 estimated in Chae et al. 2020) is displayed. . The RAR for the NMC model, with different coupling length-scales L/rs. as expected, setting the coupling length to zero (orange) amounts in recovering the RAR reproduced by the NFW model. We also display the RAR for L/rs ∼ 0.077, i.e. the average value obtained from the RC analysis of large spiral galaxies (purple). Intermediate values for L/rs describe RARs that will lie between these two extremes. For reference the RAR obtained assuming a mass-dependent scaling for the coupling length-scale L as in previous Figure is also reported (red).  Persic et al. (1996), LSB stands for the sample of low surface brightness spirals from Dehghani et al. (2020), and Dw for the sample of dwarfs by Karukes & Salucci (2017). For each bin the optical radius ropt and optical velocities vopt are reported.