Testing the Strong Equivalence Principle. II. Relating the External Field Effect in Galaxy Rotation Curves to the Large-scale Structure of the Universe

Theories of modified gravity generically violate the strong equivalence principle, so that the internal dynamics of a self-gravitating system in freefall depends on the strength of the external gravitational field (the external field effect). We fit rotation curves (RCs) from the SPARC database with a model inspired by Milgromian dynamics (MOND), which relates the outer shape of an RC to the external Newtonian field from the large-scale baryonic matter distribution through a dimensionless parameter e N. We obtain a > 4σ statistical detection of the external field effect (i.e. e N > 0 on average), confirming previous results. We then locate the SPARC galaxies in the cosmic web of the nearby universe and find a striking contrast in the fitted e N values for galaxies in underdense versus overdense regions. Galaxies in an underdense region between 22 and 45 Mpc from the celestial axis in the northern sky have RC fits consistent with e N ≃ 0, while those in overdense regions adjacent to the CfA2 Great Wall and the Perseus−Pisces Supercluster return e N that are a factor of two larger than the median for SPARC galaxies. We also calculate independent estimates of e N from galaxy survey data and find that they agree with the e N inferred from the RCs within the uncertainties, the chief uncertainty being the spatial distribution of baryons not contained in galaxies or clusters.

1. INTRODUCTION General relativity and its Newtonian limit are unique in that the internal gravitational dynamics of a selfgravitating system is unaffected by a constant external gravitational field (Will 2014).This strong equivalence principle (SEP) is generally violated by Milgromian dynamics (MOND; Milgrom 1983), which alters Newton's Corresponding author: Kyu-Hyun Chae, Harry Desmond KHC: chae@sejong.ac.kr, kyuhyunchae@gmail.comHD: harry.desmond@physics.ac.ukFL: federico.lelli@inaf.itSSM: ssm69@case.eduJMS: schombe@gmail.comlaws at low accelerations to explain the dynamics of galaxies without resorting to dark matter (DM).Specifically, MOND violates local position invariance for gravitational experiments because it is non-linear even in the non-relativistic regime, so an external acceleration is expected to affect the internal dynamics of a system.This is known as the external field effect (EFE).The most prominent observable consequence of the MOND EFE is that, whereas isolated galaxies should have asymptotically flat rotation curves (RCs), galaxies in dense environments should have declining RCs at large radius.Thus, the EFE produces a smoking gun signature of MOND as contrasted with Newtonian dynamics plus DM, as already emphasized in the first paper proposing MOND (Milgrom 1983).
Non-relativistic field theories of gravity embodying the MOND EFE include the AQUAdratic Lagrangian (AQUAL) theory (Bekenstein & Milgrom 1984) and the quasi-linear MOND (QUMOND) theory (Milgrom 2010).Relativistic theories, which reduce to either AQUAL or QUMOND in the non-relativistic limit, are under active development (see Famaey & McGaugh 2012;Milgrom 2014 for reviews).In particular, the proposal of Skordis & Zlośnik (2020) appears particularly promising as it can reproduce the angular power spectrum of the cosmic microwave background as well as the linear matter power spectrum.
There have been several attempts to test the EFE empirically (McGaugh & Milgrom 2013a,b;Hees et al. 2016;Haghi et al. 2016;Wu & Kroupa 2015;Famaey et al. 2018), all of which find better agreement with galaxy dynamics when the EFE is included.Recently, (Chae et al. 2020b, hereafter Paper I) studied the EFE in 153 rotationally supported galaxies from the Spitzer Photometry and Accurate Rotation Curves (SPARC) database (Lelli, McGaugh & Schombert 2016).They presented both individual and statistical EFE detections from fitting RCs with an AQUAL-inspired function, which contains the MOND external field strength g e as a free parameter (see also Chae et al. 2021).They also showed that the median g e agrees with an independent estimate of the external field strength from the galaxies' largescale environments.This is a surprising result from the DM perspective because the observed decline in the RCs occurs well within the tidal radius.Paper I, however, estimated the environmental external fields under the Λ cold dark matter (CDM) paradigm using DM halos as surrogates of the MOND effects; this was expedient because g e from the large-scale environment cannot be easily computed in AQUAL due to its non-linearity.
Here we use a formalism that quantitatively relates the decline of a galaxy's RC to the strength of the Newtonian external field g Ne at its position.As g Ne is sourced linearly, it can be calculated by summing the Newtonian fields from the large-scale distribution of baryons around the galaxy.Thus, the EFE can be tested directly with the observable Newtonian external field, and the baryonic content of the nearby Universe can be studied as a byproduct.The detailed calculation of the Newtonian external gravitational field produced by the baryonic large scale structure, for direct comparison with the fitted values, is the key advance of this work over Paper I.
In Section 2 we describe our formalism and RC-fitting method, with technical details given in Appendix A. In Section 3 we describe our method for estimating the environmental Newtonian fields, first with an approximate all-sky calculation (Section 3.1) and then with a more detailed calculation within the SDSS footprint (Section 3.2).We present our results in Section 4 and discuss their implications in Section 5.In Appendix B, we provide tables of the fitted parameters and the estimated Newtonian environmental field strengths.

FORMALISM AND FITTING METHOD
Our goal is to infer the external Newtonian field directly from the RCs of galaxies to compare it with the Newtonian field estimated from their cosmic environment, so we need an EFE model that is parameterized by e N ≡ g Ne /a 0 , where a 0 = 1.2 × 10 −10 m s −2 is the MOND acceleration scale (McGaugh, Lelli & Schombert 2016;Lelli et al. 2017;Li et al. 2018;Chae et al. 2020a;Li et al. 2021).In Paper I we used an AQUAL-based model that was parameterized by e ≡ g e /a 0 where g e is the MOND external field.In this model, the expected acceleration of a test particle (V 2 /R for circular motion) is given by: g = ν e (g N /a 0 )g N , where g N is the Newtonian internal gravitational field (sourced by the galaxy itself) and ν e (y) is given by Equation ( 6) of Paper I with y ≡ g N /a 0 (see also Appendix A).This equation is based on the simple interpolating function (Famaey & Binney 2005;Chae et al. 2019) and the analytic solution of the AQUAL equations in one dimension (Famaey & McGaugh 2012).
In the QUMOND formalism, one can analytically derive the ratio g/g N as a function of e N in one dimension.Incidentally, as shown in Appendix A, this ratio is equivalent to Equation (6) of Paper I with e 2 /(1 + e) = e N for non-negative values of e N .Thus, Equation (1) holds with the replacement of ν e (y) with the following where C eN ≡ 1 + |e N |/4 and D eN (y) ≡ 1 + |e N |/y.Note that e N < 0 is unphysical (similarly to e < 0), but we allow it because the measured RCs may prefer negative values empirically.For typical e N < 0.01 and the acceleration range probed by galactic RCs, the last term of Equation (2) mainly controls the EFE and scales almost linearly with e N / |e N |.Thus, we define the fitting parameter ẽ ≡ e N / |e N |.Our notations are summarized in Table 1.
Values of ẽ for SPARC galaxies may be obtained simply from e of Paper I through an analytically extended transformation such as ẽ = e/ √ 1 + e to include negative values.We prefer, however, to obtain values of ẽ directly from fitting Equation (2) to the RCs.Our procedure of obtaining probability density functions (PDFs) of ẽ and galactic parameters is the same as that of Paper I. The galactic parameters are mass-to-light ratio(s) of the disk (and the bulge if present), total gas-to-hydrogen mass ratio, distance to the galaxy normalized by the SPARC reported value, and inclination of the disk.
We also revise our selection to maximise the sample size.From the SPARC database of 175 rotationallysupported galaxies within distance D 120 Mpc, the RCs of 163 galaxies have good or acceptable quality (Q = 1 or 2).Out of these, 10 galaxies have low measured inclination angles of < 30 • .Paper I excluded these low-inclination galaxies for their analysis following Mc-Gaugh, Lelli & Schombert (2016) and thus used only 153 galaxies.Here we use all 163 galaxies including the lowinclination ones for the purpose of enlarging the sample and as our MCMC procedure allows the inclination angle to vary freely.However, including or excluding lowinclination galaxies does not affect our general results.
Examining all posterior PDFs, we find that for one galaxy (UGC 06787) the PDFs are double-peaked and the posterior distance is an order of magnitude larger than the SPARC reported value (D = 21 ± 5 Mpc).Its RC (from Noordermeer et al. 2007) has an unusual sinusoidal shape that is hard to fit in any context, either with MOND (Li et al. 2018) or with DM halos (Li et al. 2020).The galaxy has a strongly warped HI disk, so the RC shape depends on the details of the warp modeling (the adopted radial variation of position angle and inclination).Considering the unreliable value of the best-fit distance (and probably of ẽ), we drop this galaxy because we aim to compare the external Newtonian field with the location of galaxies on the large-scale structure of the Universe.Thus, our analysis is based on a starting sample of 162 SPARC galaxies.As illustrated in Figure 4 of Paper I and discussed further in Appendix B, only RCs that reach very low accelerations can probe the EFE.We select 143 galaxies as described in Appendix B. The fitted ẽ values of these galaxies are compared with the cosmic large-scale structure.
Considering the equivalence of Equation (2) and Equation (A2) with e/ √ 1 + e = ẽ, we expect the fitted value of ẽ from this work to agree with e/ √ 1 + e from Paper I. Figure 1 shows that there is a good match between the two for all galaxies except some at large ẽ.
We stress that our work is not specific to a particular MOND theory because Equation ( 2) is valid in both AQUAL and QUMOND for the idealized onedimensional case.In Appendix A, we show that Equa- tion ( 2) is a good approximation to numerical calculations in the AQUAL context.More generally, Equation (2) can be viewed simply as a fitting function to characterise non-flat outer RCs; then the extra free parameter e N is expected to be equal to the external Newtonian gravitational field (in units of a 0 ) only if the MOND paradigm is essentially correct.

THE NEWTONIAN FIELD FROM BARYONIC LARGE SCALE STRUCTURE
We estimate the environmental Newtonian field g Ne,env by summing the contributions from baryonic mass in various forms in the nearby Universe.In MOND, gravitational collapse is expected to become strongly non-linear at high z (e.g.McGaugh 2015), so the matter power spectrum at z 0 cannot be estimated analytically; one would need large-scale hydrodynamical simulations in a fully-fledged MOND cosmology.We therefore compile a census of observed baryonic structures to estimate their Newtonian accelerations directly.We use four main catalogues: the 2M++ galaxy catalogue (Lavaux & Hudson 2011), the MCXC galaxy clusters catalogue (Piffaretti et al. 2011), the Karachentsev galaxy catalogue (Karachentsev et al. 2018;Karachentsev & Kaisina 2019) and the NASA-Sloan Atlas (NSA) 1 of galaxies.2M++ and MCXC are used for the approximate all-sky calculation of g Ne,env (Section 3.1) while MCXC, Karachentsev and the NSA are used for the more detailed calculation within the SDSS footprint (Section 3.2).

An all-sky estimate of g Ne,env
We begin with an approximate all-sky calculation of g Ne,env to give a sense of its magnitude and investigate its angular variation.This is based primarily on 2M++ (Lavaux & Hudson 2011), an all-sky catalogue comprising the Two-Micron-All-Sky-Survey (2MASS) Redshift Survey (2MRS), the Six-Degree-Field Galaxy Redshift Survey Data Release 3 (6dFGRS-DR3), and the Sloan Digital Sky Survey Data Release 7 (SDSS-DR7).The catalogue reaches a depth of K S ≤ 11.5 over the full sky, and K S ≤ 12.5 in the 6dF and SDSS regions.The Galactic plane Zone of Avoidance is populated with mock galaxies by cloning large scale structure just outside the Zone to reproduce both 1-and 2-point statistics of the density field.In view of the approximate nature of this preliminary calculation, we estimate distances from the galaxy's recession velocity in the rest frame of the Cosmic Microwave Background assuming a Hubble constant of H 0 = 73 km s −1 Mpc −1 (as assumed in SPARC, Lelli, McGaugh & Schombert 2016).We estimate stellar masses assuming a mass-to-light ratio of 0.64 in the K S band (consistent with 0.5 in the Spitzer [3.6] band as assumed in SPARC; McGaugh & Schombert 2014).We follow the prescription of Lavaux & Hudson (2011) for relating apparent and absolute magnitude.We remove sources with K S > 11.5 which are preferentially seen in the deeper 6dF and SDSS fields and cut the catalogue at 200 Mpc.This leaves 54,483 galaxies out of 69,160 in total.

Adding the intra-cluster medium
We also consider the hot intra-cluster medium of galaxy clusters.MCXC is a meta-catalogue of 1743 Xray detected clusters from 12 separate catalogues, many of which derive from the ROSAT satellite (Voges et al. 1999).The catalogue provides redshift, coordinates, Xray luminosity L 500 in the range 0.1-2.4keV standardised between component catalogues at a critical overdensity of 500, and the total mass M 500 .This ΛCDM-based mass is not suitable for our purposes here.In MOND, the mass that sources g Ne,env is simply the total baryonic mass.The observed baryonic masses of galaxy clusters, however, systematically disagree with the MONDpredicted total masses (from X-ray hydrostatic equilibrium) by a factor ∼2 (Sanders 1999).Recent X-ray data show a smaller discrepancy (Ettori et al. 2019) but it remains clear that MOND has a missing mass problem in cluster centres.This discrepancy may be driven by undetected baryons (such as very cold gas, e.g.Milgrom 2008), massive sterile neutrinos (Angus et al. 2008) or a gravitating scalar field (Skordis & Zlośnik 2020).In the following, we will estimate the MOND dynamical mass and use it to calculate g Ne,env regardless of the nature of the missing cluster mass.
We begin with a scaling relation between L 500 and the observed mass of X-ray-emitting gas (M X ) from an X-ray flux-limited sample of 62 clusters (HIFLUGCS; Zhang et al. 2011): where A = 44.44,B = 1.11 and E(z) = 0.3(1+z) 3 +0.7.
We then build a scaling relation between the observed gas mass and the MOND dynamical mass (M MOND ) using the results of the hydrostatic equilibrium analysis of Angus et al. (2008) for 26 X-ray emitting clusters, finding C = 3.814 and D = 0.728 to provide a good fit to their data.
The results for the clusters' masses and hence g Ne,env are not significantly altered if we use the L 500 −M X relation of Mantz et al. (2016) instead of Equation ( 3), or if we use the results of Angus et al. (2008) to convert M 500 from MCXC directly to M MOND .The latter approach has two disadvantages: first it is more sensitive to the ΛCDM-driven choice to measure the mass within a critical overdensity of 500, and second it involves an implicit L X −M 500 scaling relation assumed by the MCXC team.M MOND is a factor of several lower than M 500 , so it is important not to use M 500 directly.

Total gravitational field
We calculate g Ne,env by linearly summing the g Ne values of each of the objects in the combined catalogue.For galactic sources -and clusters beyond R 500 -we treat the sources as point objects because the distances between the source and test points typically greatly exceed the dimensions of the source objects themselves.In case a test point is a distance d < R 500 from the centre of a cluster, we scale the g Ne,env contribution of the cluster by a factor (d/R 500 ) 3 , i.e. approximating the cluster density as uniform.This does not appreciably affect the results.As our primary interest here is in the angular variation of g Ne,env , we calculate it on a healpix (Górski et al. 2005) grid of nside=54 (∆RA, ∆Dec 1 deg) every 1.5 Mpc in line-of-sight distance out to 150 Mpc.To remove numerical artifacts we exclude source objects within 10 kpc of a test point when evaluating g Ne,env at that point.
As we discuss in detail below, many sources are neglected in this method so that g Ne,env is underestimated on average.This underestimation is not however expected to be a strong function of sky position, so the calculation should be sufficient to estimate the relative variation of g Ne,env with {RA, Dec}.

A more precise estimate within the SDSS footprint
We can calculate g Ne,env in greater detail within the SDSS footprint, where deeper photometry is available.
For this purpose we use the NSA, a value added catalogue based primarily on SDSS data.We calculate g Ne,env within two regions defined by {128 < RA/deg < 230, 0 < Dec/deg < 60} and {−25 < RA/deg < 25, −9 < Dec/deg < 30}.We use elpetro mass to estimate stellar mass and zdist with H 0 = 73 km s −1 Mpc −1 to estimate distances.We only cut sources with D > 500 Mpc, to ensure that we include sufficient long-range contributions to g Ne,env .The intracluster medium is added in the same way as in Section 3.1.1.In addition, we develop a probabilistic Monte Carlo framework to propagate uncertainties in various input quantities into g Ne,env , as described in more detail below.
In the Local Volume (D < 11 Mpc), we replace the NSA with the catalogue of Karachentsev et al. (2018) which is complete above a B-band absolute magnitude of about −12 mag and therefore contains fainter galaxies across the sky.It also provides both K-band luminosities and HI masses.We estimate stellar masses from K-band luminosities, again adopting M/L = 0.64, and cold gas masses from HI masses using the molecular gas correction of McGaugh et al. (2020): (5) If only an upper limit on M HI is recorded in the catalogue, then we estimate it from a scaling relation with M (see Section 3.2.2),with a cap at 0.8 M HI,max (this choice does not affect g Ne,env appreciably).The final Karachentsev catalogue contains 1109 galaxies to a sufficiently low mass that a correction for galaxies below the completeness limit is unnecessary.

Adding galaxies below the SDSS detection limit
The SDSS is limited to m r < 17.77 (m r < 17.6 conservatively), so galaxies of increasing stellar mass are progressively absent from the NSA at larger distances.If not corrected for, this would cause the average value of g Ne,env to fall with distance.To mitigate this selection bias, we add mock galaxies too faint to be measured by the survey.We use the triple-Schechter fit to the Li & White (2009) stellar mass function (SMF)2 to model stellar masses for mock galaxies.In particular, we take 80 logarithmically uniform bins in stellar mass M between 10 5 and 10 13 M , and 50 linear bins in distance D between 11 and 500 Mpc.In each {M , D} bin we use the SMF to calculate the expected number of galaxies.We then assign absolute r-band magnitudes to these mock galaxies using the scaling relation with E = −1.991,F = −0.265and a Gaussian scatter in M r of 0.385 dex, fitted directly from the NSA data, and convert to apparent magnitude in the same manner as for 2M++.We conservatively consider the galaxies visible to SDSS to be those with m r < 17.6.This lets us calculate the average fraction of the total stellar mass that is included in each radial bin when using the NSA catalogue directly; the remaining mass must now be added back in so that g Ne,env is not biased low at large D. Ideally a similar completeness correction would be applied to the galaxy cluster catalogue.We do not attempt this here because the cluster mass function is not well characterised, but note that the higher intrinsic brightness of galaxy clusters makes them visible to much larger distances than individual galaxies, so we expect the problem to be less severe.
While the above calculates the 1-point function of the missing mass, its 2-point function (i.e.clustering) is also relevant to the g Ne,env that it produces.To model this we consider the two most extreme scenarios: we assume either that the mock galaxies are randomly distributed in space within each distance annulus (i.e.unclustered), or that they are satellites of the galaxies included in the NSA so that they may be considered coincident with them (i.e.maximum clustering).These models bracket the full range of possible clustering strengths; we will refer to them as the "no clustering" and "max clustering" models respectively.
For the "max clustering" method we multiply the mass of each NSA galaxy by the reciprocal of the mass fraction calculated above for the corresponding distance bin, so that the correct total amount of mass is modelled in each bin.For the "no clustering" method we explicitly add mock galaxies with m r > 17.6 to our source catalogue before calculating g Ne,env ; this is done separately for each mass and distance bin described above, and the galaxies are given random positions within the bin.This primarily increases the uncertainty in g Ne,env when the random positions of the mock galaxies are marginalised over by Monte Carlo sampling, while the max clustering method primarily boosts g Ne,env .We find this boost to g Ne,env to be ∼0.1 dex, and hence not negligible.On the other hand, altering the widths of the mass or distance bins, the functional form of the SMF (within reasonable limits), or the minimum stellar mass down to which mock galaxies are included has a negligible effect on the results.Note that these methods for adding undetected galaxies preserve any over-or under-densities in galaxy number counts within the NSA as a function of distance, allowing them to manifest as a distance or other environment-dependence of g Ne,env .

Adding gas to galaxies
Galaxies contain appreciable amount of gas, either "cold" (T 10 4 K) gas, residing in star-forming atomic or molecular disks, or "hot" (T 10 5 K) gas, residing in X-ray emitting ionised halos.The gas content depends on both stellar mass and morphology, with lowmass late-type galaxies having a larger fraction of cold gas and high-mass early-type galaxies having a larger fraction of hot gas.
To model the gas content of our source galaxies, we begin by estimating the galaxy type as a function of its stellar mass.From Henriques et al. (2015), the earlytype fraction is given by f early 0.28 log(M /M ) − 2.12.We use this as a probability to assign an earlytype or late-type flag to each source galaxy, separately within each Monte Carlo realisation of the model so that the uncertainty it induces is naturally propagated into g Ne,env .For late-type galaxies, the gas mass is given by where the first and second terms in the sum consider, respectively, atomic (Lelli, McGaugh & Schombert 2016) and molecular (McGaugh et al. 2020) gas.For earlytype galaxies, the gas mass is instead given by We derived this scaling relation considering 94 earlytype galaxies with hot gas masses from the Chandra Xray observatory (Babyk et al. 2018) and stellar masses from our own WISE photometry (following the same procedures as in Lelli, McGaugh & Schombert 2016), adopting M /L W 1 = 0.7 as appropriate for early-type galaxies (Schombert et al. 2019).We take a scatter in both log(M g,cold /M ) and log(M g,hot /M ) of 0.2 dex, and add gas masses in the same way to the faint mock galaxies of Section 3.2.1.

Adding baryonic mass outside the SDSS footprint
Thus far, our catalogue of g Ne,env sources is all-sky within 11 Mpc but restricted to the SDSS footprint beyond that.Since gravity is a long-range force, this will lead to an underestimation of g Ne,env on average at D > 11 Mpc and a spurious dependence of g Ne,env and its direction on sky position within the NSA footprint.To mitigate this effect, we add a uniform grid of mass outside the SDSS footprint out to 500 Mpc.We calculate the density of this grid by integrating the SMF over all masses, including the corrections for cold and hot gas as a function of stellar mass as described above, as well as adding in the average density of cluster mass.This gives ρ tot = 10 8.549 M Mpc −3 .We use a Cartesian grid of spacing 5 Mpc in each direction, so that each grid cell has a mass of 4.4 × 10 10 M .To soften its g Ne contribution, we assume this mass to have a spherical Gaussian density profile with width σ = 1.25 Mpc, so that the gravitational field it sources at a displacement r from its centre is (9) This can have a ∼0.2 dex effect on g Ne,env and tends to reduce variation in g Ne,env with distance due to the homogeneity of the grid mass.The result is not significantly altered if the grid spacing is reduced to 2.5 Mpc, the SMF of Li & White (2009) is replaced by that of Moffett et al. (2016) or Bernardi et al. (2013), or if the grid cells are uniformly spaced in D 2 rather than D so that finer structure is resolved at smaller D. This indicates that the manner in which mass is distributed outside the SDSS footprint is not critical for the gravitational field within it, although it is important that the correct mass on average is modelled.

Adding cosmological missing baryons
The above method accounts for baryons in each of the forms in which they are readily visible: stars, cold and hot gas in galaxies, and an X-ray-emitting intra-cluster medium.However, there may be "missing baryons" not found in one of these forms.This is a necessity in ΛCDM where Ω b is known precisely from the Cosmic Microwave Background, but found to be much larger than that implied by a census of visible mass locally (e.g.Shull et al. 2012).The result for Ω b is likely to hold in a MOND cosmology as well, where at early times (e.g. during Big Bang Nucleosynthesis) the high average density of the Universe makes the MOND modification to gravity ineffective.This is however subject to the considerable uncertainty in MOND cosmology, stemming largely from the unknown relativistic parent theory; see, e.g., Famaey & McGaugh (2012).In the case of the relativistic MOND theory of Skordis & Zlośnik (2020), we expect Big Bang Nucleosynthesis to work in MOND in nearly the same way as in ΛCDM, so the missing baryons problem should be analogous in both cosmologies.
We adopt the census of baryons from Shull et al. (2012), who find that baryons in the forms we explicitly model amount to ∼ Ω b /8 (with Ω b 0.046; Macquart et al. 2020), meaning that there are potentially ∼8× more baryons than we have included so far.As an additional consistency check on our mass model, we calculate that it produces total densities of all baryonic component within the SDSS footprint that are approximately consistent with the results of Shull et al. (2012).Similarly to our method for adding faint galaxies, we consider two extreme clustering models for these missing baryons: 1) they are completely unclustered (or equivalently non-existent) so do not alter g Ne,env , or 2) they are maximally clustered with the structures we model and equally distributed between them, causing a uniform increase in the magnitude of the g Ne,env field by a factor of 8.
Our assumption is that the cosmological missing baryons are too far from the centres of galaxies and clusters to appreciably alter their measured kinematics.This is reasonable even for the "max clustering" model because the scales involved in calculating g Ne,env ( 10 Mpc) are much larger than those over which kinematics can be probed (<Mpc).These missing baryons are therefore fully independent of those potentially in-voked in Sec.3.1.1 to explain the dynamical masses of clusters.Nevertheless it is possible that the potential missing baryons inside clusters may account for some fraction of the cosmological missing baryons, leading us to overestimate g Ne,env when multiplying by a factor of 8.This effect is however small (∼ 10%), and the max clustering model is intended to give an upper bound on plausible g Ne,env values, so that it brackets the possible range along with the "no clustering" model.In fact, the difference between M X and M MOND in clusters may be due to neutrinos or a massive scalar field rather than dark baryons.

Modelling uncertainties
Finally, it is important to estimate the uncertainty in g Ne,env due to uncertainties in our baryonic mass model of the Universe.This includes uncertainties in the scaling relations that we use to determine gas and cluster masses, the locations of mock galaxies above the magnitude limit of SDSS, and the basic properties of the source objects themselves (e.g.distances and stellar masses).Our Monte Carlo framework makes it easy to propagate these uncertainties: at each test point we evaluate g Ne,env 2000 times, in each case randomly sampling from the distributions describing the input variables.This is a sufficient number for convergence of the g Ne,env PDFs.
We use a 30% uncertainty in M and M MOND , describing uncertainties in stellar mass-to-light ratios and the calculation of the MOND dynamical mass of galaxy clusters.We use a 0.2 dex scatter in cold and hot gas masses of galaxies, and 0.4 dex in cluster masses to account for the additional uncertainty in the conversion from X-ray luminosity to hot gas mass.Following Lelli, McGaugh & Schombert (2016), distance uncertainties are based on the following scheme (suitable for Hubble-flow distances): 30% for D < 20 Mpc, 25% for 20 ≤ D/Mpc < 40, 20% for 40 ≤ D/Mpc < 60, 15% for 60 ≤ D/Mpc < 80, and 10% for D > 80 Mpc.
The remaining uncertainties derive from the assignment of an early-type/late-type flag to each galaxy, and the distribution of faint mock galaxies within each distance bin in the "no clustering" model for filling in faint galaxies.We note that the "no clustering" model for cosmic missing baryons is different to that for faint mock galaxies in that it does not discretize the missing mass and assign it a random position.The uncertainty induced by cosmic missing baryons is therefore entirely in the systematic difference between the no clustering and max clustering models.We summarise the resulting g Ne,env distribution over all Monte Carlo realisations at each test point by its mean and standard deviation.
We evaluate g Ne,env at the positions of the Karachentsev and NSA galaxies themselves out to 150 Mpc to assess its trend with environment and distance, and at the positions of the SPARC galaxies for a direct comparison with g Ne,fit .We exclude from the source catalogue the test galaxy under consideration, either by simply mask-ing out the corresponding galaxy from the source array (for a Karachentsev or NSA galaxy), or, for a SPARC galaxy, by using NED to determine if it is in the NSA and an RA/Dec match with tolerance 0.1 deg to determine if it is in the Karachentsev catalogue.We consider a cylindrical space of radius R = 150 Mpc and height Z = 50 Mpc in the northern sky (declination δ > 0 • ), where the vast majority of fitted galaxies (126) are located.In Figure 2, the SPARC galaxies are projected on the equatorial plane of this cylindrical space together with galaxies from the 2M++ catalog (Lavaux & Hudson 2011).At 22 R/Mpc 45, there is an under-dense region in which SPARC galaxies tend to give lower ẽ values than the median.At 45 R/Mpc 100, instead, there are two massive structures: the CfA2 great wall (Geller & Huchra 1989) towards Right Ascension α 180 • and the Perseus-Pisces supercluster (Joeveer & Einasto 1978) towards α 45 • .Several SPARC galaxies are located in these over-dense structures and tend to give higher ẽ value than the median.This seems qualitatively consistent with an EFE interpretation of the ẽ values from RC fits.The situation is more complex at R < 22 Mpc due to a finger-like structure towards α = 190 • , which varies with height and disappears at Z > 11 Mpc.The precise geometry of this structure may affect the ẽ values in a way that cannot be simply inferred by looking at projected maps.However, note that the two "golden galaxies" NGC5033 and NGC5055 highlighted in Paper I are near the Virgo supercluster, while the control galaxy NGC6674 is in the underdense ring.The other control galaxy NGC1090 is in the southern sky, so it is not included in Figure 2.
Figure 3 shows the distribution of ẽ versus cylindrical radius R for the 126 SPARC galaxies within the cylindrical space.It is evident that galaxies in the underdense region at 22 Mpc < R < 45 Mpc have low values of ẽ consistent with zero environmental fields.The median value within this radial bin is 0.016 +0.015 −0.027 (with bootstrap uncertainties), thus consistent with no-EFE detection as expected from the large-scale distribution of galaxies.On the other hand, galaxies in a radial bin at R > 45 Mpc, in which the CfA2 great wall and the Persus-Pisces supercluster are located, display a median value of ẽ = 0.103 +0.030 −0.008 , significantly different from zero at ∼ 13σ.This shows a striking contrast with the middle bin.The majority of the galaxies in the third bin, indeed, are associated with the CfA2 great wall or the Perseus-Pisces supercluster.
These results demonstrate qualitative consistency between the Newtonian external field inferred from galactic RCs and the large-scale galaxy distribution.To make a more detailed comparison, Figure 4 shows a full-sky Mollweide projection of the mean environmental field     g Ne,env from the 2M++ plus MCXC calculation (Section 3.1) over different distance bins: 0 < D/Mpc < 11, 11 < D/Mpc < 30, 30 < D/Mpc < 60 and 60 < D/Mpc < 120.These are chosen to represent (1) the Local Volume, (2) the remainder of the overdense region we identify in the SDSS footprint (see Section 4.2), (3) the underdense region further out, and (4) the CfA2 great wall and the Perseus-Pisces supercluster regions.Note that here we use three-dimensional distance D from the Milky Way, rather than cylindrical radius R, to provide an alternative perspective.Again we find qualitative agreement between the RC-fitted field g Ne,fit and the environmental field g Ne,env in terms of their spatial variations.In particular, the Local Void (Tully et al. 2019) is clearly visible as the blue regions with low values of g Ne,env at D < 11 Mpc; galaxies within the void tend to have values of g Ne,fit below the average, as expected.

Quantitative comparison in the SDSS footprint
We now consider a quantitative comparison of ẽ with the environmental field estimated in Section 3.2.To begin, we calculate g Ne,env at the positions of the Karachentsev and NSA galaxies themselves out to 150 Mpc.This provides a fair tracing of the gravitational field over the nearby large scale structure.Figure 5 plots the mean g Ne,env values as a function of distance.We use here the "max clustering" model, which we will find below to match the fitted values better than the "no clustering" case, although the qualitative trends are the same.We see that g Ne,env falls from D = 0 to D 50 Mpc, then rises again out to ∼100 Mpc before levelling off.For both D < 40 Mpc and 80 < D/Mpc < 120, g Ne,env is significantly larger towards α 180 • than elsewhere, indicating excess structure in that direction.This is largely due to the Virgo cluster at small distance and the CfA2 great wall at D 100 Mpc.The part of the SDSS footprint in the Southern Galactic Cap (−25 • < α < 25 • , −9 • < δ < 30 • ) has a weaker field on the whole, especially at D 20 Mpc where it can be seen as a disjoint band of low-g Ne,env points.Very similar trends are seen if the distance D is replaced by the cylindrical distance R adopted in Figures 2 and 3.
We now evaluate g Ne,env at the positions of the SPARC galaxies.Restricting to the SDSS footprint retains 109 SPARC galaxies, although only 90 of these have a measured ẽ satisfying the cut x 0,3 < −10.6.Figure 6 shows log(e N,env ) for the SPARC galaxies as a function of distance, including their full uncertainties, for both the no clustering and max clustering models.Although the statistical uncertainties are considerable, the distance trend is still clearly visible, as is the factor ∼8 difference between the two clustering models driven by the influence of baryons not found in stars, cold gas, or the ICM.
Figure 7 exhibits individual values of ẽfit against ẽenv for the sample of 90 SPARC galaxies as a function of distance.Remarkably, the median value of ẽfit from rotation-curve fits is well within the allowed values from the two clustering models.This is highly nontrivial because, in a general DM context, ẽ is merely a fitting parameter to describe the outer shapes of rotation curves.It could have taken a value across orders of magnitude, so there is no a priori reason why its median value should lie within the boundaries of the ẽenv calculation from the large-scale distribution of baryonic matter.Moreover, both ẽfit and ẽenv are lower than average in the underdense region 30 < D/Mpc < 60 while clearly higher than average in the overdense region D > 60Mpc, consistent with the trend of ẽfit for all 126 galaxies within Z < 50 Mpc in the northern sky shown in Figure 3.
To check the agreement between ẽfit and ẽenv in more detail, Figure 8 shows the distribution of their differences.Remarkably, there is an excellent statistical match between ẽfit and ẽenv of the max clustering case while there is some tension with the no clustering case.This result suggests that accounting for missing baryons is important in MOND.Moreover, it agrees with previous indications that the missing baryons are likely significantly correlated with galaxies and groups (e.g.Nicastro et al. 2018;Lim et al. 2020;Das et al. 2019Das et al. , 2020)).
Finally, we check directly the correlation between ẽenv and ẽfit for the sample of 90 galaxies shown in Figures 7 and 8 using Pearson's linear correlation coefficient r.Because ẽfit has a large uncertainty ranging from 0.004 to 0.194 with a median value of 0.047 for these galaxies, we expect the correlation to be weak.In Figure 9 we compare the measured correlation coefficient with the expected distribution from a Monte Carlo simulation.This simulation assumes that ẽfit =ẽ env fundamentally, with ẽfit scattered by the error budgets of both quantities.For ẽenv we use the mean of the max clustering and no clustering values.For the uncertainty of ẽenv , we take the average of the statistical errors of the max and no clustering cases for the statistical error and one half of the difference between the max and no clustering cases for the systematic error.For the uncertainty of ẽfit , we take the measured values (see Appendix B) for the statistical error and 0.02 for the systematic error of our model (see Sec. A.2).These errors are added in quadrature for each galaxy, with the result used to scatter ẽfit in each Monte Carlo realization.
We find the measured coefficient of r = 0.057 to agree well with the distribution from the mock data, with a 27% probability of a lower r under the hypothesis ẽfit =ẽ env .The mock data have a RMS scatter of 0.077 ± 0.008 for ẽfit .This is about 3/5 of the measured scatter for this sample, indicating that the uncertainties of ẽfit and/or ẽenv may be somewhat underestimated.Regardless of the precise error model, it is clear that the theoretically expected correlation between ẽfit and ẽenv values of individual galaxies cannot be directly inferred given the current observational uncertainties.The existing data, however, are consistent with the existence of an intrinsic correlation between ẽfit and ẽenv .2) we investigate the EFE in the RCs of 143 SPARC galaxies, using a dimensionless parameter ẽ related to the external Newtonian gravitational field.We find that ẽ is well correlated with the observed large scale structure of baryonic mass distribution (Figure 2).The fitted values are consistent with zero in an underdense ring region at cylindrical radii 22 Mpc < R < 45 Mpc, while they are a factor of two higher than the average external field at or near the CfA2 great wall and the Perseus-Pisces supercluster (Figure 3).In this high density region the EFE is statistically detected at ∼13σ.The correlation between ẽfit and large-scale structure extends qualitatively across the whole sky (Figure 4).
Next, we compare the RC-fitted values of ẽ with fully independent estimates of the Newtonian fields from the baryonic mass distribution in the nearby Universe (using the 2M++, NSA, Karachentsev and MCXC catalogues).We include corrections for undetected galaxies below the magnitude limits of the surveys, for cold and hot gas in galaxies, for the intracluster medium, and for the possible presence of "missing baryons" in the warmhot intergalactic medium, with full propagation of uncertainties in a Monte Carlo framework.We find good agreement between ẽfit with ẽenv values, both in their distance dependence (Figure 7) and the average value of the Newtonian field (Figure 8).In particular, ẽfit is well matched with ẽenv when the missing baryons needed to bring the local baryon density into agreement with Big  Bang Nucleosynthesis (Ω b h 2 = 0.022, Cooke et al. 2018) are strongly clustered with galaxies and galaxy clusters.
The qualitative correlation between ẽfit and ẽenv is not very sensitive to the functional form assumed for the EFE model (Equation 2) because the relative values of inferred field strengths are similar for different models.However, the absolute value of ẽfit may vary, leading to systematically different mean and median values.For example, we fitted RCs using an alternative function based on numerical simulations (Haghi et al. 2019) and found a systematic shift of ∼0.02 in ẽ on average (Figure 10).Until the fitted values of ẽ of Equation ( 2) are thoroughly tested through realistic numerical simulations of disk galaxies, systematic shifts of order ∼ 0.02 are possible.This could have implications for the abundance and clustering of cosmic baryons.Our current values of ẽfit are within the limits set by the "max clustering" and "no clustering" models of missing baryons, although the max clustering case is preferred (Figures 7  and 8).
It may be possible for g Ne,env to exceed the prediction of the max clustering model due to the gravitating scalar fields present in some relativistic MOND theories, such as that of Skordis & Zlośnik (2020).Indeed, the EFE may manifest itself somewhat differently in this theory than in AQUAL and QUMOND because the free function in the Lagrangian (which plays the role of the interpolating function in AQUAL) has an oscillating piece related to the time-variation of the scalar field as well as a more traditional MOND term.It will therefore be important to determine the exact behaviour of the EFE in this and similar theories.
The agreement between the external fields from RC fits and large-scale structure of baryonic matter that we have discovered is predicted by MOND modified gravity theories, supplementing the evidence for the EFE (and a fortiori a violation of the SEP) presented in Paper I. This correlation is nontrivial from the ΛCDM point of view.It is not naturally expected that the density of the cosmic web would affect the dynamics of the disk in the inner region of the DM halo, never mind in a manner mimicking the MOND EFE.We are not aware of any ΛCDM simulations predicting that RCs in stronger-field regions would tend to decline while those in low density regions do not.This adds to the "small-scale problems" with ΛCDM (Bullock & Boylan-Kolchin 2017;Kroupa 2012).In this context, we note that challenges to ΛCDM cosmology have also been emphasized recently on the scales of galaxy clusters and voids (Haslbauer et al. 2020;Asencio et al. 2021).
The external Newtonian fields encapsulated in ẽ will provide a benchmark for more detailed numerical studies of the MOND EFE, taking full account of the relative orientation of the external field and the threedimensional morphology and kinematics of galactic disks.Various MOND theories may have different EFE predictions for RCs of disk galaxies under the same g Ne,env .Although our fitted values of ẽ are modeldependent, they are empirical in nature and may distinguish existing theories when accurate numerical results become available.Here we used a fitting function that approximately agrees with existing AQUAL simulations (see Appendix A), so the agreement between g Ne,fit and g Ne,env may indicate that AQUAL is close to the correct nonrelativistic theory.
Besides distinguishing between MOND formulations, our results may be used to infer the properties of "missing baryons" in a MOND context.Our preliminary result in this regard is that baryons in the warm-hot intergalactic medium needed for Ω b 0.046 (Macquart et al. 2020) are likely strongly clustered with galaxies and clusters.
In summary, the agreement between our fit results and the observed baryonic large-scale structure matches well the MOND prediction and reinforces the detection of the EFE in Paper I, pointing to a breakdown of the SEP.These results support the modified gravity hypothesis and reveal the possibility of using the internal dynamics of galaxies to study the large-scale distribution of cosmic baryons.

Chae et al. APPENDIX
A. A SIMPLE ONE-DIMENSIONAL FORMALISM FOR THE EFE In the MOND context, it is hard to quantify generally the effect of a constant external gravitational field on a galaxy disk because, for a test particle along a given orbit, the angle between the internal acceleration and the external gravitational field will vary continually, leading to a time-varying EFE.Broadly speaking, we expect that the outermost orbits will become slightly non-circular and possibly tilted with respect to the inner galaxy plane.Asymmetric, lopsided, and/or warped H I disks are frequently observed in nearby galaxies (e.g., Sancisi et al. 2008) but their relation to the EFE can only be studied with detailed numerical simulations.For example, Brada & Milgrom (2000) show that the asymmetric warp of the Milky Way may be induced by the EFE from the Magellanic Clouds, while Banik et al. (2020) reproduce the dynamical properties of the Local-Group spiral galaxy M33 as the result of the EFE from Andromeda.The main EFE, however, is a decrease in the circular velocity in the MOND regime.
In this work we are particularly interested in how the EFE depends on the external Newtonian field.Here we describe a simple one-dimensional formalism to model the EFE-induced outer decline of RCs, neglecting higher order effects due to the full three-dimensional nature of the problem, such as the development of non-circular motions, warps, and asymmetries.Incidentally, this formalism turns out to be equivalent in AQUAL and QUMOND with a heuristic MOND relation for the constant external field.
A.1.The EFE in AQUAL and QUMOND If a system is freely falling under a constant MOND external field g e and we approximate the EFE by a onedimensional AQUAL equation (see Section 6.3 of Famaey & McGaugh 2012), the gravitational acceleration g of a test particle within the system can be determined by where A e = (1 + e/2)/(1 + e), which is redefined here so that its limit tends to 1, and B e = 1 + e.For typical values e 0.1 (Paper I) we have A e ≈ B e ≈ 1.It is then apparent that the EFE-dominating last term scales approximately linearly with e.For e → 0 (i.e., truly isolated galaxies), one recovers ν 0 (y) = 1/2 + 1/4 + 1/y which is the inverse function of the simple IF.
If the external Newtonian field of the above system is g Ne and we use the QUMOND formalism, g is given by where ν(y) = ν 0 (y).Then, the ratio g/g N is given by Equation ( 2).Although it is not obvious algebraically, one can show that Equation ( 2) is equivalent to Equation (A2) for positive (physical) external field with the transformation e 2 /(1 + e) = e N .This means that although Equation ( 2) is derived directly from Equation (A3), it is also consistent with Equation (A1).
A.2.Comparison with the analytic point-mass weak-field limit and numerical simulations Our fitting function (Equation ( 2)) is a heuristic first-order model of EFE that happens to be equivalent in AQUAL and QUMOND with the correspondence ẽ = e/ √ 1 + e.For general dynamical systems, for the same Newtonian external field, these two Lagrangian theories predict different values of the ratio g/g N , which can be calculated only through numerical methods.Here we compare Equation (2) with analytic solutions of a point mass in the deep-MOND limit as well as previously published numerical results.
Far away from any finite mass distribution, the point-mass approximation will be valid regardless of the mass distribution.For the QUMOND EFE limit, we start from Equation (60) by Milgrom (2010) that gives the approximate  2019) (blue dashed curve).These are then compared with the analytic point-mass weak-field limit for different values of e, considering an average angle of 60 • between the external field and the axis of a circular orbit.In the acceleration range probed by the SPARC RCs, the difference between our model and the numerical result is relatively small in terms of e (or ẽ).See the text for further details.
analytic internal potential ψ for a point mass M under an external Newtonian field of strength g Ne = e N a 0 .Taking the gradient of ψ and considering an angle θ between g Ne and the rotation axis of a test particle on a quasi-circular orbit, the azimuthal average of the radial acceleration g = |r • ∇ψ| is given by where ν 0 (y) = 1/2 + 1/4 + 1/y and ν0 (y) ≡ d ln ν 0 (y)/d ln y which tends to −1/2 at large r.Thus, at large r the two extreme cases θ = 0 • and θ = 90 • show a difference of just 1/6 of the value.In other words, the direction of the external field has a minor effect.Similarly, for the AQUAL limit, we start from Equation (66) of Milgrom (2010) considering the MOND external field g e = e a 0 and obtain where µ 0 (x) = x/(1 + x) and L 0 (x) ≡ d ln µ 0 (x)/d ln x.
As for numerical results on the ratio g/g N , we consider a functional relation derived by Haghi et al. (2019) from AQUAL-based numerical simulations of fully pressure-supported spherical systems.In a spherical shell of radius r of a pressure-supported system, the measured velocity dispersion is a collective representation of various random orbits at different angles with respect to the external field.In this sense the velocity dispersion may capture a direction-averaged EFE.Then, for the case of an isotropic velocity distribution, the radial acceleration is given by g(r) = 3σ 2 /r, where σ is the one-dimensional velocity dispersion.Thus, comparing σ with and without an external field, one can obtain an EFE-dependent g(r)/g N (r).Haghi et al. (2019) used the so-called "standard" IF (e.g., Famaey & McGaugh 2012) for their simulations, so we transform their result to that for the simple IF by multiplying the EFE-dependent g(r)/g N (r) by the ratio of the simple IF to the standard IF.
In Figure 10 we compare our EFE model (Equation (A2), or equivalently Equation (2)) with the numerical result by Haghi et al. (2019) and the weak-field analytic expectations of QUMOND (Equation (A4)) and AQUAL (Equation (A5)).We note that QUMOND predicts 7 percent higher g in the deep MOND limit.At log(g N /m s −2 ) ≈ −10.8 Chae et al. our model matches the numerical result, but at higher (lower) accelerations our g is higher (lower) than the numerical value.This means that for a declining RC covering the low acceleration range log(g N /m s −2 ) < −10.8, our model will give a lower value of ẽ than the numerical model.For log(g N /m s −2 ) > −10.8 the opposite may occur occasionally.In the weak-field limit, both our model and the numerical result give lower g than the point-mass expectation.Our model function with e = 0.040 matches the AQUAL point-mass limit with e = 0.064, exhibiting a difference of 0.024, whereas the numerical result exhibits a smaller difference of 0.004.Thus, the difference between our model and the numerical result for log(g N /m s −2 ) < −10.8 is 0.020 in terms of the AQUAL point-mass limit.A similar median difference of ∼ 0.02 in the fitted value of ẽ is found from fitting the observed RCs.This is smaller than the typical individual uncertainties of ∼ 0.04 -0.05 in ẽ.

B. FITTED VALUES OF THE PARAMETERS AND THE NEWTONIAN ENVIRONMENTAL FIELDS
Here we provide the fitted values of ẽ and galactic parameters for all 162 galaxies considered in this work.We also describe how we select a statistical sample of 143 galaxies whose RCs reach low enough accelerations for the fitted ẽ values to be meaningful.As illustrated in Figure 4 of Paper I only RCs that reach low enough accelerations can probe the EFE from the large-scale distribution of cosmic mass.In other words, if we used only RCs in a high-acceleration range, then the g N -g relation would have little sensitivity to ẽ (or e), and consequently the median of fitted ẽ (or e) for such RCs will be biased towards zero.We therefore select the RCs that can probe a low acceleration regime.We use the parameter x 0 introduced in Paper I, which is obtained by projecting a point (log g bar , log g obs ) (here g bar means g N ) to the curve corresponding to a flat RC (the red solid curve in Figure 10).An RC having low enough values of x 0 is expected to have sensitivity to ẽ (or e), so we require a RC to have at least 3 data points (log g bar , log g obs ) satisfying x 0 < x 0,cut .Here we do this by x 0,3 < x 0,cut where x 0,3 is the 3 rd lowest value of x 0 .For the threshold we choose x 0,cut = −10.6,so that selected RCs have at least three values of 10 x0 lower than 10 −10.6 m s −2 .Figure 11 shows that, at a higher threshold than this value, the median of ẽ starts to decrease indicating that the selected sample starts to be affected by about 20 RCs in a high acceleration range with little sensitivity to ẽ.
Because x 0 is always on a flat RC by construction, when we select RCs by x 0 , there is no selection bias to a more declining (or rising) RC.This can be seen clearly in Figure 11: as the threshold x 0,cut decreases from −10.6, the sample median of ẽ remains the same.However, if one selects RCs by requiring log g obs (R) to be lower than a certain threshold, then the selection is biased in favor of more declining RCs.This is because for a given g bar (R), lower values of g obs (R), thus declining RCs, are more likely to be selected by chance due to the unavoidable observational scatter in g obs (R).Selection based on log g bar (R) is biased in the opposite sense.Namely, less declining RCs are more likely to be selected by requiring a certain number of values of x to be lower than a certain threshold.This is because for a given g obs (R), a more declining Newtonian RC g bar (R) is more likely to be selected by chance.Our selection based on x 0 is robust against possible bias.Figure 12 shows the distribution of the fitted ẽ values for the selected sample of 143 galaxies as well as for the full sample.The median values for the samples are ẽ = 0.044 +0.012 −0.009 (full sample) and ẽ = 0.053 +0.008 −0.012 (selected sample), the latter of which corresponds to g Ne 0.0028a 0 .Here the uncertainty is estimated from the Monte Carlo distribution of the medians of many bootstrap resamples of the data (this method is simply referred to as "bootstrap" in this paper).This is a ∼4.4σ detection of positive ẽ.Also, from the histogram we have 97 cases of ẽ > 0 out of 143.The null hypothesis that ẽ is as likely to be <0 as >0 is ruled out at 4.1σ confidence based on binomial statistics.For the (contaminated) full sample the significance of ẽ > 0 is 3.5σ by binomial statistics.These results reinforce the results by Paper I based on Equation (A2) and a different sample selection.
Table 2 lists the RC-fitted values of the model parameters for all 162 SPARC galaxies modeled through MCMC simulations.Values of x 0,3 (the 3 rd lowest value of x 0 as defined above) are given for the purpose of selecting galaxies with x 0,3 < −10.6 for statistical analyses.The fitted value and its uncertainty are derived from the 50 percentile and the 15.9 and 84.1 percentiles of the posterior PDF respectively.
We assign four quality flags based on the shape of the ẽ posterior: P, A, B, C. In most cases the PDF is well peaked and the 50 percentile agrees well with the best-fit.Those galaxies are assigned PDF-quality 'P'.An example (NGC 5055) is exhibited in Figure 13 (top-left panel).For the remaining 24 galaxies with quality flags A, B, or C, the PDF of ẽ extends beyond the prior range −0.5 < ẽ < +0.5, and is not well peaked in some cases.For 9 galaxies with PDF-quality 'A', the PDF is well peaked within or a little outside the prior limits.An example (F563-V2) for this type is exhibited in Figure 13 (top-right panel).For 8 galaxies with PDF-quality 'B' the PDF has a peak with a long tail in the positive direction of ẽ.An example (NGC 2976) for this type is exhibited in Figure 13 (bottom-left panel).Finally, for 7 galaxies with PDF-quality 'C', the PDF is not well peaked but stays nearly flat as ẽ increases.An example (F579-V1) for this type is exhibited in Figure 13 (bottom-right panel).For the galaxies with PDF-quality of 'A', 'B', or 'C', the lower bound is still within the prior range.For those galaxies, the given value of ẽ in Table 2 may be considered as a lower bound within the prior range.
Table 3 lists the Newtonian environmental field strengths of 109 SPARC galaxies residing in the SDSS footprint for the "max clustering" and "no clustering" cases as described in Section 3.2.

P A B C
Figure 13.Examples showing the posterior PDFs of our model parameters for four types of PDF-quality described in the text of the Appendix.For types A, B, and C, an extended prior range for ẽ ≡ eN/ |eN| is used to illustrate the behavior of its PDF.
Thus, the values of ẽ for these cases are different from those given in Table 2.

Figure 1 .
Figure 1.Fitted ẽ values based on Equation (2) from this work are compared with e values based on Equation (A2) from Paper I for 152 overlapping galaxies.
Qualitative consistency with large-scale structure

Figure 2 .
Figure2.The top left panel shows a map of galaxies in a cylindrical space of the northern sky with a radius of 150 Mpc and a height of 50 Mpc.Galaxies are taken from the 2M++ catalog (accepting only those brighter than 11.5 in K-magnitude to prevent any direction bias) and the area of each point is proportional to the stellar mass of the galaxy.SPARC galaxies with external Newtonian fields from RC fits are indicated by coloured points.External Newtonian fields higher and lower than the median (ẽ = 0.053) are shown in red and blue respectively, and the opacity of the point is proportional to |ẽ−0.053|/(uncertainty of ẽ).The central R < 45Mpc region is sliced and zoomed-in in the other panels.These panels show an overdensity-underdensity contrast in the two opposite directions of RA α = 180 • and 0 • .Two "golden" galaxies and one control galaxy from Paper I are indicated.The first ring-like region with R = (22, 45) Mpc is relatively underdense and the fitted ẽ values are mostly lower in this region.The second ring-like region of R = (45, 100) Mpc includes the CfA2 great wall and the Perseus-Pisces supercluster.The fitted ẽ values are mostly higher in this region.

Figure 3 .
Figure3.RC-fitted values of ẽ versus the cylindrical radius defined in Figure2for 126 SPARC galaxies in the Northern hemisphere.The size of a dot is inversely proportional to the uncertainty on ẽ.The light-red band shows the median value and its uncertainty for all galaxies shown in Figure12.The red bars show median values (black line) and their uncertainties in three radial bins.The number of galaxies in each radial bin and their median values with bootstrap 1σ uncertainties are indicated at the top.ẽ is systematically lower in the underdense second bin and higher in the third bin.

Figure 4 .
Figure 4. All-sky distributions of the environmental field gNe,env from 2M++ galaxies and MCXC clusters in Mollweide projection and equatorial coordinates averaged across various distance ranges.The location of SPARC galaxies with independent estimates of gNe from RC-fits are shown: red and blue points indicate ẽfit ≥ 0.053 and ẽfit < 0.053, respectively; the opacity of the point increases with |ẽ fit − 0.053|/(uncertainty in ẽfit ) as in Figure 2.

Figure 5 .
Figure5.Variation of eN,envwith distance for the galaxies in the NSA and Karachentsev catalogs.Individual galaxies are colour-coded by Right Ascension (RA).The black lines show the mean trend (solid) and standard deviation (dashed) in bins of distance.This plot assumes the max clustering model for the missing baryons (cf.Figure6).

Figure 6 .
Figure6.Variation of eN,envwith distance for the SPARC galaxies within the NSA footprint.The "max clustering" model (blue) assumes that missing baryons are effectively coincident with observed structures, while the "no clustering" model (orange) distributes them uniformly in space.See Section 3.2.1 for details.

Figure 7 .
Figure 7.The external Newtonian gravitational fields from RC-fits (dots with errorbars), parametrized by ẽ, are compared with those from the large-scale distribution of baryonic matter (blue band) as a function of distance to 90 SPARC galaxies that probe the low acceleration regime and are within the SDSS footprint.The size of each black dot is inversely proportional to its statistical uncertainty.The median value of ẽfit and the associated uncertainty is represented by the light-salmon band.The medians and associated bootstrap uncertainties in 4 distance bins are shown by red bars, with the values indicated at the top.The blue band is obtained considering the "no clustering" (dashed line) and "max clustering" (solid line) models for the missing baryons (see Section 3.2.1),interpolating between the ẽenv values recorded at the location of each SPARC galaxy.The thin blue lines show the statistical 1σ uncertainty for each clustering model.The agreement of ẽenv with the median values of ẽfit is remarkable and not expected a priori outside a MOND context. 5. DISCUSSION AND CONCLUSION Using a simple, generic fitting function (Equation2) we investigate the EFE in the RCs of 143 SPARC galaxies, using a dimensionless parameter ẽ related to the external Newtonian gravitational field.We find that ẽ is well correlated with the observed large scale structure of baryonic mass distribution (Figure2).The fitted values are consistent with zero in an underdense ring region at cylindrical radii 22 Mpc < R < 45 Mpc, while they are a factor of two higher than the average external field at or near the CfA2 great wall and the Perseus-Pisces supercluster (Figure3).In this high density region the EFE is statistically detected at ∼13σ.The correlation between ẽfit and large-scale structure extends qualitatively across the whole sky (Figure4).

Figure 8 .
Figure 8. Statistical comparison of ẽfit and ẽenv for the subsample of 90 SPARC galaxies shown in Figure 7.The max clustering and no clustering cases (see Section 3) are considered.The excellent statistical match with the max clustering model indicates that "missing" cosmic baryons are strongly correlated with observed structures.

Figure 9 .
Figure 9. Observed value of Pearson's linear correlation coefficient r between ẽenv and ẽfit for the sample shown in Figures 7 and 8 (red line) and the expected probability distribution from many Monte Carlo realizations assuming ẽfit = ẽenv and scattering by the measured uncertainties of both quantities (histogram).For ẽenv we use the average of the max clustering and no clustering cases.
A1)where g N is the internal Newtonian field and µ(x) is the MOND interpolating function (IF).If we assume the simple IF µ(x) = x/(1 + x), then the ratio g/g N is given by Equation (6) of Paper I as ν e

Figure 10 .
Figure 10.The one-dimensional EFE model from Equation (A2) (red dashed curve), or equivalently Equation (2) with ẽ = e/ √ 1 + e, is compared with the numerical simulation result by Haghi et al. (2019) (blue dashed curve).These are then compared with the analytic point-mass weak-field limit for different values of e, considering an average angle of 60 • between the external field and the axis of a circular orbit.In the acceleration range probed by the SPARC RCs, the difference between our model and the numerical result is relatively small in terms of e (or ẽ).See the text for further details.

Figure 11 .
Figure 11.Dependence of the sample median of ẽ on the sample selection cut.Each sample is selected by x0,3 < x0,cut where x0,3 is the 3 rd lowest value of x0 as defined in Paper I (see the text for the details).The vertical dashed line indicates our selection.The errorbars indicate bootstrap uncertainties.As x0,cut increases, galaxies without enough low-acceleration RC points start to be included in the sample so that ẽ is biased to a lower value.

Figure 12 .
Figure 12.Distribution of the fitted values of ẽ.The left panel shows 162 SPARC galaxies with reliable RCs (Q < 3, see Lelli, McGaugh & Schombert 2016), while the right panel shows our selected sample of 143 galaxies with RCs sensitive to the EFE, as indicated in Figure 11.The errorbars on the histogram show Poissonian uncertainties Nj, where Nj is the number of galaxies in each bin.The red curve indicates a Gaussian fit whose rms width is ∼0.11.The uncertainty on the median or on the Gaussian mean is estimated through a bootstrap method.

Table 1 .
Summary of notations N,fit / |e N,fit | ẽenv √ eN,env Note.(1) The subscript 'fit' or 'env' may be dropped when it is obvious from the context.(2) The RC-fitted quantities allow negative values.

Table 2 .
Fitted model parameters