Testing the Strong Equivalence Principle: Detection of the External Field Effect in Rotationally Supported Galaxies

The strong equivalence principle (SEP) distinguishes general relativity (GR) from other viable theories of gravity. The SEP demands that the internal dynamics of a self-gravitating system under freefall in an external gravitational field should not depend on the external field strength. We test the SEP by investigating the external field effect (EFE) in Milgromian dynamics (MOND), proposed as an alternative to dark matter in interpreting galactic kinematics. We report a detection of this EFE using galaxies from the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample together with estimates of the large-scale external gravitational field from an all-sky galaxy catalog. Our detection is threefold: (1) the EFE is individually detected at 8σ to 11σ in “golden” galaxies subjected to exceptionally strong external fields, while it is not detected in exceptionally isolated galaxies, (2) the EFE is statistically detected at more than 4σ from a blind test of 153 SPARC rotating galaxies, giving a mean value of the external field consistent with an independent estimate from the galaxies’ environments, and (3) we detect a systematic downward trend in the weak gravity part of the radial acceleration relation at the right acceleration predicted by the EFE of the MOND modified gravity. Tidal effects from neighboring galaxies in the Λ cold dark matter (CDM) context are not strong enough to explain these phenomena. They are not predicted by existing ΛCDM models of galaxy formation and evolution, adding a new small-scale challenge to the ΛCDM paradigm. Our results point to a breakdown of the SEP, supporting modified gravity theories beyond GR.


INTRODUCTION
The hypothesis that General Relativity (GR) and its Newtonian limit hold exactly in the weak gravity regime Corresponding author: Kyu-Hyun Chae KHC: chae@sejong.ac.kr, kyuhyunchae@gmail.comFL: LelliF@cardiff.ac.ukHD: harry.desmond@physics.ox.ac.ukSSM: ssm69@case.eduPL: pxl283@case.eduJMS: jschombe@gmail.comrequires that the Universe is permeated by invisible dark matter (DM).The existence of DM is a key assumption of the standard cosmological model Λ Cold Dark Matter (ΛCDM), which has been successful in explaining many cosmological observations on the largest scales of the cosmos (Peebles 2012;Frenk & White 2012).The ΛCDM paradigm, however, is facing several challenges on small scales (Bullock & Boylan-Kolchin 2017;Kroupa 2015), such as the unexpected phase-space correlation of satellite galaxies ("the satellite plane problem"; see, e.g., Kroupa et al. 2010;Müller et al. 2018) and the unexpected coupling in galaxies between the visible matter (baryons) and the observed dynamics, usually dominated by the DM halo at large radii (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017).

Chae et al.
A drastically different idea is represented by the MOND paradigm (Milgrom 1983) that modifies the standard laws of dynamics at low accelerations (weak gravitational fields) rather than assuming non-baryonic DM.Several a-priori predictions of MOND have been confirmed by later observations as reviewed by Sanders & McGaugh (2002), Famaey &McGaugh (2012), andMcGaugh (2020).The construction of a MOND cosmology remains a tall order (McGaugh 2015), but the recent relativistic MOND theory of Skordis & Zlośnik (2020) appears promising, being able to reproduce the power spectrum of the Cosmic Microwave Background as good as ΛCDM.
The relativistic theory of Skordis & Zlośnik (2020) reduces to the non-relativistic modified-gravity theory of Bekenstein & Milgrom (1984), violating the Strong Equivalence Principle (SEP) of GR: the internal dynamics of a self-gravitating body may be affected by external gravitational fields, beyond usual tidal forces.More specifically, these theories violate Local Positional Invariance (LPI) for gravitational experiments, which differentiates the SEP from the less stringent (but well tested) Einstein Equivalence Principle (EEP), containing the Weak Equivalence Principle, Lorentz Invariance, and the LPI for non-gravitational experiments only (Will 2014).
The radial acceleration relation (RAR) is of particular importance in the DM vs MOND debate (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017).This empirical relationship links the observed centripetal acceleration g obs (R) = V 2 rot (R)/R in galaxies to the expected Newtonian acceleration g bar (R) = V 2 bar (R)/R from the observed baryonic matter distribution: where ν 0 (z) is an empirical fitting function and g † is an acceleration scale.In ΛCDM the RAR must arise from the haphazard process of galaxy formation (Di Cintio & Lelli 2016;Desmond 2017;Navarro et al. 2017;Keller & Wadsley 2017) and g † is an emergent scale that may (Ludlow et al. 2017) or may not (Tenneti et al. 2018) appear in cosmological simulations.In MOND g † is a new universal constant of Nature indicated as a 0 (Milgrom 1983), while the function ν 0 (g bar /a 0 ) interpolates between the classic Newtonian regime g obs = g bar at high accelerations and the Milgromian regime g obs = √ g bar a 0 at low accelerations.While the extrapolation of Eq. (1) to large radii implies asymptotically flat rotation curves for isolated galaxies, MOND modified-gravity (Bekenstein & Milgrom 1984) predicts that galaxies in strong external fields should display a weak but distinctive decline in their outer rotation curves.This peculiar feature, linking the internal dynamics on scales smaller than 100 kpc with the cosmological environment on scales of a few Mpc, can be used to distinguish between modified gravity in MOND and standard gravity with DM.Signatures of this external field effect (EFE) have been searched for in rotationally-supported galaxies (Haghi et al. 2016;Wu & Kroupa 2015;Lelli et al. 2015) without conclusive and unambiguous evidence.
The EFE has also been investigated in pressuresupported stellar systems.Dwarf satellites of the Andromeda galaxy revealed some EFE signatures as predicted and tested by McGaugh & Milgrom (2013a,b), but the possibility of tidal interactions and out-ofequilibrium dynamics complicates the interpretation (e.g.McGaugh & Wolf 2010;Lelli et al. 2017).Several authors (Famaey, McGaugh & Milgrom 2018;Kroupa et al. 2018;Müller, Famaey & Zhao 2019;Haghi et al. 2019) proposed MOND models incorporating the EFE to explain unexpectedly low stellar velocity dispersions of a few ultra-diffuse galaxies.Globular clusters (GCs) of the Milky Way are dynamical systems subjected to external fields.MONDian kinematics for the GCs were predicted (Baumgardt, Grebel & Kroupa 2005;Haghi et al. 2009Haghi et al. , 2011)), but analyses of the observed data did not result in unambiguous signatures of the MOND EFE (Jordi et al. 2009;Frank et al. 2012).
Wide binary stars have also been used to test MOND and the EFE, with conflicting results (Hernandez et al. 2012;Pittordis & Sutherland 2019;Hernandez et al. 2019).In particular, wide binary stars from GAIA DR2 have been used to argue both for (Pittordis & Sutherland 2019) and against (Hernandez et al. 2019) the presence of the EFE, and further studies are required to provide conclusive evidence.
Here we report a robust EFE detection in rotationally supported galaxies using two complementary approaches: (1) focusing on individual galaxies where the external gravitational field is exceptionally large, (2) studying weak systematic deviations from the RAR driven by the mean gravitational field of the Local Universe.Throughout we take g † = 1.2 × 10 −10 m s −2 (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017) and use the notation x ≡ log 10 (g bar /m s −2 ) and y ≡ log 10 (g obs /m s −2 ).

DATA AND METHODOLOGY
2.1.The SPARC database The SPARC database (Lelli, McGaugh & Schombert 2016) contains 175 rotationally-supported galaxies in the nearby Universe 1 .These galaxies have stellar masses ranging from M 10 11 M to M 10 7 M and cover all Hubble types of late-type, star-forming galaxies, including low-surface brightness disk galaxies.The database provides the observed rotation velocities (V obs ) from spatially resolved HI observations and the Newtonian circular velocities from the observed distribution of stars and gas.The latter include the stellar disk contribution (V disk ) and (if present) the bulge contribution (V bul ) for a baseline mass-to-light ratios of unity, as well as the gas contribution (V gas ) for a total-to-hydrogen mass ratio of 1.33.For convenience, in this paper we redefine V gas for a total-to-hydrogen mass ratio of unity.The reported velocity V obs of a galaxy is tied to the reported inclination i obs .If the inclination is changed to i, the rotation velocity becomes The circular velocity due to the baryonic mass distribution depends on the galaxy distance D and is given by (3) where D ≡ D/D obs with D obs being the fiducial distance.In Eq. ( 3), Υ disk and Υ bul are the mass-to-light ratios of the disk and the bulge in units of the solar value M /L at 3.6µm, while Υ gas is the ratio of the total gas mass to the HI mass.When the SPARC database was published, this ratio was assumed to be 1.33 to account for the cosmic abundance of Helium from big bang nucleosynthesis.Here we consider the small amounts of Helium and metals formed via stellar nucleosynthesis during galaxy evolution (McGaugh, Lelli & Schombert 2020), so that Υ gas = X −1 where X is a function of stellar mass (M ): with M 0 = 1.5 × 10 24 M and α = 0.22.We do however allow the possibility of varying Υ gas from X −1 to consider the uncertainties in the HI flux, gas disk geometry, and the gas mass to HI mass ratio.In some cases V gas is negative at small radii, representing the fact that the Newtonian gravitational field is not oriented towards the center when a large fraction of the gas disk lies in the outer regions.To account for the cases of negative V gas we write V gas |V gas | rather than V 2 gas in the last term of Equation (3), although this detail has negligible effects on our study.

The external field effect
Empirically, the observed centripetal acceleration (g obs = V 2 obs /R) is related to the Newtonian baryonic acceleration (g bar = V 2 bar /R) via the RAR ν 0 (g bar /g † ) of Eq. (1) with a free parameter g † (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017).In a MOND framework, g † = a 0 is a fundamental constant of Nature (Milgrom 1983) and Eq. ( 1) can be obtained by modifying either inertia (Newton's second law of dynamics) or gravity (the Poisson's equation) at the non-relativistic (x0, y0) Figure 1.The external field effect in the weak-field limit of the radial acceleration relation.Eq. ( 5) is overlaid on the RAR for various values of e in Eq. ( 5).Values of e > 0 correspond to the MOND EFE, while e < 0 is unphysical from the MOND point of view.Values of e ≈ 0.033 corresponds to the average prediction for 153 SPARC galaxies based on their gravitational environments (Desmond et al. 2018).The heat map shows the original SPARC mass models (Lelli, McGaugh & Schombert 2016) with fixed stellar mass-to-light ratios for the same galaxies.
level (Famaey & McGaugh 2012).In MOND modifiedinertia theories Eq. (1) holds exactly for any circular orbit (Milgrom 1994), while in MOND modified gravity theories holds only for highly symmetric mass distributions (such as spheres) and represents a first-order approximation for actual disk galaxies (Brada & Milgrom 1995).In all these scenarios, however, Eq. ( 1) is strictly valid only for isolated systems, when the external field effect (EFE) is negligible.
To build a general fitting function that approximates the EFE, we start from the nonlinear MOND modified Poisson's equation (Bekenstein & Milgrom 1984) in the one-dimensional case.If we assume a uniform external gravitational field g ext (Famaey & McGaugh 2012) and the so-called Simple interpolating function (IF) (Famaey & Binney 2005), we have with Chae et al.
where z ≡ g bar /g † , A e ≡ e(1+e/2)/(1+e), B e ≡ (1+e), and e ≡ g ext /g † .For e = 0, ν e (z) is reduced to the Simple IF ν 0 (z) = 1/2 + 1/4 + 1/z.Equation ( 6) is based on the footnote to Eq. ( 59) of Famaey & McGaugh (2012), but we corrected a small typo and rearranged it.Note here that the Simple IF allows the convenient analytic form of Equation ( 6) with e > 0 while it is only subtly different (Chae et al. 2019) in the (EFE-irrelevant) high acceleration limit from the function used by Mc-Gaugh, Lelli & Schombert (2016) and Lelli et al. (2017) to fit the SPARC galaxies.Our results on the EFE detection are not affected by the choice of the Simple IF.
Then, the expected circular velocity is given by Although ν e (z) (Eq.6) is based on idealized assumptions, it captures the basic feature of the EFE: a systematic downward deviation from ν 0 (z) (Eq. 1) when e > 0 as z → 0. Eq. ( 6) also allows for upward deviations when e < 0, which seem unphysical but may be preferred by the data at the empirical level.These features are illustrated in Fig. 1.MOND with the EFE predicts that the RAR must be a family of functions rather than a universal function.This also means that if galaxies in different environments are tried to be fitted with a single functional form of Eq. (1), then there will arise some small intrinsic scatter of g † due to the EFE.Most importantly, regardless of its MOND origin, Eq. ( 6) may be considered a mere fitting function that improves over Eq. ( 1) by adding the free parameter e, which has no a-priori knowledge of the external gravitational field in which galaxies reside.

MCMC simulations
In our Bayesian analysis the posterior probability of parameters β = {β k } is defined by where Pr(β k ) is the prior probability of parameter β k and χ 2 is given by is the reported error of V obs (R j ) for the reported inclination i obs .As in earlier studies of the RAR using SPARC galaxies (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017), we use only 153 galaxies with i obs ≥ 30 • and Q ≤ 2 (a quality cut on the rotation curve).The parameters β in Eq. ( 8) are given by β = {Υ disk , Υ bul , Υ gas , D, i, e} for the case of using Eq. ( 5) with a fixed g † = 1.2 × 10 −10 m s −2 .The priors on these parameters are summarized in Table 1.The mean values and standard deviations of Υ disk and Υ bul are motivated by state-of-the-art stellar population synthesis models for star-forming galaxies (Schombert, McGaugh & Lelli 2019).The mean value of Υ gas is given by Eq. ( 4), while the standard deviation is motivated by the typical error on the HI flux calibration, but it could also represent variations in the assumed gas disk thickness and/or the mean gas-to-HI mass ratio.The mean values and standard deviations of D and i consider the baseline SPARC values and their fiducial errors.For e we adopt an uninformative uniform prior covering a reasonably broad range.
The posterior probability density functions (PDFs) of the model parameters are derived from MCMC simulations through the public code emcee (Foreman-Mackey et al. 2013).These simulations represent an extension to the previous SPARC analysis (Li et al. 2018) including the EFE parameter e.We choose N walkers = 10000 and N iteration = 6000.We discard models up to N iteration = 500 and thin the rest by a factor of 50 as the auto-correlation lengths for the parameters are < 100.The posterior PDFs of x = log 10 g bar (R) and y = log 10 g tot (R) follow from the posterior PDFs of the parameters i, log 10 D, log 10 Υ disk , log 10 Υ bul , and log 10 Υ gas .

The environmental gravitational field
We estimate the environmental gravitational field g env due to the large-scale distribution of matter at the positions of the SPARC galaxies.We perform this calculation within the standard ΛCDM context (Desmond et al. 2018).A similar calculation is not feasible in a MOND context due to the strong non-linearities in the theory and the lack of a proper MOND cosmology.The ΛCDM calculation, however, is a good first-order approximation for MOND and other modified gravity theories (Desmond et al. 2018), up to some systematic uncertainty due to the unknown relation between g env in these theories.We use g env primarily for the purpose of picking out extreme cases with exceptionally high or low g env (which should remain true in a relative sense in any cosmological scenario) and to check that the maximumlikelihood values of e from fitting Eq. ( 5) are sensible in an order-of-magnitude fashion.
Our calculation of g env starts with the total dynamical masses of the galaxies in the all-sky 2M++ survey (Lavaux & Hudson 2011) using abundance matching.We then use N-body simulations in ΛCDM to populate the surrounding regions with halos hosting galaxies too faint to be recorded in 2M++, using statistical correlations between halo abundances and properties of the galaxy field.Finally, we add mass in long-wavelength modes of the density field according to the inferences of the BORG algorithm (Lavaux & Jasche 2016) applied to the 2M++ catalog.We use the final density field to calculate a posterior distribution for g env at the position of each SPARC galaxy, fully propagating uncertainties in the input quantities.We define e env ≡ g env /g † , and find values in the range 0.01 e env 0.1 with a mean of 0.033 among the SPARC galaxies: typical values are in the range 0.02 − 0.05.

Individual galaxies
The estimated values of the environmental gravitational field strength e env ( § 2.4) span almost one order of magnitude, thus there is sufficient dynamic range to check whether the rotation-curve shapes depend or not on the large-scale environment.Among the SPARC galaxies whose rotation curves (RCs) reach g obs < g † , NGC5033 and NGC5055 live in exceptionally dense environments with e env ≈ 0.1, while NGC1090 and NGC6674 are exceptionally isolated with e env ≈ 0.01.The former two represent "golden galaxies" for the EFE to be detected, while the latter two are control targets for the null detection.
We fit the RCs using the EFE-incorporated RAR fitting function (Eq 5) with a free external field g ext parameterized by e = g ext /g † ( § 2.2).The case of e = 0 implies flat RCs and reduces exactly to the original RAR (Eq.1).Fig. 2 shows the MCMC results for the RCs of four galaxies in the two extreme environments.The 'corner' plots showing the posterior PDFs of the parameters for these galaxies can be found in the appendix.
For NGC5055, the detailed shape of the RC is very well fitted with a positive e but poorly fitted with e = 0. We find e = 0.054 ± 0.005: this is an 11σ detection.Remarkably, this value is consistent within 2σ with e env = 0.094 +0.089 −0.022 that is independently determined from the large-scale environment.The Bayesian information criterion (BIC ≡ −2 ln L max + k ln N where L max is the maximum likelihood, k is the number of free parameters, and N is the number of the fitted rotation velocities) for e = 0 relative to the free e case is ∆BIC = 144, indicating very strong evidence for e > 0 based on the conventional criterion of ∆BIC > 10 for strong evidence.
For NGC5033, the overall fit is also improved by freeing up e since ∆BIC = 83.9.We find e = 0.104 +0.013 −0.012 .This is an 8σ detection, in excellent agreement with e env = 0.102 +0.086 −0.021 from the large-scale environment.The observed properties of this galaxy, however, are not as robust as those of NGC5055.The rotation velocities at R < 60 arcseconds (about 5 kpc) are probably underestimated due to beam-smearing effects in the HI data, although our results on e are not affected by these data points.Moreover, while the distance of NGC5055 is robust because it is based on the tip magnitude of the red giant branch (D = 9.90 ± 0.30), that of NGC5033 is very uncertain because it is estimated using Hubble flow models (D = 15.70 ± 4.70).Interestingly, our MCMC result for NGC5033 predicts a relatively large distance (D = 23.5 +2.0 −1.8 Mpc) with e > 0 but a low one (D = 13.0 +0.7 −0.6 Mpc) with e = 0. Hence, future observations can provide a key independent test.
In striking contrast to the highest e env sample, the galaxies in the lowest e env sample show no strong evidence for e > 0 based on ∆BIC (or any other widelyused statistic).These two galaxies are similar to the golden galaxies in morphology, mass, and size.The only noticeable difference is that they are unusually isolated.
The fitted e values are consistent with the independent e env values within about 2σ.

Statistical approach
Since the EFE has subtle effects on rotation-curve shapes, positive values of e are detected with high statistical significance only in individual galaxies where e env is exceptionally large (like NGC5055 and NGC5033).The EFE, however, should also imprint a statistical signature in the low-acceleration portion of the RAR (see Fig. 1).

The systematic trend in the low-acceleration portion of the RAR
We use 153 galaxies from the SPARC database ( § 2.1).Fig. 3 (top panels) shows the RAR for 2696 points having accuracy in V rot better than 10%.In the top left panel we first show the original SPARC mass models (Lelli, McGaugh & Schombert 2016) with fixed massto-light ratios at 3.6µm of Υ disk = 0.5M /L for the disk and Υ bulge = 0.7M /L for the bulge (Lelli et al. 2017).The MCMC mass models obtained here with varied mass-to-light ratios ( § 2.3) are shown in the top right panel.We divide the data points into bins perpendicular to the best-fit curve assuming Eq. (1).Each data point (x, y) is projected onto the point (x 0 , y 0 ), so the orthogonal residual ∆ ⊥ encodes any possible systematic deviation from the e = 0 case.
The data show a small systematic deviation from Eq. (1) for x 0 −11.This trend is present, though weakly, in the original SPARC mass models with fixed mass-to-light ratios for the disks and bulges.The MCMC models in the middle column show a stronger effect.The systematic deviation is weak in absolute   5) (left panels).The colored bands show the 1σ confidence limits for the rotation curve (red) and the separate contributions of gas disk (green), stellar disk (blue), and stellar bulge (orange) if present.For the "golden galaxies" subjected to the strongest environmental gravitational fields, the fit is improved dramatically with e > 0, resulting in 11σ and 8σ individual detections of the EFE.For the galaxies subjected to the weakest fields, the EFE is not detected as expected.In all cases, the fitted values of e are fully consistent with the independent values of eenv from the large-scale galaxy environment within ∼ 2σ.∆BIC indicates evidence by the Bayesian Information Criterion.
terms (0.05-0.08 dex for the lowest x 0 bin) but at least 4 times larger than the bootstrap error of the median in the bin.This demonstrates that Eq. ( 1) does not fully capture the trends in the RAR.Introducing e as an additional free parameter, we obtain a better fit and find e ≈ 0.02 − 0.04 in close agreement with the independent estimate of e env 0.033 from an all-sky galaxy catalog ( § 2.4).

The statistical detection of the EFE
The systematic trend in the RAR also implies that the fitted individual values of e of Eq. ( 5) will be systematically displaced from the non-EFE case e = 0.The posterior PDFs of e are quite broad with a typical standard deviation of ∼ 0.04 (see the appendix for examples).Nevertheless, the statistical distribution of the fitted values will have a signature.Because e was allowed to take any value (positive or negative), this distribution provides a blind test of MOND EFEs ( § 2.2).
Fig. 4 shows the distribution of the orthogonal residual ∆ ⊥ and the fitted median value of e from the MCMC simulations with Eq. ( 5).From Fig. 1 it is expected that data points at high enough accelerations do not have any sensitivity to e. Indeed, for data points with −10.3 < x 0 < −9 the distribution of ∆ ⊥ gives a null result.Similarly, for galaxies with −10.3 < x 0 < −9, the distribution of e MLE gives a null result.
Data points at low enough accelerations will have sensitivity to e and distributions with non-zero mean value are expected from Fig. 3.For data points with x 0 < −11.3 the distribution of ∆ ⊥ has a mean of −0.061 ± 0.008 (a bootstrap error) which is statistically significant at more than 7σ.For a much larger number of data points with x 0 < −10.3, ∆ ⊥ has a smaller deviation of −0.035 ± 0.003, but the statistical significance of the deviation is more than 11σ.
Fig. 5 shows individual e values and their uncertainties for a subset of 113 galaxies with median x 0 < −10.3.Due to the large uncertainties on e, some galaxies can occasionally return negative values.However, the median value of e is 0.052 ± 0.011 (bootstrap error), which represents ≈ 5σ detection of positive e.This value is statistically consistent with the median environmental gravitational field for these galaxies ( e env = 0.034 ± 0.001 (bootstrap error)).Furthermore, based on the robust binomial statistic with equal probabilities for e > 0 and e < 0, 78 cases of e > 0 out of 113 is 4σ away from the expected mean of 56.5 cases.
Fig. 6 further shows the distribution of the individual difference e − e env .It has a broad distribution due to the large uncertainty in e but is clearly consistent with zero: e − e env = 0.011 ± 0.013.It is intriguing that the mere fitting parameter e returns, on average, the same value of the mean environmental gravitational field of the nearby Universe, computed in a fully independent way.

Statistical properties of the posterior parameters of the galaxies
Fig. 7 shows the distribution of the parameters from the MCMC simulations with Eq. ( 5) for all 153 selected galaxies.The distribution of the distances is consistent with the SPARC reported values with an rms scatter of 0.02 dex (5 percent).This is smaller than typical measurement uncertainties of ∼ 14 percent.The posterior inclination angles are also consistent with the SPARC reported values with an rms scatter of 2.1 • , smaller than typical measurement uncertainties of ∼ 4 • .The distributions of the mass-to-light ratios (Υ disk and Υ bul ) for the disk and the bulge are consistent with the estimates from infrared studies, i.e., Υ disk = 0.5 M /L and Υ bul = 0.7 M /L , with an rms scatter of 0.14 dex.If anything Υ bul might be 0.6, a little smaller than 0.7.Finally, the distribution of Υ gas is in excellent agreement with X −1 from Eq. ( 4), giving a mean value of 1.38 ± 0.04 which is intermediate between a metal-poor dwarf galaxy with X −1 = 1.34 and a metal-rich giant spiral with X −1 = 1.42.

Analysis of possible systematic effects
One may wonder whether the systematic deviations from Eq. ( 1) are due to some systematic uncertainties.There are three main observational effects that may systematically affect the low acceleration portion of the RAR: galaxy distances, the thickness of the gas disk, and possible variations of M /L in the stellar disk with radius.To mitigate the first two uncertainties, the left panel of Fig. 8 considers data points from galaxies with accurate distances based on the tip magnitude of the red giant branch, Cepheids, or Supernovae (Lelli, McGaugh & Schombert 2016), as well as low gas contributions (f gas = M gas /M bar < 0.4).Compared with Fig. 3 in the main manuscript, it is clear that the scatter is smaller and the median trend is consistent with the full dataset.
The thickness of the gas disk is a concern because the EFE is detected in the galaxy outskirts, where the gas contribution becomes non-negligible or even dominating in some cases.Recent studies (Bacchini et al. 2019) suggest that gas disks may become thicker at large radii: this would systematically decrease V gas , hence g bar , moving points to the left of the RAR.Therefore, we repeat the MCMC fits considering gas disks that are three times thicker than assumed in the SPARC database.This is a very extreme scenario because not all galaxies will have such thick gas disks.Our goal is simply to provide an upper bound on the possible impact of this effect.Fig. 8 (middle panel) shows that there is still a significant systematic deviation from Eq. ( 1) even when we consider very thick gas disks.
Negative gradients of M /L with R could also systematically decrease V disk , hence g bar , moving points to the left of the RAR.While we are treating the bulge separately in the most massive spirals (Sa to Sb), the In the top panels, the Newtonian acceleration from the baryons g bar is plotted against the observed acceleration g obs for a sample of 153 SPARC galaxies.The typical error bars are indicated in the bottom.The data are fitted using Eq. 1 (red solid curve with g † = 1.2 × 10 10 m s −2 ) corresponding to e = 0, and using the additional free parameter e accounting for the EFE (green dashed curve).The black dots show the median values within the bins orthogonal to Eq. (1) (red dotted lines).The inset illustrates how orthogonal residuals are calculated.The bottom panels show the orthogonal residuals versus x0: the deviation at x0 < −11 represents a statistical detection of the EFE.The inset zooms into this interesting region.The left column shows the original SPARC mass models with fixed stellar mass-to-light ratios, while the right column shows the MCMC results with varied stellar mass-to-light ratios and considering the EFE.In both cases, the fitted e value is remarkably similar to eenv ≈ 0.033 from the large-scale mass distribution in the nearby Universe.The median value of e = 0.052 ± 0.011 is statistically consistent with eenv = 0.034 ± 0.001 (see also Supplementary Materials).The individual galaxies considered in Fig. 2 are indicated: for the golden galaxies at exceptionally high eenv values, e is significantly different from zero at 8σ (NGC5033) and 11σ (NGC5055).Big dots indicate galaxies with accurate distances.
stellar disk may potentially display a radial variation of its stellar populations.At 3.6 µm these variations have a relatively weak effect (Schombert, McGaugh & Lelli 2019), but we nevertheless repeat the MCMC fits considering a linear decrease in Υ disk by a factor of 2 from the center to the outermost observed radius.Again, this is an extreme scenario since most stellar disk are likely not showing such strong radial gradients in Υ disk .Fig. 8 (right panel) shows that there is still a significant systematic deviation from Eq. (1).

Comparison with previous results
Only a few attempts have been made so far to detect the EFE from the RCs of galaxies (Haghi et al. 2016;Wu & Kroupa 2015).In particular, Haghi et al. (2016) considered the RCs of 18 galaxies taken from the literature available at that time.These galaxies are known to have relatively nearby massive neighbors.Eleven of them are also included in our sample of 153 galaxies The distribution of e derived from MCMC fits to the rotation curves is compared with that estimated from the observed environments of the galaxies (eenv).There is good agreement up to the large uncertainties on the fitted values.
studied here.They are DDO 154, IC 2574, NGC 2998, NGC 3198, NGC 3521, NGC 3769, NGC 4100, NGC 4183, NGC 5033, NGC 5055, and NGC 5371. Haghi et al. (2016) obtained values of e ranging from about 0.1 to 0.6 with a median of ∼ 0.3 and a typical uncertainty of ∼ 0.1 for these 11 galaxies.Their values are systematically higher than our values ranging from about −0.1 to 0.3 with a median of ∼ 0.075 and a typical uncertainty of ∼ 0.04.This is primarily due to the fact that the disk models of Haghi et al. (2016) are based on a baryonic mass profile that declines more slowly than observed at large radii, requiring a larger EFE in the MOND context (a deficit of DM in the ΛCDM context).
There have also been indications of the EFE in pressure-supported galaxies (McGaugh & Milgrom 2013a,b;Famaey, McGaugh & Milgrom 2018;Kroupa et al. 2018).Pressure-supported galaxies are analyzed through their observed line-of-sight velocity dispersions.Because their stellar orbits are complex and not observed directly, a robust kinematic analysis to infer the EFE is challenging.However, McGaugh & Milgrom (2013a,b) have found that the observed velocity dispersions of the dwarf galaxies of the Andromeda galaxy are consistent with a MOND theory with EFE.More recently, galaxies that appeared to have too low observed velocity dispersions and thus lack dark matter in the ΛCDM context (van Dokkum et al. 2018(van Dokkum et al. , 2019) ) may well be explained by the MOND EFE (Famaey, McGaugh & Milgrom 2018;Kroupa et al. 2018;Müller, Famaey & Zhao 2019;Haghi et al. 2019). .Fitted parameters for the 153 SPARC galaxies.The distributions of the fitted parameters from our MCMC simulations using Eq. ( 5) are compared with the SPARC measured or assumed values.
Galaxies of similar properties but subjected to different external gravitational fields show noticeably different rotation-curve behaviors at large radii (i.e. at very low accelerations).Two galaxies in the strongest environmental fields show declining RCs in the outer parts, while two similar galaxies in the weakest environmental fields have flat RCs.The connection between internal dynamics and large-scale environment is corroborated by a statistical analysis of the entire SPARC sample.At accelerations 10 times lower than g † , the RAR is not fully described by a simple function of g bar /g † (Eq. 1) but requires an EFE-incorporated generalized function with an additional free parameter e (Eq.6).Moreover, rotation-curve fits with Eq. ( 6) give a mean value of e that is indistinguishable from the mean environmental gravitational field at the location of SPARC galaxies, computed in a fully independent fashion from the average distribution of mass in the nearby Universe.These results are summarized in Figs. 3 and 5.Note that these results of fitting Eq. ( 6) to RCs are fully empirical, independent of any theoretical interpretation.
Can these results be explained in the standard ΛCDM framework?For the two golden massive galaxies subjected to strong large-scale gravitational field g env , declining RCs are observed over a radial range of about 30 -50 kpc, which are less than ∼ 15% of the virial radius of the DM halo.Clearly, this is not the decline that should occur in the outer parts of ΛCDM halos, where the density profile decreases as r −3 , since we are probing the inner parts of the halo where the density profile goes approximately as r −2 , leading to flat RCs.
Thus, the only remaining option is represented by tidal forces.We calculated the expected tidal radii in ΛCDM using the formalism of King (1962), taking the source of the tidal field to be the nearest 2M++ galaxy to the SPARC galaxy in question.We assume the source and test galaxies to have NFW (Navarro, Frenk & White 1997) halos following the M -M vir relation of Kravtsov et al. (2018) and the M vir −concentration relation of Diemer & Kravtsov (2015).We find the tidal radii to be much larger than the last measured points of the RCs, so the galaxies themselves are effectively shielded against large-scale tides.
The agreement between the MOND fitting parameter e (Eq.5) and the environmental gravitational field e env is an unpredicted result from the ΛCDM point of view.In principle, the baryon plus DM combination can combine to produce a declining rotation curve within tens of kpc as found here (i.e. e > 0).For that matter, however, there is no a priori reason that the degree of declin-   3 in the main paper.The left panels show a subset of data points with sub-dominant gas contribution fgas < 0.4 and accurate distance measurements.The middle panels show the MCMC results assuming 3× thicker gas disks for all galaxies.The right panels show the MCMC results assuming a radial M /L gradient in the stellar disks of all galaxies.Our conclusions on the EFE detection hold in all cases.
ing must agree with the strength of the environmental gravitational field.There could have been an order-ofmagnitude difference between e and e env .Yet, we are seeing an interesting coincidence between the two.
Moreover, a downward deviation in the RAR near a tenth of g † is not predicted by current ΛCDM state-ofthe-art simulations or semi-analytical models (Di Cintio & Lelli 2016;Desmond 2017;Navarro et al. 2017;Keller & Wadsley 2017;Tenneti et al. 2018) with some predicting the opposite trend (Ludlow et al. 2017;Fattahi et al. 2018;Garaldi et al. 2019).To the best of our knowledge, there is no reported scenario in which the DM-baryon coupling in the outskirts of the disks depends on the external gravitational field from the large-scale galaxy environment in the manner found here.
The empirical evidence is fully consistent with the EFE predicted by MOND modified gravity (Bekenstein & Milgrom 1984).More generally, our results suggest a violation of the SEP in rotationally-supported galaxies.While in GR the internal dynamics of a gravitationallybound system is not affected by a uniform external field, our analysis indicates that external fields do impact the internal dynamics.Our results are encouraging for modified gravity as an alternative (or modification) to the DM hypothesis and the standard ΛCDM cosmological model.They also highlight the path for future theoretical investigations of relativistic theories of gravity beyond GR (see, e.g., Skordis & Zlośnik 2020), possibly leading to a new cosmological model.

CONCLUSIONS
In this paper we provide observational evidence for the existence of the EFE (or a phenomenon akin to it) predicted by MOND modified gravity (Bekenstein & Milgrom 1984).We use accurate rotation curves and mass models from the SPARC database (Lelli, McGaugh & Schombert 2016) and detect the EFE in three separate ways: 1.The EFE is individually detected in "golden" galaxies subjected to exceptionally strong external gravitational fields.The detection is highly significant (11σ in NGC5055 and 8σ in NGC5033) and the best-fit values of the external gravitational fields are fully consistent with the independent estimates from the large-scale distribution of mass at the galaxies' location.Conversely, the EFE is not detected in control galaxies residing in the weakest external gravitational fields, as expected.
2. The EFE is statistically detected at more than 4σ through a blind test using 153 SPARC galaxies.The mean value of the external gravitational field among the SPARC galaxies is again consistent with the independent estimate from the average distribution of mass in the nearby Universe.
3. The EFE also manifests as a small ( 0.05 dex), downward deviation from the empirical RAR occurring around 0.1g † .This behavior is not predicted by any of the existing galaxy-formation models in ΛCDM that were proposed to "naturally" reproduce the RAR.In contrast, this downward deviation is predicted by the MOND modified gravity at the right acceleration scale.
Our results suggest a breakdown of the SEP: the internal dynamics of a gravitational system in free-fall is affected by a uniform external gravitational field.This sheds new light on the dark-matter problem and paves the way for relativistic theories of modified gravity in the weak-field regime of gravity g 10 −10 m s −2 .

Figure 2 .
Figure2.Detection of the EFE in individual galaxies.The observed rotation curves (points with errorbars) are fitted using Eq.(1) with no EFE (right panels) and a generalized equation considering the EFE (Equation5) (left panels).The colored bands show the 1σ confidence limits for the rotation curve (red) and the separate contributions of gas disk (green), stellar disk (blue), and stellar bulge (orange) if present.For the "golden galaxies" subjected to the strongest environmental gravitational fields, the fit is improved dramatically with e > 0, resulting in 11σ and 8σ individual detections of the EFE.For the galaxies subjected to the weakest fields, the EFE is not detected as expected.In all cases, the fitted values of e are fully consistent with the independent values of eenv from the large-scale galaxy environment within ∼ 2σ.∆BIC indicates evidence by the Bayesian Information Criterion.

Figure 3 .
Figure3.EFE detection in the low-acceleration portion of the RAR.In the top panels, the Newtonian acceleration from the baryons g bar is plotted against the observed acceleration g obs for a sample of 153 SPARC galaxies.The typical error bars are indicated in the bottom.The data are fitted using Eq. 1 (red solid curve with g † = 1.2 × 10 10 m s −2 ) corresponding to e = 0, and using the additional free parameter e accounting for the EFE (green dashed curve).The black dots show the median values within the bins orthogonal to Eq. (1) (red dotted lines).The inset illustrates how orthogonal residuals are calculated.The bottom panels show the orthogonal residuals versus x0: the deviation at x0 < −11 represents a statistical detection of the EFE.The inset zooms into this interesting region.The left column shows the original SPARC mass models with fixed stellar mass-to-light ratios, while the right column shows the MCMC results with varied stellar mass-to-light ratios and considering the EFE.In both cases, the fitted e value is remarkably similar to eenv ≈ 0.033 from the large-scale mass distribution in the nearby Universe.

Figure 4 .Figure 5 .
Figure 4. Distributions of ∆ ⊥ and e from fitting Eq. (5) to the SPARC galaxies.The top panels show the distributions of orthogonal residuals ∆ ⊥ for three acceleration bins from the MCMC results shown in the middle column of Fig.3.The mean of the distribution is displaced from zero for lower acceleration bins, indicating declining RCs.The bottom panels show the distributions of the e values fitted to the individual galaxies binned by the median values of x0 within the galaxies.As expected, for the galaxies in the high acceleration bin (−10.3 < x0 < −9.0), the data do not have any sensitivity to e and so the distribution has a mean of ∼ 0. For lower acceleration bins the distributions are shifted to positive e with high statistical significance, indicating a preference for the EFE.The broad distributions are due to the broad individual posteriors on e.

Figure 6 .
Figure6.Comparison of external field strength estimates from kinematics versus analyzing the galaxies' environments.The distribution of e derived from MCMC fits to the rotation curves is compared with that estimated from the observed environments of the galaxies (eenv).There is good agreement up to the large uncertainties on the fitted values.
Figure7.Fitted parameters for the 153 SPARC galaxies.The distributions of the fitted parameters from our MCMC simulations using Eq.(5) are compared with the SPARC measured or assumed values.

Figure 8 .
Figure8.Testing systematic uncertainties in the EFE detection.This figure has the same format as Fig.3in the main paper.The left panels show a subset of data points with sub-dominant gas contribution fgas < 0.4 and accurate distance measurements.The middle panels show the MCMC results assuming 3× thicker gas disks for all galaxies.The right panels show the MCMC results assuming a radial M /L gradient in the stellar disks of all galaxies.Our conclusions on the EFE detection hold in all cases.

Figure 13 .
Figure 13.Parameter corner plot for NGC2955.The posterior PDFs of the parameters for a normal galaxy NGC2955 produced from MCMC simulations using Eq.(5).

Figure 14 .
Figure 14.Parameter corner plot for NGC6195.The posterior PDFs of the parameters for a normal galaxy NGC6195 produced from MCMC simulations using Eq.(5).

Table 1 .
Summary of prior constraints on the model parameters