<figure><img src="image_1.png" /></figure>

HAL
open science

# Fifth force constraints from the separation of galaxy mass components

Harry Desmond, Pedro G. Ferreira, Guilhem Lavaux, Jens Jasche

► To cite this version:

Harry Desmond, Pedro G. Ferreira, Guilhem Lavaux, Jens Jasche. Fifth force constraints from the separation of galaxy mass components. Physical Review D, 2018, 98 (6), pp.064015. 10.1103/Phys-RevD.98.064015. hal-01851327

HAL Id: hal-01851327

https://hal.science/hal-01851327v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fifth force constraints from the separation of galaxy mass components

Harry Desmond* and Pedro G. Ferreira

*Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, United Kingdom*

Guilhem Lavaux

*Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis boulevard Arago, 75014 Paris, France and Sorbonne Universités, Institut Lagrange de Paris (ILP), 98 bis boulevard Arago, 75014 Paris, France*

Jens Jasche

*The Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE 106 91 Stockholm, Sweden and Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, D-85748 Garching, Germany*

<figure><img src="image_3.png" /></figure>

(Received 3 July 2018; published 12 September 2018)

One of the most common consequences of extensions to the standard models of particle physics or cosmology is the emergence of a fifth force. While generic fifth forces are tightly constrained at Solar System scales and below, they may escape detection by means of a screening mechanism which effectively removes them in dense environments. We constrain the strength $\Delta G/G_N$ and range $\lambda_C$ of a fifth force with Yukawa coupling arising from a chameleon- or symmetron-screened scalar field—as well as an unscreened fifth force with differential coupling to galactic mass components—by searching for the displacement it predicts between galaxies’ stellar and gas mass centroids. Taking data from the Alfalfa survey of neutral atomic hydrogen (HI), identifying galaxies’ gravitational environments with the maps of [H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, Mon. Not. R. Astron. Soc. 474, 3152 (2018)] and forward modeling with a Bayesian likelihood framework, we find, with screening included, $6.6\sigma$ evidence for $\Delta G > 0$ at $\lambda_C \approx 2$ Mpc. The maximum-likelihood $\Delta G/G_N$ is 0.025. A similar fifth force model without screening gives no increase in likelihood over the case $\Delta G = 0$ for any $\lambda_C$. Although we validate this result by several methods, we do not claim screened modified gravity to provide the only possible explanation for the data: this conclusion would require knowing that the signal could not be produced by “galaxy formation” physics. We show also the results of a more conservative—though less well-motivated—noise model which yields only upper limits on $\Delta G/G_N$, ranging from $\sim 10^{-1}$ for $\lambda_C \approx 0.5$ Mpc to $\sim \text{few} \times 10^{-4}$ at $\lambda_C \approx 50$ Mpc. Corresponding models without screening receive the somewhat stronger bounds $\sim \text{few} \times 10^{-3}$ and $\sim \text{few} \times 10^{-4}$ respectively. We show how these constraints may be improved by future galaxy surveys and identify the key features of an observational program for directly constraining fifth forces on scales beyond the Solar System. This paper provides a complete description of the analysis summarized in [H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, arXiv:1802.07206].

DOI: 10.1103/PhysRevD.98.064015

## I. INTRODUCTION

While ΛCDM has so far passed most of its tests with flying colors, the possibility remains that it requires modification either at levels of sensitivity currently unattainable or in regions of parameter space so far unexplored. In the past few decades, this has motivated a burgeoning program aimed at formulating theoretical alternatives to

concordance cosmology, developing phenomenological frameworks to explore their consequences and devising novel experimental and observational probes.

A generic feature of fundamental extensions to the standard model is the introduction of new dynamical degrees of freedom. These may be produced by substituting dimensional parameters by dynamical fields (e.g., masses [1], dark energy [2] or the gravitational constant [3,4]) and also describe the effects of higher derivatives in the Lagrangian as well as extra dimensions. Any extension

<footnote>*harry.desmond@physics.ox.ac.uk</footnote>
to the Einstein-Hilbert action will necessarily involve new fields [5], and high-energy completions of known physics produce a plethora of new degrees of freedom which may operate on a wide range of scales. Given the ubiquity of novel fields in both gravitational and particle-based extensions to the standard model, seeking new physics is largely tantamount to devising methods to detect additional scalars, vectors or tensors.

Any new field couples naturally to the Ricci scalar R in the gravitational action and hence contributes both to the Universe's energy-momentum and to the expression for the force. The simplest and most common case is a single scalar field with standard kinetic energy, potential V(φ), and coupling φ²R. This modifies the potential generated by a point mass M to

$$ \Phi_{\text{tot}} = -\frac{G_N M}{r} \left( 1 + \frac{\Delta G}{G_N} e^{-mr} \right) = \Phi_N - \frac{\Delta GM}{r} e^{-mr}, \quad (1) $$

where G_N is the bare (Newtonian) gravitational constant, Φ_N is the Newtonian potential, and m and ΔG are new parameters describing the mass of the field and the strength of the new interaction respectively. In theoretical terms, m is set by the potential as m ~ d²V/dφ² while ΔG depends on the strength of the nonminimal coupling and the background value of the field relative to the Planck mass. The bare mass may be altered to the effective mass mₐff given screening (see below), which sets the range of the field as λ_C ≈ 1/mₐff. Equation (1) describes a violation of the inverse-square law by the new field and the introduction of a fifth force (F₅) which adds to the standard gravitational force with a Yukawa form. General relativity (GR) is recovered in the limit ΔG → 0 and also in the case where the field's mass is sufficiently large for the fifth force to be strongly Yukawa suppressed (m → ∞, λ_C → 0). In the case where the field is light (m → 0, λ_C → ∞), the sole effect of the new field is to renormalize the gravitational constant to G ≡ G_N + ΔG. Note that the existence of fifth forces requires only the presence of scalar fields and hence does not necessarily require modifications to the standard model: the Higgs field e.g., leads to a fifth force with a subnuclear range.

The empirical program of probing fifth forces (or, relatedly, violations of the inverse-square law) has been ongoing for around two decades; see [6] for a review. At small scales (~10⁻⁵ – 10⁻² m), the strongest limits are set by the Eöt-Wash experiment which involves precise measurement of the force between two nearby objects [7,8], supplemented by high-frequency torsion oscillators [9] and microcantilevers [10] at even shorter distances. These constraints are ΔG/G_N ≲ 10⁻⁶ at λ_C = 10⁻⁵ m and ΔG/G_N ≲ 10⁻² at λ_C = 10⁻² m. At astrophysical scales, constraints derive from the Earth-Moon distance probed by lunar laser ranging (ΔG/G_N ≲ 10⁻¹⁰ at λ_C = 10⁹ m) and planetary orbits (ΔG/G_N ≲ few × 10⁻⁹ m).

at λ_C = 10¹¹ m) [11,12]. No direct constraint exists at larger scales due to the difficulty of predicting precisely the motions of masses beyond the Solar System, even under the assumption of standard gravity.

Unlike an equivalence principle test, a fifth force search is restricted to ranges λ_C comparable to the size of the system under investigation. Ranges much smaller than the system will go unnoticed as the scalar field will not couple its various parts; conversely, if the range is much larger, the entire system simply feels a renormalized Newton's constant G and the inverse-square law is satisfied. This renders tests which rely on the distance dependence of the force between masses, such as Eöt-Wash and lunar laser ranging, inoperable. However, because a light scalar field affects timelike but not null geodesics it can be probed by means of the Eddington light-bending parameter γ = (1 - ΔG/G_N)/(1 + ΔG/G_N) in the parametrized post-Newtonian framework. γ has been constrained to ≲ 10⁻⁵ by the Cassini satellite [13], requiring ΔG/G_N ≲ 10⁻⁵. This by itself is sufficient to force a putative scalar into astrophysically and cosmologically uninteresting regions of parameter space, as its contribution to the Klein-Gordon equation cannot be more than a fraction O(10⁻³) that of standard gravity.

Nevertheless, there exist ways in which a fifth force may escape these constraints and still have an impact at larger scales. In a number of theories (e.g., generalized scalar-tensor theories and massive gravity) the properties of F₅ depend on the local environment: in high-density regions such as the Milky Way its strength or range becomes small, while it remains operative in the lower density environments of dwarf galaxies and the Universe as a whole. This phenomenon is called “screening” (see [14,15] for reviews) and comes in three basic flavors. The chameleon mechanism [16] arises when the effective mass mₐff becomes dependent on the local density (and thus on ∇²Φ), such that in denser regions mₐff becomes large and the fifth force short range, while in less dense regions (or on cosmological scales) mₐff → 0 and the fifth force reappears. In kinetic screening (e.g., K-mouflage [17]), the scalar interaction is governed by a kinetic function ~(∂φ)², so that φ becomes nonlinear, and hence effectively decouples, when its gradient is large. Finally, in the Vainshtein [18] and symmetron [19] mechanisms, the coupling of the field to matter depends on the local environment, such that ΔG → 0 near massive bodies. In each case, Solar System tests would be expected to probe the screened regime and hence yield the GR result.

In the presence of screening, new scalar fields are possible with Compton wavelengths ranging from microscopic to cosmological values. This renders observationally viable a broad class of gravitational actions, including those that may be expected to form the effective low-energy limits of UV-complete theories. Further, screening holds opens the possibility for cosmic acceleration to receive a
dynamical explanation, potentially circumventing the cosmological constant problem [20], and sanctions the search for fifth force phenomenology on all scales. Crucially, however, tests of screening require unscreened objects, which are the only ones to manifest $F_5$.

The aim of this paper is to constrain fifth forces in galaxies. Not only will this provide direct information on scales much larger than the Solar System (~0.4–50 Mpc), it will also enable constraints in the presence of screening. This is possible due to the wide range of densities and environments of galaxies in the local Universe, including those with masses much lower than the Milky Way and in environments much sparser, which would remain unscreened even for small background values of a novel scalar field.

We restrict ourselves here to chameleon and symmetron screening,<fnref n="1" /> which are conceptually simpler and more amenable to testing through the signal we wish to explore. In these mechanisms, whether or not an object is screened depends on the value of its Newtonian potential $\Phi$ [22]. For a single spherical source this is deducible analytically as a consequence of the fact that chameleon and symmetron charge is contained in a thin shell near the object’s surface (the “thin-shell effect” [16,23]). For multiple sources (and taking into account the possibility of environmental as well as self-screening), numerical simulations show that $\Phi = \Phi_{in} + \Phi_{ex}$ (where $\Phi_{in}$ is the Newtonian potential generated by the object itself at its surface, and $\Phi_{ex}$ that due to its surroundings) remains the best proxy for the degree of screening, as measured by differences between dynamical mass, which is subject to the fifth force, and lensing mass which is not [24–27]. In addition to fifth force strength $\Delta G$ and range $\lambda_C$, then, generic chameleon and symmetron models come also with a critical value of $\Phi$, $\Phi_c$, at which screening kicks in. To good approximation, objects with $|\Phi| > |\Phi_c|$ are in sufficiently dense environments to be screened ($m_{eff} \to \infty$), while those with $|\Phi| < |\Phi_c|$ are not ($m_{eff} \to 0$).<fnref n="2" /> For viable and interesting theory parameters, $|\Phi_c|$ takes values in the range ~$10^{-5} - 10^{-8}$ ($c \equiv 1$), putting some galaxies in the screened and others in the unscreened regime.

Thin-shell screening possesses a range of intragalaxy signals, all of which derive from a common basic effect. Within an unscreened galaxy, gas and dark matter interact via a fifth force with neighboring unscreened mass, which leads to an effective increase in Newton’s constant $\Delta G = 2\beta^2 G$ if the scalar is light and coupled to matter by $\nabla^2\phi = 8\pi\beta G\rho$. (If the galaxy is marginally screened there is an addition

thin-shell factor, but this is not important for the cases of practical interest here.) Main sequence stars, however, are themselves massive objects with surface Newtonian potentials $|\Phi_{in}| \sim 10^{-6}$, and hence for $|\Phi_c| \lesssim 10^{-6}$ they self-screen and do not feel $F_5$. The effective equivalence principle violation [28] caused by the differential sensitivity of stars, gas and dark matter to the fifth force has four observable consequences [29]: (1) a displacement between the stellar and gas disk in the direction of $F_5$, (2) warping of the stellar disk, (3) enhancement of the rotation velocity of the gas relative to the stars, and (4) asymmetry in the stellar disk’s rotation curve. We focus in this paper on the first. By collating neutral atomic hydrogen (HI) and stellar (optical) data for a large sample of galaxies, calculating their gravitational environments and forward modeling the expected signal for given theory parameters $\{\Delta G, \lambda_C, \Phi_c\}$, we statistically constrain a fifth force screened by any thin-shell mechanism. We also provide constraints in the case without screening but where the fifth force is simply coupled differently to stars, gas and dark matter.

The structure of this paper is as follows. In Sec. II we derive the form of the HI-optical offset for a given coupling of $F_5$ to stars and gas, and in Sec. III we describe the survey data that we use to search for it. Section IV details our method: first our mapping of the potential and fifth force acceleration fields across the local Universe; second the calculation from these of the signal expected for given $\{\Delta G, \lambda_C, \Phi_c\}$; and third the likelihood framework we use to constrain the theory parameters by Markov chain Monte Carlo (MCMC). Section V presents our results and forecasts improvements from future data. In Sec. VI we provide further information on our noise model, document the systematic uncertainties in our analysis, compare our work to the literature and offer suggestions for future study. Finally, Sec. VII concludes.

## II. THE SIGNAL: GAS-STAR OFFSETS IN THE PRESENCE OF A FIFTH FORCE

Consider an unscreened galaxy in a chameleon- or symmetron-screened theory subject to an external Newtonian acceleration $\vec{a}$ and a Newtonian acceleration $\vec{a}_5$ due to unscreened matter within $\lambda_C$. Take $|\Phi_c| \lesssim 10^{-6}$ so that stars self-screen. Under the step-function approximation of Eq. (1), unscreened objects with separation $r < \lambda_C$ interact through standard gravity and $F_5$, while those with $r > \lambda_C$ interact only through the former. Thus the unscreened gas and dark matter will effectively couple to the mass generating $\vec{a}_5$ with enhanced gravitation constant $G_N + \Delta G$ and hence feel a total acceleration

$$ \vec{a}_{g,DM} = \vec{a} + \frac{\Delta G}{G_N} \vec{a}_5. \qquad (2) $$

In order for the stars to remain associated with the gas and dark matter, their insensitivity to $\vec{a}_5$ must be

<footnote n="1">¹Similar behavior is expected for the environmentally dependent dilaton [21].</footnote>

<footnote n="2">²An alternative screening criterion (e.g., [24]) is $|\Phi_{in}| > |\Phi_c|$ or $|\Phi_{ex}| > |\Phi_c|$. As this is more difficult to satisfy than $|\Phi_{in}| + |\Phi_{ex}| > |\Phi_c|$, our screened fractions may be considered upper bounds, making our constraints (which benefit from more unscreened galaxies in which the predicted signal is >0) conservative in this regard.</footnote>
compensated by an acceleration due to a separation between their center of mass and that of the dark matter halo. If the stellar centroid lags behind the halo center in the direction of $\vec{a}_5$ by a distance $r_*$, its total acceleration will be

$$ \vec{a}_* = \vec{a} + \frac{G_N M(<r_*)}{r_*^2} \hat{r}_*, \quad (3) $$

where $M(<r_*)$ is the dark matter plus gas mass enclosed by a sphere of radius $r$ around the halo center. Requiring $\vec{a}_{g,DM} = \vec{a}_*$ so that all parts of the galaxy move together, and distinguishing explicitly between the cases where the galaxy itself is screened or unscreened, we find that the equilibrium offset $\vec{r}_*$ satisfies

$$ \frac{M(<r_*)}{r_*^2} \hat{r}_* = \vec{a}_5 \frac{\Delta G}{G_N^2}, \quad |\Phi| < |\Phi_c| \\ \vec{r}_* = 0, \quad |\Phi| > |\Phi_c| \quad (4) $$

where $\Phi$ is the potential of the galaxy as a whole. This is the equation we will solve in Sec. IV, as a function of $\Delta G$, $\lambda_C$ and $\Phi_c$, to calculate the expected signal for the observed galaxies.

Although we have couched our discussion in terms of thin-shell screening, the phenomenon we describe would also result from a generic fifth force, screened or unscreened, that couples differently to different galaxy mass components. In that case, $\Delta G$ in Eq. (4) would be related to that in Eq. (1) by a factor of the normalized coupling coefficient of $F_5$ to gas and dark matter minus that to stars.

### III. OBSERVATIONAL DATA

Our observational catalog of HI-optical displacements is derived from the Alfalfa survey [30–35], a second-generation blind HI survey of ~30 000 extragalactic sources conducted with the Arecibo Observatory. Alfalfa is sensitive to objects of mass $10^6 \lesssim M_{HI}/M_\odot \lesssim 10^{10.8}$ and covers over 7000 deg$^2$ out to $z \simeq 0.06$. More than 98% of HI detections were associated with optical counterparts (OCs) by actively examining DSS2 and SDSS images, in addition to objects in the Nasa Extragalactic Database (NED).

We cut the sources designated as “priors” (quality flag 2), which have signal-to-noise ratio $s \lesssim 6.5$ but are nevertheless matched to OCs with optical redshifts consistent with those derived from HI, as well as those estimated to be high velocity clouds in the Milky Way (quality flag 9). We also remove galaxies with an offset angle between HI and OC centroids of >2 arcmin, which is likely caused by a poor signal-to-noise ratio. These are similar cuts to those employed in previous work along these lines [36]. Objects at larger distance $D$ have larger physical uncertainty in HI-optical displacement for given angular uncertainty $\theta$, and we will find in Sec. V that for realistic $\theta$ values

(see Sec. IV E) objects with $D > 100$ Mpc do not provide further constraining power on $\{\Delta G, \lambda_C, \Phi_c\}$. We therefore cut the catalog at 100 Mpc, yielding 10 822 galaxies in the final sample.

### IV. METHOD

We construct a Bayesian framework for constraining $\{\Delta G, \lambda_C, \Phi_c\}$ by comparing forward-modeled $\vec{r}_*$ values to those observed in Alfalfa galaxies. To do so, we derive or specify probability distributions for the input parameters necessary to solve Eq. (4) and use Monte Carlo sampling to propagate their uncertainties into the predicted $\vec{r}_*$. We thus create a likelihood function for that signal, and score models by applying Bayes’s theorem and a MCMC algorithm to it.

To simplify our inference, we first approximate the exponential in Eq. (1) by a step function at $\lambda_C$,<fnref n="3" /> and second impose the relation between $\lambda_C$ and screening threshold $\Phi_c$ implied by a specific chameleon-screened theory, the Hu-Sawicki [37] model of $f(R)$ gravity [38,39]. In this case, both $\lambda_C$ and $\Phi_c$ are set by the background value of the scalar field, which is equal to the derivative of the function $f$ at the cosmological value $R_0$ of the Ricci scalar $R$, $f_{R0} \equiv df/dR|_{R_0}$ (also the background scalar field value $\phi_0$):

$$ \lambda_C = 32 \sqrt{f_{R0}/10^{-4}} \text{ Mpc}, \quad (5a) $$

$$ |\Phi_c| = \frac{3}{2} f_{R0}, \quad (5b) $$

so that

$$ |\Phi_c| = \frac{3}{2} \times 10^{-4} \left( \frac{\lambda_C}{32 \text{ Mpc}} \right)^2. \quad (5c) $$

Qualitatively, these relations are also of more general applicability with $\lambda_C$ interpreted in terms of the self-screening parameter $\phi_0/(2\beta M_{pl})$, often called $\chi_c$. Although in the most general screening theory $\lambda_C$ and $\Phi_c$ are likely unrelated, we will find the observational data to possess insufficient information for a three-dimensional inference and hence impose Eq. (5) by fiat. We leave it to future modeling with better data to open the parameter space further. The model that we seek to test, then, is fully specified by $\Delta G$ and $\lambda_C$.

For given $\{\Delta G, \lambda_C\}$, we derive the gravitational inputs to the predicted $\vec{r}_*(\Phi$ and $\vec{a}_5$) in Secs. IVA and IV C and the

<footnote n="3">³As our gravitational maps become inaccurate beyond ~200 Mpc, this approximation is necessary for larger fifth force ranges because we cannot include contributions from sources $\gtrsim 2\lambda_C$ from test points ~100 Mpc from us. With future maps extending to higher redshift it will be preferable to treat Eq. (1) more exactly.</footnote>
galactic input [$M(<r)$] in Sec. IV D. Section IV B details our sampling method, Sec. IV F describes our likelihood function, and Sec. IV G presents convergence and consistency tests.

### A. Mapping $\Phi$ and $\vec{a}$

We determine an object's degree of screening by means of its Newtonian potential $\Phi$. This is the sum of an internal component $\Phi_{in}$, due to the object's own mass, and an external component $\Phi_{ex}$ due to its environment. $\Phi_{in}$ can be estimated either from the object's mass and size or from its characteristic velocity. The simplest estimate<fnref n="4" /> is

$$|\Phi_{\text{in}}| = \frac{GM_{\text{vir}}}{R_{\text{vir}}} = V_{\text{vir}}^2. \qquad (6)$$

$\Phi_{ex}$ on the other hand—along with the acceleration $\vec{a}$, which has only an external component<fnref n="5" />—is a nonlocal function of the full density field in space. The remainder of this subsection is devoted to calculating it.

Our starting point is the method of [40] which sums the contributions from three sources: mass associated with survey galaxies, halos hosting galaxies above the survey magnitude limit, and matter situated outside collapsed structures. Briefly, the procedure for determining these three parts is as follows (further details in that paper):

(1) Estimate the total dynamical mass distributions of galaxies in 2M++ [41], an all-sky redshift survey complete to $K = 12.5$ out to $z \approx 0.05$. This is achieved by associating a NFW halo to each 2M++ galaxy using the technique of inverse abundance matching (AM) [42–44]—which assumes a nearly monotonic relation between galaxy luminosity (in this case K-band) and a function of halo mass $M_{vir}$ and concentration $c$, known as the “proxy”—to assign galaxy luminosities to halos in an N-body simulation. We adopt the particular AM model of [45], which uses the proxy $V_{vir}(V_{max}/V_{vir})^\alpha$, with best-fit value $\alpha = 0.6$, and a Gaussian uncertainty

$\sigma_{AM}$ with best-fit value 0.16 dex. This reproduces the observed clustering of galaxies from the simulated clustering of halos (in our case ROCKSTAR [46] halos in the DARKSKY-400 [47] box). The first part of $\Phi_{ex}$ and $\vec{a}$, denoted by subscript “vis”, are then as follows:

$$\Phi_{\text{vis}} = -\sum_{i} \frac{GM_{\text{vir},i}}{r_i} \frac{\ln(1 + \frac{c_i r_i}{R_{\text{vir},i}})}{\ln(1 + c_i) - \frac{c_i}{1+c_i}}, \qquad (7)$$

$$\vec{a}_{\text{vis}} = -\sum_{i} \frac{GM_{\text{vir},i} \frac{1}{r_i} \ln(1 + \frac{c_i r_i}{R_{\text{vir},i}}) - \frac{c_i}{1+c_i} \frac{R_{\text{vir},i}}{R_{\text{vir},i}}}{r_i \ln(1+c_i) - \frac{c_i}{1+c_i}} \hat{r}_i, \qquad (8)$$

for source halo $i$ at distance $r_i$ from the test point. Correct $\Phi_{vis}$ and $\vec{a}_{vis}$ for halos hosting galaxies below the 2M++ detection threshold. We supply K-band luminosities for DARKSKY-400 halos by AM, and calculate their apparent magnitude to partition them into subsets visible and invisible to the 2M++ survey. We calculate the “correction factors” $\vec{c} = \{c_\Phi, c_a, c_\theta\}$, at the locations of DARKSKY-400 objects, required to convert $\Phi_{vis}$ and $\vec{a}_{vis}$ to $\Phi_{hal}$ and $\vec{a}_{hal}$, sourced by all halos with $M_{vir} > 7.63 \times 10^{10} M_\odot$ (the 1000-particle resolution limit of the simulation [48]):

$$\Phi_{\text{hal}} \equiv c_{\Phi} \Phi_{\text{vis}}; \quad |\vec{a}_{\text{hal}}| \equiv c_a |\vec{a}_{\text{vis}}|; \quad c_{\theta} \equiv \frac{\vec{a}_{\text{hal}} \cdot \vec{a}_{\text{vis}}}{|\vec{a}_{\text{hal}}||\vec{a}_{\text{vis}}|}. \qquad (9)$$

To apply this correction in the real Universe, we correlate these factors in the simulation with the observational proxies $d$ (distance to the object), $N_{10}$ (number of objects visible to 2M++ within 10 Mpc) and $\Phi_{vis,10}$ (potential due to objects visible to 2M++ within 10 Mpc). We calculate $\{d, N_{10}, \Phi_{vis,10}\}$ for each real test galaxy and assign it the $\vec{c}$ values of a randomly chosen DARKSKY-400 halo in the same bin of proxy values.

(3) Add the contributions $\Phi_{sm}$ and $\vec{a}_{sm}$ from the remaining matter, not associated with resolved halos, by estimating the mass in linear and quasilinear modes of the density field (subscript “sm” denotes smooth). This is achieved by means of a Bayesian reconstruction of the $z \lesssim 0.05$ density field with the Bayesian Origin Reconstruction from Galaxies (BORG) algorithm [49–53], which propagates information on the number densities and peculiar velocities of 2M++ galaxies assuming $\Lambda$CDM cosmology (with best-fit Planck parameter values) and a bias model. The parameters of the bias model are not known a priori but fitted directly to the data. We use the brand-new particle mesh BORG algorithm [54], which improves in accuracy on the

<footnote n="4">⁴This is how $\Phi_{in}$ was calculated in the simulations which established $\Phi_{in} + \Phi_{ex}$ to be a useful proxy for chameleon screening [24,26,27]. By Eq. (7), however, the Newtonian potential near the centers of halos where galaxies live may be smaller by an order of magnitude or more. Switching to this would reduce screened fractions and boost the predicted signal; using Eq. (6) is therefore conservative. Nevertheless, given more information on the internal dynamics of test galaxies it would be preferable to develop a more precise estimator for the degree of self-screening.</footnote>

<footnote n="5">⁵We assume that galaxies are small and uniform enough, relative to the fifth force range, for tidal fields within them to be unimportant. This may incur error for $\lambda_C$ at the small end of our range of interest (400 kpc): at this point the acceleration of one mass component within a galaxy may receive a relatively significant contribution from others, making $\vec{a}_5$ dependent on the structure of the galaxy itself. We will see in any case that the $\Delta G$ constraints for $\lambda_C \sim 400$ kpc are relatively weak.</footnote>
previous 2-loop perturbation theory determination used in [40].

For our purposes, the output of this algorithm is a posterior distribution of density fields for the local $(677.7h^{-1}$ Mpc)$^3$ volume, with a grid scale of $\Delta r = 2.65h^{-1}$ Mpc. We approximate each field as a set of point masses at the grid cell centers. Although sufficient for $\lambda_C \gtrsim 10$ Mpc, when several grid cells are contained within $\lambda_C$ around each test point, this discretization will generate excessive shot noise at smaller $\lambda_C$ due to the relative sparsity of source masses. Test points equidistant from neighboring source masses may enclose none within $\lambda_C$ and hence artificially receive $\Phi_{ex} = \vec{a} = r_* = 0$. To rectify this, we define a minimum value for the radius out to which source objects are considered, $R_{max}$, of $R_{thresh} = 10$ Mpc $\approx 2.5\Delta r$, which we adopt for $\lambda_C < R_{thresh}$. To convert $\Phi_{sm}(R_{thresh})$ and $\vec{a}_{sm}(R_{thresh})$ to $\Phi_{sm}(\lambda_C)$ and $\vec{a}_{sm}(\lambda_C)$, we use the fact that the region within $\lambda_C$ is small and of nearly uniform density to Taylor expand $\rho$ to lowest order around the position of the test point. Putting the origin at that point,

$$ \Phi(R_{max}) = G \int_{0}^{R_{max}} dV \frac{\rho(\vec{r})}{r} \approx G \int_{0}^{R_{max}} dV \frac{\rho_0}{r} \propto R_{max}^2, \quad (10) $$

$$ \begin{aligned} \vec{a}(R_{max}) &= G \int_{0}^{R_{max}} dV \frac{\rho(\vec{r})}{r^2} \hat{r} \\ &\approx G \int_{0}^{R_{max}} dV \frac{\rho_0 + \vec{\partial}_r \rho|_0 \cdot \vec{r}}{r^2} \hat{r} \\ &= G \vec{\partial}_r \rho|_0 \cdot \int_{0}^{R_{max}} dV \frac{\vec{r}}{r^2} \hat{r} \propto R_{max}^2. \end{aligned} \quad (11) $$

Thus both $\Phi$ and $\vec{a}$ scale with the surface area of the region which sources them.<fnref n="6" /> For $\lambda_C < R_{thresh}$ we therefore set

$$ \Phi_{sm}(\lambda_C) = \frac{\lambda_C^2}{R_{thresh}^2} \Phi_{sm}(R_{thresh}), \quad (12) $$

$$ \vec{a}_{sm}(\lambda_C) = \frac{\lambda_C^2}{R_{thresh}^2} \vec{a}_{sm}(R_{thresh}) \quad (13) $$

and estimate the total (external) potential and acceleration fields as

$$ \Phi_{ex} = \Phi_{hal} + \Phi_{sm}, \quad \vec{a} = \vec{a}_{hal} + \vec{a}_{sm}. \quad (14) $$

## B. Monte Carlo sampling

The calculations of [40] used approximate maximum likelihood estimates of the input model parameters, resulting in one $\Phi_{ex}$ and one $\vec{a}$ field (for given $R_{max}$) over space. However, as anticipated in that work, it is desirable to marginalize over the uncertainties in those parameters to generate a set of plausible gravitational fields, weighted by their respective likelihoods. We implement this here by repeating our calculation of $\Phi$ and $\vec{a}_5$ $N_{MC} = 1000$ times, in each case drawing randomly and independently from the probability distributions of the input parameters to generate a set of final fields. This generates a full probability distribution for $\Phi_{ex}$ and $\vec{a}_5$ at each test galaxy, and hence, by further propagation with Eq. (4), for the expected signal $r_*$. This will allow us to derive posteriors on $\{\Delta G, \lambda_C\}$ by means of a full likelihood framework.

Each of the three steps in Sec. IVA naturally lends itself to this approach by supplying a probability distribution for the relevant parameters, as follows:

(i) AM with nonzero scatter imposes an incomplete constraint on the galaxy-halo connection, in the sense that many halos may plausibly be associated with a given galaxy. We therefore repeat the inverse-AM step 200 times to generate a set of catalogs relating 2M++ galaxies to DARKSKY-400 halos.<fnref n="7" />

(ii) Given values for the observational proxies $d, N_{10}$ and $\Phi_{vis,10}$, a range of correction factors $\vec{c}$ are possible, with probability distributions given by their relative frequencies over resolved DARKSKY-400 objects. Each random draw from their joint distribution corresponds to the association of a given real galaxy with a different simulated object at fixed $\{d, N_{10}, \Phi_{vis,10}\}$. Thoroughly sampling this distribution therefore marginalizes over the possible contributions to $\Phi$ and $\vec{a}$ from halos of mass $M_{vir} > 7.63 \times 10^{10} M_{\odot}$ hosting galaxies that are unresolved by 2M++.

(iii) Uncertainties in the inputs to the BORG algorithm— which include galaxy number densities and bias model—produce uncertainty in the corresponding smooth density field. Because the particle-mesh BORG run that we use produces $\lesssim 10$ independent samples, we marginalize over 10 realizations from the burnt-in part of the chain, chosen further apart from one another than the correlation length.

These complete the first three rows of Table I, which lists the probability distributions of all the parameters of our model. The fourth describes an additional scatter of 10% we impose between the observed and true values of halo

<footnote n="6">⁶We have verified by greatly increasing the resolution of the grid in small subvolumes that these scalings provide estimates of $\Phi_{ex}$ and $\vec{a}$ accurate to within 10% for $\lambda_C$ as low as $\Delta r/20$.</footnote>

<footnote n="7">⁷Note that we do not marginalize over the uncertainties in the AM parameters themselves, which therefore constitute a systematic uncertainty of our analysis. However, due to the tight clustering constraints on these parameters [45], this uncertainty is almost certainly subdominant to the statistical AM uncertainty. We collect and discuss our systematic uncertainties in Sec. VI B.</footnote>
<table><caption>TABLE I. Parameters forming the input to our calculation of the gas-star offset $\vec{r}_*$ for each Alfalfa galaxy, and the distributions from which they are drawn. For each Monte Carlo realization of the fifth force model we sample randomly and independently from these distributions to build an empirical likelihood function $\mathcal{L}(\vec{r}_*|\lambda_C, \Delta G)$. Subscripts *t* and *o* denote true and observed values respectively. Further information is provided in Sec. IV.</caption><thead><tr><th>Parameter</th><th>Description/origin of uncertainty</th><th>Parent distribution/sampling procedure</th></tr></thead><tbody><tr><td>P(M<sub>vir</sub>, c|L; α, σ<sub>AM</sub>)</td><td>Stochasticity in galaxy-halo connection</td><td>200 mock AM catalogs at fixed {α, σ<sub>AM</sub>}</td></tr><tr><td>P(c&#773;|d, N<sub>10</sub>, Φ<sub>vis,10</sub>)</td><td>Location of unresolved halos</td><td>Random draw from frequency distribution in N-body simulation</td></tr><tr><td>Smooth density field</td><td>Galaxy number densities, bias model, etc.</td><td>Ten independent draws from BORG posterior</td></tr><tr><td>P(M<sub>vir,t</sub>, c<sub>t</sub>, r&#773;<sub>t</sub>|M<sub>vir,o</sub>, c<sub>o</sub>, r&#773;<sub>o</sub>)</td><td>Halo properties</td><td>Uncorrelated 10% Gaussian uncertainty</td></tr><tr><td>P(Φ<sub>in,t</sub>|Φ<sub>in,o</sub>)</td><td>Internal potential</td><td>10% (with NSA data) or 30% (without) Gaussian uncertainty</td></tr><tr><td>P(early|M<sub>*</sub>)</td><td>Probability of being early type</td><td>Linear fit to [55], Fig. 5</td></tr><tr><td>P(M<sub>*</sub>, R<sub>eff</sub>, R<sub>d,gas</sub>|M<sub>HI</sub>)</td><td>Galaxy properties from Alfalfa data</td><td>Scalings of [56] with scatters as quoted there</td></tr><tr><td>P(Σ<sup>0</sup><sub>D,t</sub>|Σ<sup>0</sup><sub>D,o</sub>)</td><td>Central dynamical surface density</td><td>0.5 (with NSA data) or 1 (without) dex Gaussian uncertainty</td></tr></tbody></table>

mass, concentration and position, and the remainder will be described later in this chapter. By giving conservative uncertainties to these parameters, we intend the width of the fifth force part of our likelihood function—and hence our inference—to be conservative also.

### C. Relating $\vec{a}$ to $\vec{a}_5$

Section IVA describes the calculation of $\Phi$ and $\vec{a}$ from all mass within a distance $\lambda_C$ of a given test point. $\Phi$ is ready for use in Eq. (4), but the acceleration of interest, $\vec{a}_5$, is due to the fifth force alone, which couples only to unscreened mass. We must therefore determine which source masses are screened and which unscreened, in order to include only the latter. This section gives our method.

First, we calculate $\Phi_{ex}(\lambda_C)$ at the position of each 2M++ galaxy by the method of Sec. IVA. We estimate $\Phi_{in}$ by plugging into Eq. (6) the properties of the halo associated with the galaxy in the inverse-AM step (with a random scatter given by the fifth row of Table I). The halo is unscreened and hence contributes to $\vec{a}_5$, if $|\Phi_{ex}(\lambda_C)| + |\Phi_{in}| < |\Phi_C|$. As $\Phi$ depends on the values of the probabilistic input parameters described in Sec. IV B, we recalculate it for each of the $N_{MC}$ realizations of the underlying model and hence derive a probability $p_j(\lambda_C)$ for halo j to be unscreened. We repeat this also using solely the external part of the potential (i.e., setting $\Phi_{in} = 0$) to make $p_{j,ex}(\lambda_C)$, which we will use below.

With the screening properties of the 2M++ halos determined, we turn next to the remaining mass component: the smooth density field. To ascertain where this is screened, we reconstruct the complete $p_{ex}(\vec{r}; \lambda_C)$ field by linearly interpolating its values $p_{j,ex}(\lambda_C)$ at the positions of the 2M++ halos, which, by virtue of the high completeness of that survey, provide a fair sampling of the entire $z < 0.05$ volume. We then simply evaluate $p_{ex}$ at the position $\vec{r}_k$ of each point mass k representing the smooth density field and assume that the structures in which this mass is located are sufficiently diffuse never to self-screen. $p_{ex}(\vec{r}_k; \lambda_C)$ thus provides the probability for portion k of the smooth density field to be unscreened.

We are now in position to calculate the fifth force acceleration $\vec{a}_5$ for each Alfalfa galaxy. For each Monte Carlo realization, we begin by setting each source halo and portion of the smooth density field to be screened or unscreened according to its probability $p_j(\lambda_C)$ or $p_{ex}(\vec{r}_k; \lambda_C)$ as calculated above.<fnref n="8" /> We calculate $\vec{a}_5(\lambda_C)$, separately for each realization, by evaluating the sum in Eq. (8) over unscreened halos and adding the $GM/r^2$ contribution from each unscreened portion of the smooth density field within $\lambda_C$. As the screened stars do not contribute to $\vec{a}_5$ even in unscreened galaxies, however, we first subtract $M_*$ from $M_{vir}$ in the halo contribution. In practice this modification is minimal since stars comprise $\lesssim 5\%$ of total halo mass.

While calculating $\vec{a}_5$ itself requires sources at most $\lambda_{C,max}$ beyond the maximum distance of any test galaxy, determining whether those sources are themselves screened requires sources a further $\lambda_{C,max}$ out. Since the 2M++ catalog and BORG-reconstructed density fields become unreliable beyond ~200 Mpc, and, as we find below, the fixed angular uncertainty in the Alfalfa HI-optical offset means that test points beyond ~100 Mpc bring little further constraining power to $\Delta G$ and $\lambda_C$, we take $\lambda_{C,max} = 50$ Mpc.<fnref n="9" /> Our minimum $\lambda_C$ value is set by the point at which the predicted $r_*$ values become so small that the corresponding $\Delta G/G_N$ constraint becomes

<footnote n="8">⁸By doing this independently for each object we assume the screening fractions to be uncorrelated. Although not strictly true, this is a good approximation: variations in most of the input parameters affect regions of size $< \lambda_C$, while typical objects are separated by much larger distances.</footnote>

<footnote n="9">⁹By Eq. (5c), $\lambda_C = 50$ Mpc corresponds to $|\Phi_C| = 2.3 \times 10^{-4}$. This is larger than the potential of typical main sequence stars ~$10^{-6}$, suggesting that in this model stars would not be sufficiently dense to screen themselves. Thus both gas and stars would feel $F_5$ and the signal $r_*$ would go away. (This model would also leave the Milky Way, with similar $\Phi$, unscreened.) Nevertheless, we include this model to explore more fully the constraining power of our signal: as in the case when screening is switched off, strictly our constraints require $F_5$ to couple differently to stars and gas.</footnote>
uninteresting. We will show in Sec. V B that $\Delta G/G_N$ cannot be constrained to better than $\mathcal{O}(1)$ for $\lambda_C \simeq 400$ kpc, which is the value we therefore adopt for $\lambda_{C,min}$. This is also the approximate scale at which tidal effects within galaxies may be expected to become important.

### D. Calculating $\vec{r}_*$

Armed with a determination of $\Phi$ and $\vec{a}_5$ for each Alfalfa galaxy, we are now nearing our goal of using Eq. (4) to calculate the expected signal $\vec{r}_*$ as a function of $\Delta G$ and $\lambda_C$. The remaining unknown is $M(<r_*)$, the total enclosed halo plus gas mass within the HI-optical separation $r_*$. Three classes of method present themselves for finding this: one could (1) attempt to fit halo profiles to the galaxies with empirical methods such as AM, (2) use dynamical information at small radius, or (3) employ empirical scalings between mass and light at galaxies' centers. Which of these is most precise depends on the information at one's disposal, as well as the size of the $r < r_*$ region relative to the entire galaxy or halo. We shall see that realistic $\vec{a}_5$ and interesting $\Delta G/G_N$ values within our range of $\lambda_C$ produce $r_*$ values of order 10–100 pc, at least an order of magnitude smaller than galaxy size or halo scale radius. This makes methods such as the first, which rely on global fits to galaxies' luminosity or velocity profiles, unreliable. While central kinematics would provide the most convenient and precise means of calculating $M(<r_*)$, this is not available for most Alfalfa detections, some of which are not even clearly associated with a concentration of optical light. We are therefore forced to resort to the third method.

We estimate $M(<r_*)$ using the observed relation between central baryonic and dynamical surface mass densities reported in [57,58] and parametrized by [59]

$$ \Sigma_D^0 = \Sigma_M S(\Sigma_B^0 / \Sigma_M), \quad (15) $$

where

$$ S(y) = y/2 + y^{1/2}(1 + y/4)^{1/2} + 2\sinh^{-1}(y^{1/2}/2), \quad (16) $$

$\Sigma_D^0$ is the central dynamical surface density; $\Sigma_B^0$ is the central baryonic surface density; and $\Sigma_M$ is a constant, $1.37 \times 10^8 M_\odot/\text{kpc}^2$. This relation captures at small $r$ the empirical trend of increasing dark matter dominance for decreasing baryonic surface density [60–63]: $\Sigma_D^0 \rightarrow \Sigma_B^0 + \Sigma_M$ for $\Sigma_B^0 \gg \Sigma_M$, and $\Sigma_D^0 \rightarrow (\Sigma_M \Sigma_B^0)^{1/2}$ for $\Sigma_B^0 \ll \Sigma_M$. Although possessing some scatter and hence providing a noisier determination of $\Sigma_D^0$ than would direct dynamical data, this method has the advantage of requiring neither spectroscopic information nor assumptions about the form of the overall dark matter distribution.

We take $\Sigma_D^0$ to be the projection of a core with uniform density $\rho_0$ in the very central $r < r_*$ regions of halos, an assumption supported observationally at radii much less than the halo scale $r_s$ [64,65]. Hence

$$ M(<r_*) = \frac{4}{3}\pi r_*^3 \rho_0 = \frac{4}{3}\pi r_*^3 \frac{\Sigma_D^0 - \Sigma_*^0}{2h}, \quad (17) $$

where $h$ is the scale height of the projection and we subtract $\Sigma_*^0$ from $\Sigma_D^0$ to derive the gas plus dark matter mass within the position of the stellar centroid. For realistic galaxy parameters we will find $\rho_0 \sim 10^9 M_\odot \text{kpc}^{-3}$, with a standard deviation between galaxies a factor of a few larger. Under the assumption of constant density the halo acts as a spring, with the restoring force on the stellar centroid proportional to its distance from the halo center. As this force is required to provide a fixed acceleration equal to $\vec{a}_5$, the stiffer the spring (i.e., the greater the central density) the smaller the displacement. Plugging into Eq. (4) and solving for $\vec{r}_*$ yields

$$ \vec{r}_* = \frac{3 \Delta G}{2\pi G_N^2 \Sigma_D^0 - \Sigma_*^0} \frac{h}{\vec{a}_5}, \quad (18) $$

with the following components in the right ascension (J2000, RA, $\alpha$) and declination (DEC, $\delta$) directions:

$$ r_{*,\alpha} = \vec{r}_* \cdot \hat{\alpha}, \quad r_{*,\delta} = \vec{r}_* \cdot \hat{\delta}. \quad (19) $$

These are the quantities we will compare between model and observations. As our observed displacements are restricted to the plane of the sky, when we speak of $\vec{r}_*$ hereafter we will typically mean $r_{*,\alpha}\hat{\alpha} + r_{*,\delta}\hat{\delta}$.

Given Eq. (15), our task now is to deduce $\Sigma_B^0$, the central surface density of total (star plus gas) baryonic mass. This is known to greater precision the more structural galaxy information is available, and we therefore seek to augment our knowledge of the Alfalfa galaxies, which besides HI and OC positions have measured linewidths and total HI masses only, with more detailed optical data. For this purpose we cross-correlate the catalog with the NASA Sloan Atlas (NSA),<fnref n="10" /> a compilation of measurements and derived quantities for nearby galaxies with state-of-the-art sky subtraction and photometric determinations [66]. After cleaning the data by removing galaxies with unreliable size or redshift measurements [67], we take from the NSA stellar mass $M_*$ (derived from kcorrect [68] with a Chabrier [69] initial mass function), Sérsic index SERSIC_N $\equiv n$, aperture velocity dispersion $\sigma$, Sérsic half-light radius along the major axis SERSIC_TH50, and apparent Sérsic minor-to-major axis ratio SERSIC_BA $\equiv (b/a)_{\text{obs}}$.

We associate each Alfalfa galaxy with the nearest NSA galaxy, provided it is within 2 Mpc. This accounts approximately for offsets in the distance estimates recorded in the two catalogs, although we caution that this cross-identification is likely to contain some errors. We show in Sec. V A, however, that the width of the fifth force part of the likelihood function, which measures the uncertainty in

<footnote n="10">¹⁰www.nsatlas.org.</footnote>
the predicted signal and is sensitive to assumptions of galaxy structure, is subdominant to that of the Gaussian component which accounts for measurement uncertainty and non-fifth-force physics. For current data, therefore, precision estimation of $\Sigma_B^0$ is not critical.

Due to the only partly overlapping footprints of the two surveys and the fact that HI and optical flux are often poorly correlated or anticorrelated, our NSA identification succeeds for only 22% of Alfalfa galaxies, and for the remainder we must use less precise means for determining $M(<r_*)$. The following subsections describe our procedure for the cases in which NSA information is or is not available.

# 1. With NSA information

First, we use the observed minor-to-major axis ratio $(b/a)_{obs}$ to estimate the intrinsic axis ratio $(b/a)$. We do this by randomly assigning an inclination $i$ for the galaxy and setting [70]

$$ \left(\frac{b}{a}\right)_{\text{int}}^2 = 1 - \frac{1 - (b/a)^2_{\text{obs}}}{\sin(i)^2} \quad (20) $$

with the sole requirement that $(b/a)_{int} \ge 0.15$, the lowest axis ratio recorded in the NSA. We then calculate the angular diameter distance $d_A$ to each galaxy from the redshift and an assumed cosmology {$h, \Omega_\Lambda, \Omega_m, \Omega_k$} = {0.7, 0.7, 0.3, 0} and calculate the major-axis length as $R_{\text{eff}}^{\text{maj}} = d_A \times \text{SERSIC\_TH50}$. This is related to the circularized, $R_{\text{eff}}$, and minor-axis, $R_{\text{eff}}^{\text{min}}$, half-light radii as

$$ R_{\text{eff}} = (b/a)^{1/2} R_{\text{eff}}^{\text{maj}} \quad \text{and} \quad R_{\text{eff}}^{\text{min}} = (b/a) R_{\text{eff}}^{\text{maj}}. \quad (21) $$

Next, we analytically deproject the spherically symmetric 3D Sérsic profile [71] to derive the stellar surface density and hence calculate the central value $\Sigma_*^0$ [72,73],

$$ \Sigma_*^0 = \frac{M_*}{2\pi n b_n^{-2n} \Gamma(2n) R_{\text{eff}}^2}, \quad (22) $$

where $b_n \equiv 2n - 1/3 + 0.009876/n$. We then calculate the central gas density $\Sigma_g^0$. Although the total HI mass $M_{\text{HI}}$ is measured in Alfalfa, the distribution of this mass is not. We take $M_{\text{gas}} = 1.33 M_{\text{HI}}$, to account for cosmological helium, and assume it to be distributed in an exponential disk with effective radius given by the scaling relation of [56],

$$ \log_{10}(R_{\text{eff,g}}/\text{kpc}) = \log_{10}(0.92 R_{\text{eff}}^{\text{maj}}/\text{kpc}), \quad (23) $$

with a Gaussian scatter of 0.25 dex. Setting $\Sigma_B^0 = \Sigma_*^0 + \Sigma_g^0$ then lets us calculate $\Sigma_D^0$ by Eq. (15). We estimate the scale height $h$ of the dynamical mass distribution as $R_{\text{eff}}^{\text{min}}$.

Finally, we estimate the internal potential as $\Phi_{\text{in}} = -\sigma^2$, with an uncertainty of 10%.

# 2. Without NSA information

If NSA information is not available, it is necessary to estimate not only $h$ and $R_{\text{eff,g}}$, but also $M_*$, $R_{\text{eff}}$ and the characteristic velocity that sets $\Phi_{\text{in}}$. We take $\Phi_{\text{in}} = -V_{\text{max}}^2$, with $V_{\text{max}}$ estimated from the full width at half maximum $W_{50}$ of the radio detection

$$ V_{\text{max}} = \frac{W_t}{2\sqrt{1 - (b/a)^2}}, \quad (24) $$

where

$$ W_t = \sqrt{W_{50}^2 + 4\sigma_d^2} \quad (25) $$

is a turbulence-correction linewidth [74], $\sigma_d = 8$ km/s is a characteristic dispersional velocity [75,76] and $(b/a) = 0.2$ is the assumed intrinsic axis ratio of the system. We assign $V_{\text{max}}$ a 15% Gaussian uncertainty. We then estimate the stellar properties of the galaxy using empirical scaling relations with HI properties [56]:

$$ \log_{10}(M_*/M_{\odot}) = 1.89\log_{10}(M_{\text{gas}}/M_{\odot}) - 8.12 \quad (26) $$

with 0.3 dex scatter, and

$$ \log_{10}(R_{\text{eff}}^{\text{maj}}) = \log_{10}\left(\frac{1}{2} R_0 (M_*/M_0)^\alpha (1 + (M_*/M_0)^\gamma)^{\frac{\beta-\alpha}{\gamma}}\right) \quad (27) $$

with 0.2 dex scatter. To account for first order for the size difference between “early”- and “late”-type galaxies, we assign each Alfalfa galaxy a probability $P_{\text{early}}$ of being early type, derived by interpolating the relation derived by [55] (their Fig. 5) with a linear relation as a function of $\log_{10}(M_*/M_{\odot})$. We use this probability to randomly set the galaxy to be early or late type and assign the corresponding {$\alpha, \beta, \gamma, \log_{10}(M_0), \log_{10}(R_0)$} from Table 1 of [56]. We then use Eq. (23) to estimate $R_{\text{eff,g}}$ and proceed as in Sec. IV D 1.

# 3. Sampling and scatter in $\Sigma_D^0$

As in Sec. IV B, we sample each probability distribution in Secs. IV D 1 and IV D 2 randomly and independently for each Monte Carlo realization of our model, effectively marginalizing over them in our determination of the likelihood distribution of $\vec{r}_*$. This fills in rows 6 and 7 of Table I.

Besides having relatively large scatter, however, these scalings may be expected to provide good fits only to subsets of the total galaxy population, and there is no guarantee that they will accurately describe the Alfalfa sample. In order to ensure that our constraints on $\Delta G$ are nevertheless conservative, we therefore scatter $\Sigma_D^0$ by an
additional 0.5 dex in the case where NSA is available and 1 dex when it is not (final row of Table I). This propagates directly into $r_{*,\alpha}$ and $r_{*,\delta}$ by Eq. (18) and guarantees that the galaxies with NSA information are prioritized in the inference. We find that increasing the scatter of $\Sigma_D^0$ in the model increases predicted $r_*$ values on the whole and hence reduces inferred values of $\Delta G/G_N$. Given that this scatter is itself significantly uncertain, our results for $\Delta G/G_N$ in Sec. V B ought not be considered more than order of magnitude estimates. In Sec. VI B we will show explicitly the result of implementing different scatters for $\Sigma_D^0$. On the other hand, global variations in $\Sigma_D^0$ scatter are scale independent and hence introduce no bias in the $\lambda_C$ dependence of any quantity.

Ideally, unknowns such as this would be treated as free parameters in the inference and marginalized over in the determination of $\Delta G/G_N$ and $\lambda_C$. This may become possible with data of greater quality and quantity in the future.

## E. Modeling the noise

The constraining power of the gas-star offsets for fifth force parameters depends crucially on the angular measurement uncertainty $\theta$. As the OC centroid is typically known to 3 arcsec or better, this uncertainty is dominated by the 21 cm radio detection of HI. This depends on the telescope pointing, data reduction pipeline and signal-to-noise ratio $s$ of the measurement and is not therefore well known a priori. We will determine it empirically from the data itself. Our fiducial method is to set $\theta_\alpha$ and $\theta_\delta$ equal to the dispersion in the RA and DEC components of the measured HI-OC offset itself in bins of $\log_{10}(s)$, as recommended by [32]. In particular, we create 50 uniform bins between the minimum and maximum values 0.66 and 3, and in each bin we calculate separately the standard deviation of $r_{*,\alpha}$ and $r_{*,\delta}$. These have a mean of 18 arcsec over the entire sample and a standard deviation among galaxies of 4 arcsec.

This method guarantees that the measured offsets will be consistent with 0 to ~1$\sigma$ and hence requires that any potential fifth force contribution to the signal be relatively small. Although we will show this to be the case in Sec. V A for all reasonable values of $\Delta G/G_N$, a more self-consistent approach is to parametrize the HI-OC uncertainties and fit for them simultaneously with $\Delta G$ and $\lambda_C$. We present three variants of this approach, each assuming that $\theta$ is correlated only with $s$. In the first, we assume a universal value $\bar{\theta}$ for both $\theta_\alpha$ and $\theta_\delta$ and fit for $\{\Delta G, \lambda_C, \bar{\theta}\}$. In the second and third we model $\theta$ as a polynomial in $\log_{10}(s)$: in one case we fit

$$ \theta = A + B \log_{10}(s), \quad (28) $$

taking A and B as free parameters, and in the other

$$ \theta = C + D \log_{10}(s) + E \log_{10}(s)^2, \quad (29) $$

for $\log_{10}(s) < 1.6$, and $\theta = F$ otherwise, for C, D, E and F free. Each of $A-F$ are given flat priors, with A, C, F > 0 to force $\theta > 0$. We will show in Sec. V B 1 that in each case $\theta$ converges to values consistent with those of the first method and produces an almost identical $\Delta G/G_N$ posterior. This is because by far the greater part of the signal derives from noise for the $\Delta G/G_N$ values of interest in our analysis. Marginalizing over $\bar{\theta}$, $A-B$ or $C-F$ broadens the $\Delta G/G_N$ posteriors only slightly. As discussed further in Sec. VI B, operationally $\theta$ also includes the contribution to the signal from other physics, under the assumption that it supplies a Gaussian component to the likelihood with width a function only of $s$. The degeneracy between measurement uncertainty and additional physics in the determination of $\theta$ will become important for future HI surveys with greater spatial resolution (see Sec. V D).

We will also show results for a choice we call “conservative,” in which the $\theta$ values of the fiducial method are doubled. This yields robust upper limits on $\Delta G/G_N$ and is the choice we focused on in [77]. Note however that this model produces significantly lower overall likelihood values than the fiducial one due to overprediction of the dispersions of the HI-OC displacements at each $s$ and in a model comparison sense is therefore clearly disfavored.

## F. Likelihood model

Over the $N_{MC}$ realizations of our probabilistic model, and for given $\Delta G$ and $\lambda_C$, we have now derived a distribution of predicted $r_{*,\alpha}$ and $r_{*,\delta}$ values for each Alfalfa galaxy. Heuristically these constitute the likelihood distribution of the signal, and by scoring it against the observations and applying Bayes’s theorem one can derive constraints on $\Delta G$ and $\lambda_C$ by MCMC. In practice we proceed as follows.

We begin by constructing a template signal, assuming that each test galaxy is unscreened and taking $\Delta G/G_N = 1$ and one of 20 values for $\log_{10}(\lambda_C/Mpc)$ uniformly spaced between $\log_{10}(0.4)$ and $\log_{10}(50)$. For each test galaxy, we find the minimum ($r_0$) and maximum ($r_1$) of $r_{*,\alpha}$ and $r_{*,\delta}$ separately over the $N_{MC}$ realizations of the model. We then partition the results of these realizations into $N_{bins} = 10$ uniform bins between these limits and define $P_j \equiv N_j/N_{MC}$, where $N_j$ is the number of realizations falling in bin $j$. We take this histogram to represent the unscreened component of the model’s likelihood function in $\{r_{*,\alpha}, r_{*,\delta}\}$ space. We calculate also the fraction $f$ of realizations in which the galaxy in question is unscreened. Given that the model predicts $r_* = 0$ for the screened remainder, and treating the orthogonal RA and DEC components as independent, the overall likelihood $\mathcal{L}_i(r_{*,\alpha}, r_{*,\delta})$ for test galaxy $i$ to have $\bar{r}_*$ components $r_{*,\alpha}$ and $r_{*,\delta}$ is given by

$$ \mathcal{L}_i(r_{*,\alpha}, r_{*,\delta} | \Delta G, \lambda_C) = \mathcal{L}_i(r_{*,\alpha} | \Delta G, \lambda_C) \times \mathcal{L}_i(r_{*,\delta} | \Delta G, \lambda_C), \quad (30) $$
where

$$ \mathcal{L}_i(r_{*,\alpha}|\Delta G, \lambda_C) = (1-f)\delta(r_{*,\alpha}) + f\sum_{j=0}^{N_{\text{bins}}} P_j \delta(r_{*,\alpha} - \Delta G(r_{0,\alpha} + j\Delta r_{\alpha})). \quad (31) $$

Here $\Delta r_{\alpha} = (r_{1,\alpha} - r_{0,\alpha})/(N_{\text{bins}} - 1)$ is the width of each bin, and we have approximated the unscreened component of the likelihood as a discrete sum of delta functions $\delta$ at the midpoints of the bins. $f, P_j, r_0$ and $r_1$ are implicit functions of $\lambda_C$: by Eq. (5c) and the procedure of Sec. IV, a larger $\lambda_C$ corresponds to a larger $|\Phi_c|$ and hence a larger unscreened fraction $f$ and signal limits $r_0$ and $r_1$. Since $r_*$ is proportional to $\Delta G$ [Eq. (4)], this parameter simply multiplies the bin positions. The likelihood for the DEC component is defined analogously.

We convolve this with a Gaussian probability distribution function for these true values given observed values $r_{*,\alpha,\text{obs}}$ and $r_{*,\delta,\text{obs}}$ and their measurement uncertainties $\sigma_i$. Dropping the dependence on $\Delta G$ and $\lambda_C$ for brevity,

$$ \mathcal{L}_i(r_{*,\alpha,\text{obs}}) = \int dr_{*,\alpha} \mathcal{L}_i(r_{*,\alpha}) \times \mathcal{L}_i(r_{*,\alpha,\text{obs}}|r_{*,\alpha}), \quad (32) $$

where

$$ \mathcal{L}_i(r_{*,\alpha,\text{obs}}|r_{*,\alpha}) = \frac{\exp\left\{-\frac{(r_{*,\alpha,\text{obs}} - r_{*,\alpha})^2}{2\sigma_{i,\alpha}^2}\right\}}{\sqrt{2\pi}\sigma_{i,\alpha}}. \quad (33) $$

The observational uncertainty is given by

$$ \sigma_{i,\alpha} = d_{A,i} \cos(\delta_i) \theta_i, \quad (34) $$

where $d_{A,i}$ is the angular diameter distance to galaxy $i$, $\delta_i$ is the galaxy's declination, and $\theta_i$ is the angular uncertainty assigned by either the fiducial or conservative method of Sec. III.<fnref n="11" /> The corresponding DEC uncertainty is given by

$$ \sigma_{i,\delta} = d_{A,i} \theta_i. \quad (35) $$

Performing the convolution in Eq. (32) then yields

$$ \begin{aligned} \mathcal{L}_i(r_{*,\alpha,\text{obs}}|\Delta G, \lambda_C) = & (1-f) \frac{\exp\{-\frac{r_{*,\alpha,\text{obs}}^2}{2\sigma_{i,\alpha}^2}\}}{\sqrt{2\pi}\sigma_{i,\alpha}} \\ & + f \sum_{j=0}^{N_{\text{bins}}} P_j \frac{\exp\{-\frac{(r_{*,\alpha,\text{obs}} - \Delta G r_j)^2}{2\sigma_{i,\alpha}^2}\}}{\sqrt{2\pi}\sigma_{i,\alpha}}, \end{aligned} \quad (36) $$

where $r_j \equiv (r_{0,\alpha} + j\Delta r_{\alpha})$ is the position of the jth bin center. Treating the galaxies as independent, the likelihood of the entire data set $d$ is then

$$ \mathcal{L}(d|\Delta G, \lambda_C) = \prod_i \mathcal{L}_i(r_{*,\alpha,\text{obs}}|\Delta G, \lambda_C) \mathcal{L}_i(r_{*,\delta,\text{obs}}|\Delta G, \lambda_C). \quad (37) $$

Finally, we apply Bayes's theorem to deduce the probability of $\Delta G$ and $\lambda_C$ themselves:

$$ P(\Delta G, \lambda_C | d) = \frac{\mathcal{L}(d | \Delta G, \lambda_C) P(\Delta G, \lambda_C)}{P(d)}, \quad (38) $$

where $P(\Delta G, \lambda_C)$ is the prior distribution of the model parameters and $P(d)$ is the constant probability of the data for any $\{\Delta G, \lambda_C\}$. We take a uniform prior in $\log_{10}(\lambda_C/\text{Mpc})$ in the range $\log_{10}(0.4) - \log_{10}(50)$, and the improper prior $\Delta G \ge 0$ flat in $\Delta G$.<fnref n="12" /> We will find our likelihood function to provide clear upper bounds on $\Delta G$, obviating the need to impose such a bound a priori.

Deriving the likelihood function directly from the results of the Monte Carlo realizations in this way removes the need to assume a particular functional form for it. This reduces bias in the inference, and we have checked explicitly that our method accurately reconstructs $\Delta G/G_N$ when mock data generated with a particular value are fed in (Sec. V C). Since several of the most important input probability distributions in Table I are non-Gaussian, we find that assuming the overall likelihood to be Gaussian (i.e., using the mean and standard deviation of the $N_{\text{MC}}$ predicted $r_{*,\alpha}$ and $r_{*,\delta}$ values rather than the full histogrammed distribution) biases inference of a known input $\Delta G/G_N$ by a factor of $\sim 2$.

## G. Convergence and consistency tests

We have taken 10 realizations of the BORG smooth density field and 200 of the AM galaxy-halo connection. We verify that these numbers are sufficient to fully propagate these probability distributions into the predicted signal by directly investigating the impact of changing them. While with fewer realizations the spread in $r_{*,\alpha}$ and $r_{*,\delta}$ in the signal templates rises with the number of realizations, indicating additional uncertainty in the marginalization, this effect levels off near the values chosen. Using a larger number of realizations would not therefore appreciably impact the results. We have also checked that using a different set of 10 smooth density fields and 200 AM galaxy-halo connections does not alter any of our conclusions.

Similarly, we have used a finite number of Monte Carlo draws from the model to estimate the likelihood function.

<footnote n="11">¹¹As explained in Secs. IVE and VD, $\theta$ also models a Gaussian component of the likelihood model due to non-fifth-force physics. We refer to it as measurement uncertainty here simply for brevity.</footnote>

<footnote n="12">¹²We find a Jeffrey's prior to yield almost identical results.</footnote>
With $N_{\text{MC}} \lesssim 1000$, repetitions of the entire procedure yield constraints on $\Delta G/G_N$ and $\lambda_C$ (Sec. V) that differ at the $\gtrsim 20\%$ level. That this is not true with $\gtrsim 1000$ realizations indicates that in this case the probability distributions of the input parameters of Table I have been thoroughly sampled and their uncertainties fully marginalized over in the final constraints. We check also that varying our estimates of $\Phi_{\text{in}}$, $M_{\text{vir}}$, $c$, $\Sigma_D^0$, $M_*$, $R_{\text{eff}}$ and $R_{\text{d,gas}}$ within a factor of 2 does not significantly affect the results; this change is subdominant not only to the measurement uncertainties, but also to the other theoretical uncertainties deriving from the smooth density field, galaxy-halo connection and location of low-mass halos. (The scatter in $\Sigma_D^0$ is however important; see Sec. VI B.)

We check that our inference is not driven by potential outliers by excising up to 10% of points with the most extreme values of $\langle r_{*,\alpha} \rangle$, $\langle r_{*,\delta} \rangle$, $r_{*,\alpha,\text{obs}}$ or $r_{*,\delta,\text{obs}}$, where angle brackets denote a modal average over the $N_{\text{MC}}$ model realizations. This has little impact on the results. In Sec. V C we will show explicitly that jackknife or bootstrap resamples of the Alfalfa data set behave similarly.

Finally, we check the convergence of our MCMC chains (run with the affine-invariant Python sampler `EMCEE`) by visual inspection and by means of the Gelman-Rubin statistic. The low dimensionality of the inference combined with the relatively simple form of the likelihood function and unimodality of the posteriors makes exploring the parameter space relatively easy. With 10 walkers, we find 300 burn-in steps to suffice for 1D inference ($\Delta G/G_N$ at fixed $\lambda_C$), and 600 burn-in steps to suffice for 2D inference (both $\Delta G/G_N$ and $\lambda_C$).

We caution that errors in the probability distributions of the model parameters (including the implicit choice of delta-function priors for parameters which should in reality have some uncertainty) may lead to systematic error. We discuss this further in Sec. VI B.

# V. RESULTS

## A. Comparison of observed and predicted signal

To begin our investigation of the presence of a fifth force signal in the Alfalfa data, we show in Fig. 1(a) (red) the distribution of the predicted values of $\langle r_* \rangle = \langle (r_{*,\alpha}^2 + r_{*,\delta}^2)^{1/2} \rangle$ for fiducial model parameters $\Delta G/G_N=1$ and $\lambda_C = 5$ Mpc. For each galaxy this is the mode of the binned distribution of the $N_{\text{MC}}$ model realizations, i.e., $\langle r_* \rangle = r_0 + j_0 \Delta r$ where $j_0$ maximizes $P_j$. For comparison, we show in green the case in which screening is switched off. This is tantamount to setting $|\Phi_c| \to \infty$, implying $f=1$ for all Alfalfa galaxies and that all matter within $\lambda_C$ of each test point is fully unscreened. Thus any pair of masses within that separation interacts with a fifth force, leading to an enhanced signal. We show the analogous distribution for the real data in Fig. 1(b). In Fig. 1(c) we show the uncertainty on the predicted $r_*$, defined as the

minimal width enclosing 68% of the model realizations. As in Fig. 1(a) we show the results both with and without screening. Figure 1(d) shows the distribution of uncertainties $\theta_i$ in the HI centroids in the real data, assigned by the fiducial method.

Figure 2(a) shows the correlation of $\langle r_* \rangle$ with $\Phi$. The green points are for the case in which screening is switched off, the red points for when it is included—with the signal of screened ($f < 0.5$) galaxies set to 0—and the blue points are the same as the red except without the $f < 0.5$ points being lowered. The black dashed line shows the screening threshold $\Phi_c$ for this choice of $\lambda_C$, and the median error bars in both directions, enclosing 68% of the model realizations, are shown in the top left. Figure 2(b) shows the analogous correlation with $a_5$ with identical color scheme, and Figs. 2(c) and 2(d) show the corresponding correlations for the real data. In Figs. 3(a) and 3(b) we plot the RA and DEC components of the observed and predicted signals directly against each other.

We draw four inferences at this stage:

(1) For this choice of model parameters, the predicted signals are typically $O(1)$ kpc while the observed ones are $O(10)$ kpc with a narrower distribution. This reflects the strength of the restoring force on the stellar centroid due to its offset from the halo center [derived from the central density in Eq. (17)], on the one hand, and the relatively large observational uncertainties in the HI centroid position on the other. This indicates that were the observed $r_*$ values due entirely to a fifth force, the required $\Delta G/G_N$ would be $O(10)$. If competitive constraints are to be placed on $\Delta G$ with this signal, therefore, it must primarily be by means of the correlation of the relative directions and magnitudes of $\vec{r}_*$ and $\vec{a}_5$ across the galaxy sample. Assuming $\Delta G/G_N \lesssim 1$, modified gravity can account for at most $\sim 10\%$ of the HI-OC offset.

(2) For $\Delta G/G_N = 1$, the uncertainty on the predicted signal (width of the likelihood function) is typically $\sim 10^{-2} - 10^{-1}$ kpc, while that of the observations is again $O(10)$ kpc. As these two terms contribute equally to the likelihood in Eq. (36), this implies that the strength of the constraints we set in Sec. V B is determined predominantly by the precision $\theta$ with which the HI centroid of the Alfalfa galaxies is located. As the width of the likelihood function scales with $\Delta G$ [Eq. (36)], this is even more pronounced at the small $\Delta G/G_N$ values [$O(10^{-2})$] of interest there.

(3) Both $\langle r_* \rangle$ and $\sigma(r_*)$ are slightly larger for the unscreened than screened model. This is because $\vec{a}_5$ is larger in the unscreened case due to contributions from additional masses in the vicinity of a given test point. Note however that this does not take into account the fact that in the model with screening
<figure><img src="image_4.png" /><figcaption>FIG. 1. Histograms of the expectation (upper) and uncertainty (lower) of the fifth force predicted (left) and observed (right) gas-star offset $r_*$ over all Alfalfa galaxies. For the fifth force prediction we take the mode of the 1000 Monte Carlo realizations for each galaxy of the screened model with $\lambda_C = 5$ Mpc and $\Delta G/G_N = 1$, given 10 bins for both the RA and DEC components ($r_{*,\alpha}$ and $r_{*,\delta}$) between their minimum and maximum values. The uncertainty is the minimal width enclosing 68% of the realizations. For the observed signal the uncertainty is assigned using the fiducial method of Sec. IV E. $r_{*,\text{obs}}$ is typically several times larger than $\langle r_* \rangle$, and the corresponding uncertainty two orders of magnitude larger, even for the relatively extreme case in which the fifth force is as strong as gravity ($\Delta G = G_N$). The majority of the signal must therefore be due to non-fifth-force effects.</figcaption></figure>

the test galaxies themselves may be screened and hence have a predicted $r_*$ of 0 (the results here are only for the unscreened fraction). In terms of our inference this will prove to be the most significant difference between the two models.

(4) By construction $\langle r_* \rangle$ shows clear and characteristic dependence on $\Phi$ and $a_5$. The predicted signal is proportional to $a_5$ by Eq. (4), generating a positive correlation whether or not screening is included. Without screening, the positive correlation between $|\Phi|$ and $a_5$ (which are sourced by the same mass distribution) simply induces a positive correlation also between $r_*$ and $|\Phi|$. The introduction of screening has two effects. First, some source masses surrounding test galaxies in high-density regions become screened and hence no longer contribute

to $\tilde{a}_5$, causing the predicted signal to drop: this is the difference between the blue and green points at $|\Phi| > |\Phi_c|$ in Fig. 2(a). That this offset is small indicates that in this regime $\Phi$ is set mostly by the test galaxies' own masses rather than their environment, so that most surrounding mass is unscreened and hence contributes to $a_5$ and $r_*$ even when screening is included. Second, test galaxies with $|\Phi| > |\Phi_c|$ become themselves screened, causing their predicted signal to fall to 0 (red). These galaxies may span the entire range of $a_5$ [Fig. 2(b)], although they are more likely to be found at higher values; the largest predicted signals in the case with screening therefore tend to be produced in galaxies very close to being screened. On the contrary, $r_{*,\text{obs}}$ shows no apparent correlation with $\Phi$ or $a_5$.
<figure><img src="image_5.png" /><figcaption>FIG. 2. Correlations of predicted (upper) and observed (lower) signal with Newtonian potential $\Phi$ (left) and magnitude of fifth force acceleration $a_5$ (right). The prediction is calculated as in Fig. 1, again for the case $\lambda_C = 5$ Mpc, $\Delta G/G_N = 1$, and $\langle\Phi\rangle$ and $\langle a_5 \rangle$ are the modal averages over the 1000 Monte Carlo model realizations. In the top row, the red points are for the full chameleon- or symmetron-screened model, the blue points are the same as the red except without the $r_*$ values of screened galaxies set to 0, and the green points show the case where screening is entirely switched off. The black dashed line shows the threshold $|\Phi_c|$, above which galaxies are screened. The median sizes of the error bars in the four directions—minimal 68% bounds for the prediction, $1\sigma$ for the observations—are shown by the red crosses in the top left corners (the uncertainties in $\langle r_* \rangle$ are the same size as the red point). While the predictions show a clear cutoff in $r_*$ for $|\Phi| > |\Phi_c|$ and a positive correlation with $a_5$ [Eq. (4)], this is not apparent in the observations.</figcaption></figure>

This again implies that the greater part of the observed signal must be of non-fifth-force origin.

These estimations of the fifth force effect are relatively coarse. A precise quantitative determination of the preference of the data for $\Delta G > 0$—and of the constraints on $\Delta G$ and $\lambda_C$ that it implies—requires the full likelihood formalism of Sec. IV. It is to this that we now turn.

## B. Constraints on $\Delta G$ and $\lambda_C$

### 1. Effect of noise model

From Sec. VA it is clear that for any reasonable value of $\Delta G/G_N$ modified gravity accounts for at most a small fraction of the total gas-star offset, justifying the fiducial

method of Sec. IV E for setting the measurement uncertainties $\theta_i$. This ought to manifest in an agreement between the results of that model and those of the others in that section where $\theta$ is parametrized and fitted for. Here we show this explicitly for the screened model with $\lambda_C = 1.8$ Mpc.<fnref n="13" /> Figure 4 shows the posteriors of $\{\Delta G/G_N, \bar{\theta}\}$ in the case where $\theta = \bar{\theta}$ arcsec is a free parameter, and $\{\Delta G/G_N, A, B\}$ and $\{\Delta G/G_N, C, D, E, F\}$ for the models of Eqs. (28) and (29) respectively. These are to be compared with Fig. 5, which shows the $\Delta G/G_N$

<footnote n="13">¹³We focus on this value of $\lambda_C$ because it maximizes the overall likelihood, as will be shown in Sec. VB 3.</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 3. Correlation of the observed and predicted components of $r_*$ in the RA ($r_{*,\alpha}$, left) and DEC ($r_{*,\delta}$, right) directions, with the predicted signal again calculated as in Fig. 1 and error bars as in Fig. 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. For the screened model with $\lambda_C = 1.8$ Mpc, posteriors of $\Delta G/G_N$ and additional parameters specifying the angular width $\theta_i$ of the Gaussian, non-fifth-force component of the likelihood function for galaxy i. (a) $\theta_i = \bar{\theta}$ arcsec for all galaxies. (b),(c) $\theta_i$ determined as a function of signal-to-noise ratio $s$ of the detection according to Eqs. (28) and (29) respectively. In each case the chains converge to the same solution for $\Delta G/G_N$ as when $\theta$ is set a priori to the standard deviation of the measured $r_*$ in bins of $\log_{10}(s)$ (Fig. 5), indicating that our inference is robust to variations in the noise model.</figcaption></figure>

posterior for the fiducial noise model in which $\theta$ is fixed a priori to equal the standard deviation of $r_{*,\text{obs}}$ in bins of $\log_{10}(s)$. The same solution for $\Delta G/G_N$ is picked out in each case, and in Fig. 4(a) it is apparent that the maximum-likelihood $\bar{\theta}$ is close to the average scatter in $r_{*,\text{obs}}$ across the entire data set (18 arcsec). The model of Eq. (29) picks out almost the same $\theta(s)$ as Eq. (1) in [32]. We have checked that these results hold at other $\lambda_C$ values also. For simplicity we therefore restrict ourselves from now on to the fiducial model in which $\theta$ is fixed from the outset.

## 2. 1D inference without screening

For our first main analysis we switch screening off, fix $\lambda_C$ at one of its 20 values between 400 kpc and 50 Mpc and

constrain $\Delta G$ alone. We find the $\Delta G/G_N$ posterior to be peaked at or near 0 for all $\lambda_C$, indicating no significant deviation from GR. We show in green the $1\sigma$ upper limit on $\Delta G/G_N$ as a function of $\lambda_C$ in Fig. 6, three example posteriors in Fig. 7, and the maximum increase in log-likelihood over the case $\Delta G = 0$ in Fig. 8. As may be expected, the constraints are tighter for larger $\lambda_C$: $\Delta G/G_N \lesssim 10^{-4}$ ($1\sigma$) for $\lambda_C = 400$ kpc, and $\lesssim$few $\times 10^{-3}$ for $\lambda_C = 50$ Mpc. This is because there are more source masses within a larger $\lambda_C$, leading to an enhanced $a_5$ and hence predicted $r_*$, and therefore requiring $\Delta G/G_N$ to be smaller for consistency with the fixed observations. Assuming there is indeed no preference for $\Delta G > 0$ in this model, the bumps in the green lines of Figs. 6 and 8 indicate the noise level of our experiment.
<figure><img src="image_8.png" /><figcaption>FIG. 5. Posterior of $\Delta G/G_N$ for $\lambda_C = 1.8$ Mpc under the fiducial noise model, with and without screening.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. $1\sigma$ upper limit on $\Delta G/G_N$ ($\Delta G_{1\sigma}$), as a function of $\lambda_C$, using the fiducial uncertainties for the model without screening (green) and the conservative uncertainties for the model with screening (red). Switching to the conservative uncertainties makes only a small ($\sim 0.1$ dex) difference in the unscreened case. The results with screening are weaker at small $\lambda_C$ first because a sizable fraction of the test galaxies are screened, and hence have predicted $r_*=0$, and second because fewer neighboring masses contribute to $\bar{a}_5$.</figcaption></figure>

A limit of particular interest for the case without screening is $\lambda_C \to \infty$, which corresponds to a light scalar field as found e.g., in Brans-Dicke theory. Extrapolating our constraints to high $\lambda_C$, our analysis may be expected to require $\Delta G/G_N \lesssim 10^{-4}$ in this case, which corresponds to $\omega \gtrsim 500$ in Brans-Dicke. Although not competitive with Solar System results [78], this constraint derives from a very different gravitational regime. Were $f(R)$ gravity (which requires $\Delta G/G_N = 1/3$) not to employ screening, it would be ruled out to $1\sigma$ across our parameter space, requiring $f_{R0} < 10^{-8}$.

It is important to remember, however, that the physical interpretation of these no-screening results is complicated by the fact that the signal itself relies on a difference between the fifth force interactions of stars and gas, which most naturally arises when the former are screened and the latter not. Absent screening entirely, the constraints here would require the scalar field simply to couple to the gas and dark matter but not the stars. Our analysis can however be readily extended to the case of arbitrary differential coupling to galaxy mass components. In that case $\bar{a}_5$ in Eq. (4) must be multiplied by the difference in the relative coupling coefficient of the stars and gas, and the constraints on $\Delta G/G_N$ are simply weakened by that factor.

### 3. 1D inference with screening

With screening switched on, $\lambda_C$ values in the range $\sim 0.4-4$ Mpc prefer $\Delta G > 0$. We show in red in Fig. 7(a) the range of $\Delta G/G_N$ posteriors, and in Figs. 8(a) and 8(b) the maximum log-likelihood increase over GR and best-fit $\Delta G/G_N$ values respectively as a function of $\lambda_C$. The preference for a screened fifth force reflects a typically positive correlation between $\bar{a}_5$ and $\vec{r}_{*,\text{obs}}$. In particular, the data would seem to favor $\Delta G/G_N \sim \mathcal{O}(10^{-2} - 10^{-1})$ at $\lambda_C \sim 1$ Mpc ($f_{R0}, \Phi_c \approx 10^{-7}$), with the maximum-likelihood value $\lambda_C = 1.8$ Mpc, $\Delta G/G_N = 0.025$ (also shown in Fig. 5). This corresponds to $|\Phi_c| = 4.7 \times 10^{-7}$ by Eq. (5c), and is not ruled out by any previous analysis. If such a force were realized in nature, neighboring $\lambda_C$ values would also be

<figure><img src="image_10.png" /><figcaption>FIG. 7. Posteriors on $\Delta G/G_N$ for three values of $\lambda_C$ under the fiducial noise model, with and without screening. In the left panel, the posterior for the model without screening is at very small $\Delta G/G_N$. For $\lambda_C \lesssim 3$ Mpc a nonzero $\Delta G$ is statistically preferred, but only when screening is included.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 8. Left: Maximum increase in log-likelihood over the case $\Delta G = 0$, for any $\Delta G/G_N$, as a function of $\lambda_C$ in the range 0.4–50 Mpc. We use the fiducial uncertainties and show the results both with (red) and without (green) screening. The model without screening achieves maximum log-likelihood values at most a few larger than GR, which we believe to indicate the noise level of our inference. The model with screening however achieves a maximum increase of 16 for $\lambda_C = 1.8$ Mpc, with a coherent trend over similar models differing only minorly in the fifth force range. Right: Best-fit $\Delta G/G_N$ value corresponding to the maximum-likelihood point at each $\lambda_C$. Larger $\Delta G/G_N$ values are inferred at smaller $\lambda_C$ where the radius is lower, within which mass contributes to the fifth force, and hence $\vec{a}_5$ and the predicted $r_*$ are smaller.</figcaption></figure>

expected to be preferred over GR simply because for a given galaxy $\vec{a}_5$ is a continuous function of $\lambda_C$, and hence the predicted $r_*$ values over nearby $\lambda_C$'s are correlated. This is the qualitative trend exhibited in Fig. 8(a); in Sec. V C we show this trend to be quantitatively as expected also. That no model without screening achieves a meaningful increase in likelihood over GR indicates both that our noise model is robust and that the phenomenon we are observing depends crucially on the environment and internal structure dependence of screening. This seems hard to mock up by other means. We test our possible detection in detail in following sections.

The results for the conservative measurement uncertainties (twice fiducial) are shown as the red curve in Fig. 6. That the $\Delta G/G_N$ constraint is degraded when screening is introduced reflects the fact that when source objects are screened they contribute neither to $a_5$ nor therefore to $r_*$, while when test objects are screened they automatically receive $r_* = 0$. The weaker the predicted signal at fixed $\Delta G$, the larger $\Delta G$ may be before it becomes statistically discrepant with the observations. The effect of screening is greater when $|\Phi_c|$ is lower, which by Eq. (5c) is at lower $\lambda_C$. Conversely, for $\lambda_C \gtrsim 10$ Mpc the majority of source and test points are unscreened in any case, so removing screening has little effect. Note that for $\Delta G/G_N = 1/3$, as in $f(R)$, our conservative error choice requires $\lambda_C \lesssim 500$ kpc, which corresponds to $f_{R0} \lesssim 2 \times 10^{-8}$.

### 4. 2D inference with screening

Finally, we perform the full 2D inference for both $\Delta G$ and $\lambda_C$ in the presence of screening. Under the fiducial

noise model the posterior simply picks out the maximum-likelihood solution $\lambda_C = 1.8$ Mpc, with corresponding $\Delta G/G_N = 0.025 \pm 0.003$ [Figs. 9(a)–9(c)]. Under the conservative uncertainties [Figs. 9(d)–9(f)] the preference for low $\lambda_C$ is simply a volume effect: by allowing $\Delta G/G_N$ to span a wider range, a smaller $\lambda_C$ boosts the contribution from the prior. Even with these artificially inflated uncertainties, however, we are able to place the strong constraint $\Delta G/G_N < 0.95$ ($3\sigma$) for any $\lambda_C$ in the range 0.4–50 Mpc.

### C. Validation

Figure 8(a) appears to show a clear preference for $\Delta G > 0$ for $\lambda_C \simeq 1–4$ Mpc in the case where screening is included. Here we document the tests we have performed to validate this detection, including the use of mock data with a $\Delta G/G_N$ value injected by hand, jackknife and bootstrap resampling, and rotation of the measured displacement on the plane of the sky.

#### 1. Mock data with $\Delta G/G_N = 0$

We begin by generating 100 mock data sets with $\Delta G = 0$, so that $r_{*,\text{obs},i,\alpha/\delta} \sim \mathcal{N}(0, \theta_{i,\alpha/\delta})$. By giving the offset a random direction on the plane of the sky we ensure its complete noncorrelation with $\vec{a}_5$ over the sample. We repeat our inference for each mock data set with the maximum-likelihood model $\lambda_C = 1.8$ Mpc, relaxing the prior $\Delta G \ge 0$ and refitting $\theta_{i,\alpha/\delta}$ from the spread in the $r_{*,\alpha/\delta}$ values in bins of $\log_{10}(s)$ (Sec. IV E). In each case we calculate the magnitude of the deviation (in $\sigma$) of the best-fit $\Delta G/G_N$ value from 0 by dividing the median of the
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 9. Posteriors of $\Delta G/G_N$ and $\lambda_C$ for the full 2D inference, using either the fiducial (top) or conservative (bottom) noise model. The left and center panels show the marginalized constraints while the right panels show the full 2D posteriors. With the fiducial uncertainties the MCMC picks out the maximum-likelihood point of Fig. 8(a), $\lambda_C = 1.8$ Mpc; as this corresponds to a single bin the posterior of $\lambda_C$ is flat. With the conservative uncertainties the preference for low $\lambda_C$ is purely a volume effect arising from the uniform prior, since lower $\lambda_C$ permits a larger range of $\Delta G/G_N$.</figcaption></figure>

posterior by the standard deviation. The histogram over all mock data sets is shown in Fig. 10. As expected this forms roughly a standard normal distribution, confirming that the offset from $\Delta G = 0$ inferred from the real data—$6.6\sigma$— would be highly unlikely in the frequentist sense to arise from the noise alone.

## 2. Mock data with maximum-likelihood $\Delta G/G_N$

Next we generate mock data with a fifth force signal injected by hand. This provides an independent test of the validity of the fifth force interpretation by checking for biases in the reconstructed value of $\Delta G/G_N$. We generate 250 mock data sets at the maximum-likelihood point of our analysis—$\lambda_C = 1.8$ Mpc, $\Delta G/G_N = 0.025$—by scattering the expected $r_*$ values$^{14}$ by the fiducial Gaussian noise in bins of $s$. We refit each one with the $\lambda_C = 1.8$ Mpc model (including recalculation of $\theta_{i,\alpha/\delta}$, as above) and calculate

<figure><img src="image_18.png" /><figcaption>FIG. 10. Distribution of deviations of best-fit values of $\Delta G/G_N$ from 0 (median of $\Delta G/G_N$ posterior divided by standard deviation) derived by fitting a model with $\lambda_C = 1.8$ Mpc to mock data sets with $r_{*,\text{obs},\alpha/\delta}$ randomly scattered around 0 by the uncertainties. This forms a roughly standard normal distribution, as expected for a healthy noise model. By contrast, the deviation for the real data at this $\lambda_C$ is $6.6\sigma$, indicating that it is very unlikely to have been generated by a model with noise only.</figcaption></figure>

<footnote>$^{14}$As in Sec. VA, we take this to be the mode of the 1000 Monte Carlo model realizations in the case that the galaxy is likely unscreened ($f > 0.5$) and 0 otherwise. We have checked that using the mean instead does not lead to appreciably different results.</footnote>
<figure><img src="image_19.png" /><figcaption>FIG. 11. Results of refitting a screened model with $\lambda_C = 1.8$ Mpc to 250 mock data generated by that model, differing only in noise under the fiducial noise model. Panel (a) compares the distribution of best-fit $\Delta G/G_N$ values to the input value (0.025), demonstrating negligible bias. Panel (b) shows that the increase in log-likelihood over the case $\Delta G = 0$ is strongly correlated with the $\Delta G/G_N$ value inferred. This increase is however considerably smaller than that achieved in the real data [shown in panel (c)], indicating that the model is not fully sufficient to account for the observations.</figcaption></figure>

both the maximum-likelihood $\Delta G/G_N$ and the increase in log-likelihood $\Delta \log(\mathcal{L})$ that value achieves over the GR case $\Delta G = 0$. The results are shown in Fig. 11. We find the reconstructed $\Delta G/G_N$ values to be centered around the input value (vertical red line), indicating that our likelihood function picks out a known truth without bias. (We have checked that this holds also for mock data generated by models of different $\lambda_C$.) The considerable variation between mock data sets is due to the dominance of the random noise.

We show in Fig. 11(b) that the $\Delta \log(\mathcal{L})$ values are strongly correlated with best-fit $\Delta G/G_N$, in the sense that mock data sets in which a stronger signal is inferred achieve a larger likelihood increase over GR. The recovered $\Delta \log(\mathcal{L})$ values are however a factor ~4 smaller on average than that in the real data [Fig. 11(c)]. This indicates that the relative contribution of signal and noise to our mocks does not perfectly reflect that in the observations: the mock samples have a relatively weaker signal than Alfalfa. This almost certainly derives from an inadequacy in our noise model, most likely its failure to account for a correlation of the non-fifth-force component of $\vec{r}_*$ with anything other than $s$. This will need to be improved in the future when better priors on both the HI centroid measurement uncertainty and dependence of $\vec{r}_*$ on “galaxy formation” physics are available. Other factors which may contribute to disagreement in Fig. 11(c) include an inaccurate $\Phi_c - \lambda_C$ relation [Eq. (5c)] strictly holds only for Hu-Sawicki $f(R)$] and hence to a bias in galaxies’ screening fractions as a function of fifth force range or to inaccuracies in the galaxy and halo properties input into the model. We discuss these issues further in Sec. VI.

We have checked also that when the inference is repeated using models with a range of $\lambda_C$ values for the same mock data, the greatest $\Delta \log(\mathcal{L})$ is achieved at the input value $\lambda_C = 1.8$ Mpc. The peak in $\Delta \log(\mathcal{L})$ is however narrower with $\lambda_C$ than in the data, so that neighboring $\lambda_C$ models have a discrepancy in $\Delta \log(\mathcal{L})$ even greater than a factor of 4.

### 3. Bootstrap and jackknife resampling

Next we repeat our analysis with 350 bootstrap- or jackknife-resampled mock data sets drawn from the Alfalfa data. For the jackknife case we retain a random 70% of the galaxies. Again for the model with $\lambda_C = 1.8$ Mpc, we calculate the discrepancy in $\sigma$ between the $\Delta G/G_N$ posterior and $\Delta G = 0$, and we plot the resulting histograms, along with the corresponding maximum-likelihood $\Delta G/G_N$ values, in Fig. 12. The vertical red lines indicate the results in the full Alfalfa sample. On the whole the mock data sets have a smaller significance due to the reduced statistics, especially in the jackknife case where the sample size is lower. Nevertheless, the majority of resamples in both cases have a discrepancy with $\Delta G = 0$ of >3$\sigma$ and a best-fit $\Delta G/G_N$ value close to 0.025, indicating that our conclusions do not depend sensitively on peculiar features of the Alfalfa sample.

### 4. Rotation of the signal

Finally, we repeat the analysis with the predicted HI-OC displacements (i.e., direction of $\vec{a}_5$) rotated through 90° or 180° on the plane of the sky. This is to check that the directions of the observed and predicted signal are indeed positively correlated, as must be the case if part of the signal has a fifth force origin. The inference of $\Delta G > 0$ should be weakened when one of the vectors is turned through 90° and eliminated when turned through 180°. Figure 13 shows that this is indeed the case (cf. the unrotated result in Fig. 5). We have checked that similar results hold for different $\lambda_C$.

In a companion piece we present a similar detection from a largely independent data set by analyzing warps in galactic disks [79].

### D. Forecasting constraints from future surveys

We now investigate the potential of our constraints for further improvement in the future, focusing solely on the
<figure><img src="image_20.png" /><figcaption>FIG. 12. Left and center: Distributions of deviations as in Fig. 10 (for the model with $\lambda_C = 1.8$ Mpc), but for data sets derived by bootstrap or jackknife resampling the Alfalfa sample. Each jackknife data set contains 70% of the full sample, and the vertical red lines show the deviation in the real data (6.6σ). Although on the whole both sets imply $\Delta G > 0$ with lower significance than the full data set due to reduced statistics, the fact that the majority of resamples achieve a $\gtrsim 3\sigma$ detection indicates that the preference for a screened fifth force is not a peculiar property of the Alfalfa sample. Right: Maximum-likelihood $\Delta G/G_N$ values inferred from the bootstrap and jackknife resamples. These cluster around the Alfalfa result, shown as the vertical red line.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIG. 13. $\Delta G/G_N$ posteriors at $\lambda_C = 1.8$ Mpc for the case in which the predicted $\vec{r}_*$ of each galaxy is turned through (a),(b) 90° in either direction on the plane of the sky or (c) 180°. In the first and third cases the preference for $\Delta G/G_N > 0$ is eliminated, indicating no alignment or antialignment of the measured $\vec{r}_*$ with the rotated $\vec{a}_5$. In the second a residual correlation remains, although it is weaker than the unrotated case as evidenced by the lower maximum-likelihood $\Delta G/G_N$ and tail towards small values relative to Fig. 5. This demonstrates that the observed $\vec{r}_*$ points on the whole in the direction of $\vec{a}_5$.</figcaption></figure>

conservative error choice which yields only upper bounds on $\Delta G/G_N$. To begin, we recall from Sec. VA (Fig. 1) that the uncertainties in the measured $\vec{r}_*$ are larger by over two orders of magnitude on average than the theoretical uncertainties due to the potential and acceleration fields and test galaxy structure. This implies that the best way to strengthen the constraints is to improve the HI survey; only once the HI uncertainties have been greatly reduced will significant further improvement come from increasing the precision of $\Phi$, $\vec{a}_5$ and $M(<r_*)$.

Our constraint on fifth force strength $\Delta G/G_N$ therefore depends principally on two parameters, the angular uncertainties $\theta_i$ in the gas-star offsets and the number $N_{\text{gal}}$ of galaxies observed. While the latter is a property purely of the galaxy survey, the former encapsulates not only the measurement uncertainty in the HI centroid position but also the contribution to the likelihood from non-fifth-force physics. The relative importance of these is not currently

known: it may be determined either by increasing the HI resolution—if the HI-OC offset scales in the same way it derives predominantly from measurement uncertainty—or by estimating the magnitude of offset that various phenomena in $\Lambda$CDM may be expected to induce. We intend to pursue the second avenue in future work by means of hydrodynamical simulations. In this section, however, we adopt the most optimistic assumption that the entire HI-OC offset is due to measurement uncertainty and may therefore be arbitrarily reduced by future surveys. In practice constraints will be weaker than this due to underestimation of the scatter in HI-OC offset for given HI spatial resolution.

We calculate the $\Delta G/G_N$ constraint as a function of $\theta_i$ and $N_{\text{gal}}$ by treating them as variables in the generation of mock data sets with $\Delta G = 0$, as described in Sec. VC. To begin, we multiply the conservative $\theta_i$ values by a universal scaling factor $10^{-3} < \Theta < 1$. This implies an average HI uncertainty of $\langle \theta \rangle = 36 \Theta$ arcsec. Since we only have
<figure><img src="image_22.png" /><figcaption>FIG. 14. Contour plot forecasting the $1\sigma \log_{10}(\Delta G/G_N)$ constraints obtainable with a data set of $N_{\text{gal}}$ galaxies and average angular HI resolution $\langle\theta\rangle = \Theta \times 36$ arcsec, for the screened model with $\lambda_C = 5$ Mpc. We assume the most optimistic case for the noise, in which the entire non-fifth-force part of the signal derives from measurement uncertainty in the position of the HI centroid.</figcaption></figure>

locations for existing Alfalfa galaxies, we consider the dependence of $\Delta G$ on $N_{\text{gal}}$ for $N_{\text{gal}} < 10, 822 \equiv N_{\text{Alf}}$, the number of galaxies in the Alfalfa data set. We retain a fraction $10^{-3} < f_N < 1$ of these galaxies, chosen at random, and sample $\Theta$ and $f_N$ at eight logarithmically equally spaced intervals; to forecast the constraints achievable for $N_{\text{gal}} > N_{\text{Alf}}$ we will extrapolate these results.

For a given mock data set specified by $\Theta$ and $f_N$, we calculate the $1\sigma$ limit on $\Delta G/G_N$ at the fiducial value $\lambda_C = 5$ Mpc. To sample both the noise in the mock data and the specific galaxies retained we repeat this procedure

10 times to find an average $\Delta G/G_N$ constraint, $\Delta G_{1\sigma}$. We show the result as a contour plot in Fig. 14, and in Figs. 15(a) and 15(b) we show the variation of $\Delta G_{1\sigma}$ with $\Theta$ and $N_{\text{gal}}$ separately (for a fixed value of the other parameter) for three example values.

We find $\Delta G_{1\sigma}$ to have approximately power-law dependence on $\langle\theta\rangle$ and $N_{\text{gal}}$, and therefore we fit to it a function of the form

$$ \Delta G_{1\sigma} \approx a \left(\frac{10^3}{N_{\text{gal}}}\right)^b \left(\frac{\langle\theta\rangle}{1 \text{ arcsec}}\right)^c, \quad (39) $$

finding best-fit values $\{a, b, c\} = \{8.6 \times 10^{-4}, 0.91, 1.00\}$. We show these fits as the solid lines in Fig. 15. $\Delta G_{1\sigma}$ therefore falls more rapidly with $N_{\text{gal}}$ than predicted by the scaling $\Delta G_{1\sigma} \sim 1/\sqrt{N_{\text{gal}}}$, indicating greater sensitivity to sample size than would apply for the case of $N_{\text{gal}}$ random draws from a stochastic model under GR. The dependence on $\Theta$ in linear. This fit may be used to forecast constraints on $\Delta G/G_N$ for any future survey: e.g. (again subject to the proviso that the entire scatter in HI-OC displacement is due to measurement uncertainty), $\{\langle\theta\rangle, N_{\text{gal}}\} \approx \{10^{-1} \text{ arcsec}, 10^8\}$, achievable with the “mid” configuration of SKA1 [80,81], should probe $\Delta G/G_N$ to the $10^{-9}$ level. Comparing to Fig. 4 of [6], we see that at this stage fifth force constraints by our method will be comparable to lunar laser ranging and planetary probes at much smaller scales; they will also be competitive with those from proposed Solar System tests such as laser ranging to Phobos and optical networks around the Sun [82]. Even before SKA, its pathfinders APERTIF [83] and ASKAP [84] will provide significant further constraining power over Alfalfa. We caution however that further modeling work will be necessary to extend our gravitational maps to

<figure><img src="image_23.png" /><figcaption>FIG. 15. Variation of $1\sigma \Delta G/G_N$ constraint, $\Delta G_{1\sigma}$, achievable with a data set of average angular resolution $\langle\theta\rangle$ and size $N_{\text{gal}}$, for fixed values of the other parameter as indicated in the legend. The solid lines show the fits of Eq. (39).</figcaption></figure>
the distance that this $N_{\text{gal}}$ requires ($z \sim 0.5$) and potentially to model time variation in parameters such as $f_{R0}$. Given future data sets with more constraining power it will also be desirable to drop or generalize the constraint of Eq. (5c) in order to treat more general cases of screening.

As we account for the overall statistical impact of sample size but not a systematic shift in gravitational environment or galaxy mass, accuracy of our extrapolated constraints requires the locations and structures of galaxies of lower HI luminosity, observed by future surveys, to be similar on average to those already measured by Alfalfa. As fainter galaxies are less likely to self-screen and may also tend to live in sparser environments, our forecasts are conservative in this regard.

# VI. DISCUSSION

## A. Noise model

The evidence for a screened fifth force depends crucially on the non-fifth-force part of the likelihood function, modeled by a Gaussian with angular width $\theta_i$. This is poorly known a priori. While average values of $\theta_i$ may be fitted as part of the model, leading to results consistent with our fiducial analysis, correlations with variables other than signal-to-noise ratio could bias our inference. It is however hard to see how this could eliminate the fifth force signal we have inferred: this would require that $\theta_i$ correlate not only with galaxies' environments similarly to $\vec{a}_5$, but also with their degree of screening (recall that no model without screening achieves an increase in maximum-likelihood over GR) and internal structures in the manner of $\vec{r}_*$ in Eq. (4). Indeed, our fiducial errors are already conservative in setting the signal-to-noise ratio of the detection of any offset between stars and gas to ~1, a valuable safeguard against erroneously imputing the signal a modified gravity origin in a framework that neglects more mundane phenomena. No additional effect is therefore necessary to explain the observations: the analysis could easily prefer $\Delta G = 0$ for all $\lambda_C$. While our conservative error choice is useful for sidestepping the question of $\theta_i$ and hence achieving robust upper bounds on $\Delta G/G_N$, its inadequacy is clear from its severe overprediction of the dispersion of the measured displacements and resultant low overall likelihood.

We have been forced to derive the measurement uncertainties empirically due to insufficient prior information on the precision of the HI centroid location: it is simply too difficult to propagate uncertainties in the telescope pointing, data reduction pipeline and HI-optical cross-correlation, which themselves are poorly known. Future HI surveys, such as those afforded by interferometric observatories such as SKA, may provide better control of these. However, even if information of this nature cannot be used to set informative priors on the HI uncertainties, the increased precision in the centroiding will greatly improve

the strength of the test performed here. As demonstrated in Sec. V D, simply repeating our analysis on a SKA-class sample would probe $\Delta G/G_N$ to $O(10^{-9})$ for $\lambda_C \sim O(1-10 \text{ Mpc})$, easily sufficient to confirm or rule out our putative detection regardless of the precise manner in which uncertainties are modeled.

## B. Galaxy formation physics and other systematics

We have been careful to incorporate known uncertainties in the input parameters of our model into probability distributions for those parameters and hence marginalize over them in the prediction of $\vec{r}_*$. Provided the true parameter values are contained within the priors, these uncertainties cannot cause systematic error in our inference but will simply inflate our posteriors and lead to conservative bounds. Nevertheless, there remain a number of inputs, some of them implicit, whose true values may not lie within our priors. We discuss these here in roughly decreasing order of importance.

Our model assumes that the true $r_{*,\alpha}$ and $r_{*,\delta}$ values (after accounting for noise) arise purely from fifth force effects. It is however highly likely that these could be nonzero even in the case $\Delta G = 0$. This may arise from a number of baryonic processes within a galaxy that affect gas and stars differently, e.g., hydrodynamical drag, ram pressure and the transfer of energy and momentum by stellar feedback. A realistic model must therefore contain a contingent for $r_* \neq 0$ when $\Delta G = 0$ by introducing further parameters describing these baryonic effects. This is present to first order in our model, as the principal component of the likelihood function is simply a Gaussian in $r_{*,\alpha}$ and $r_{*,\delta}$ with a width $\theta$ that correlates with the signal-to-noise ratio $s$ of a galaxy's detection (Sec. IV E). Absent modified gravity, larger HI-optical offsets are likely to be caused by baryonic physics, which a larger uncertainty $\theta$ would effectively absorb: in each bin of $s$ the baryonic contribution to the signal adds to the measurement uncertainty in quadrature. The information we glean on modified gravity comes primarily from the relative directions of $\vec{r}_*$ and $\vec{a}_5$, the correlation of their magnitudes across the sample, and the dependence of the screened fraction on the gravitational environment, quantified here by $\Phi$.

Nevertheless, it is important to check that baryonic effects could not generate an environment dependence of $\vec{r}_*$ similar to the effect of screening. This may be done by applying our inference framework to mock data sets generated by hydrodynamical simulations in which these effects are included. An informative test however requires better than $O(\text{kpc})$ resolution in the central regions of halos, which rules out the majority of current-generation cosmological simulations [85–88]. Conversely, zoom simulations do not have the halo numbers required for statistically significant results. Only a few present simulation sets may therefore be of use (e.g., [89–92]). As these
simulations differ in the subgrid models required at the small scales of interest here, it will ultimately be necessary to check several of them. Although the manner in which hydrodynamics affects the internal and environmental properties of halos is at least qualitatively known (e.g., [93–95]), the specific signal in which we are interested has not yet been investigated.

If such effects in ΛCDM are not sufficient to account for the correlations we attribute here to a screened fifth force, our constraints on ΔG/G$_N$ must be further validated by examining the impact of baryonic physics in a Universe governed by modified gravity. While such simulations have begun in recent years, and confirm that baryonic and gravitational effects on large scale structure and galaxies’ internal dynamics are nonlinearly coupled and hence not fully orthogonal (e.g., [96,97]), they focus on fifth force strengths at least an order of magnitude larger than those that constitute our maximum-likelihood models.

We have assumed the model specified by {$ΔG, λ_C$} to produce a cosmology consistent with ΛCDM, specifically in terms of its prediction for the smooth density field at $z = 0$ and the properties of the halo population. While chameleon- or symmetron-screened modified gravity does alter the expansion history of the Universe and growth rate of structure, and hence the halo mass function and internal structure of halos, its effects in the region of parameter space of most interest here ($ΔG/G_N \lesssim 0.1$, $λ_C \lesssim 5$ Mpc) are small [98–101]. They are likely subdominant to the uncer- tainties we do model in the smooth density field, galaxy-halo connection and properties of the individual halos them- selves. Our method should not therefore be thought a means of testing modified gravity in cosmology—and has little if any relevance to cosmic acceleration—but rather as simply a means of seeking fifth forces, of gravitational or nongravita- tional origin, on the scales of galaxies and their environ- ments. Screened modified gravity also affects stellar luminosities and hence galaxy mass-to-light ratios [102] as well as the relative kinematics of stars and gas [103], but again these effects are small for $λ_C \sim 1$ Mpc and would not be expected to skew our inference.

As discussed in [40], our incorporation of both halos and a non-halo-based mass contribution leads us to overesti- mate the total density on $\gtrsim 10$ Mpc scales by $\sim 10\%$. This will cause a similar overestimation of $|Φ|$, biasing the screened fraction of our test points high and the template signal low. This leads to a weakening of the ΔG constraints, making our limits conservative. While this effect is partly mitigated by a corresponding overestimation of $\vec{a}_5$, we showed in [40] that the acceleration is more sensitive than the potential to mass on smaller scales around the test point, which is less susceptible to double counting.

Although we marginalize over the possible galaxy-halo connections for fixed AM proxy and scatter, we do not marginalize over the distributions of these parameters themselves or otherwise consider the possibility that this

model inadequately describes our source galaxy popula- tion. Given that AM parameters are fairly well constrained by clustering studies, however, this extra uncertainty is unlikely to be significant.

Factors that affect the overall magnitude of the predicted signal [e.g., $Σ_D^0$, $h$ and $|\vec{a}_5|$ in Eq. (18)] are degenerate with the inferred value of ΔG/G$_N$. As mentioned in Sec. IV D 3, a particularly important quantity is the scatter in central dynamical surface density $Σ_D^0$ around the result of Eq. (15), which affects the magnitude of the predicted $r_*$ and hence our ΔG/G$_N$ constraints by a relatively large factor. In Fig. 16 we plot the scatter in $Σ_D^0$ (for the galaxies without NSA information; the scatter for those with NSA informa- tion is pinned to half this value) against the maximum- likelihood value of ΔG/G$_N$ inferred under the fiducial uncertainties for the screened model with $λ_C = 1$ Mpc. The slope of the best-fit line is -3.8, indicating that a change to the $Σ_D^0$ scatter of 0.26 dex is sufficient to shift the best-fit ΔG/G$_N$ by an order of magnitude. Were the evidence for ΔG > 0 confirmed by independent signals or data, it would become important to refine model inputs such as these to better determine ΔG/G$_N$. Consistency between signals under a generic fifth force scenario could even be used to constrain them.

By contrast, the $λ_C$ dependence of the maximum- likelihood and ΔG/G$_N$ constraint is degenerate only with factors that scale with $λ_C$ in the same way as the predicted $\vec{r}_*$. As this condition is far harder to satisfy due to the

<figure><img src="image_24.png" /><figcaption>FIG. 16. Dependence of maximum-likelihood ΔG/G$_N$ on the scatter in central dynamical surface density $Σ_D^0$ [Eqs. (15)–(18)], in the case where NSA information is not available, for the screened model with $λ_C = 1$ Mpc. Our fiducial scatter is 1 dex, as quoted in Table I. The ΔG/G$_N$ constraint is highly sensitive to this choice and differs by more than an order of magnitude across the range of a priori reasonable values. This makes our inference of ΔG/G$_N$ quite uncertain and the formal uncertainties in Figs. 5 and 7 potentially misleading. The trend of best-fit ΔG/G$_N$ with $λ_C$ however is not sensitive to this uncertainty because it derives from the scale dependence of the predicted signal rather than its overall normalization and scatter.</figcaption></figure>
peculiar dependence of $\vec{r}_*$ on $\Phi$ and $\vec{a}_5$ (and the further dependences of those on both halo structure and environment), we consider the trend of $\Delta \log(\mathcal{L})$ with $\lambda_C$ to be considerably more robust than the constraint on $\Delta G/G_N$ itself. For example, altering the $\Sigma_D^0$ scatter within the range of Fig. 16 does not significantly affect Fig. 8(a).

### C. Comparison with the literature

We begin by comparing our work with previous studies at the scales of interest here (0.4–50 Mpc; Sec. VIC 1) before broadening our discussion to include tests of related models at larger and smaller scales (Sec. VIC 2).

#### 1. On astrophysical scales

Much of the groundwork for constraining modified gravity with intragalaxy data was laid in [29], which provided the first calculations of the functional forms and sizes of the expected signals. These authors estimated the magnitude of $\vec{r}_*$ by taking the external field $\vec{a}_5$ to be generated by a single halo towards which the test galaxy is falling. Assuming the test halo has constant density, they argue that $r_*$ could typically be expected to be ~1 kpc for $\Delta G/G_N = 1$. This is quantitatively confirmed by Fig. 1(a). We note however that our estimates of $r_*$ employ significantly more complex and realistic assumptions than those of [29]. First, the fifth force accelerations for our test galaxies are derived from their full environments rather than a single nearby object (which is itself liable to be screened for lower values of $|\Phi_c|$). Second, we take only the central regions of the test halos to be cored and calculate their densities from the central baryon surface density rather than by averaging over the entire halo. As halo density presumably falls with $r$ outside the core, this leads us to assign a larger central density than [29], reducing the expected $r_$. However, the significant scatter that we apply to the central densities (1 dex in most cases) causes many halos within each model realization to have lower density and hence greater predicted $r_*$.

The study most similar to ours is that of [36] (specifically their Sec. IV), whose authors also searched for evidence of screening in the HI-optical offsets of Alfalfa detections. Our work expands on this in three main ways.

(1) In [36], the fifth force acceleration and central mass $M(<r)$, and hence the expected $\vec{r}_*$, are not known. This precludes forward modeling the displacements from theory parameters and hence constraining them quantitatively. In addition, the direction of $\vec{r}_*$ on the plane of the sky is of no use if the direction of $\vec{a}_5$ is unknown. These authors must therefore rely instead on the average difference in $r_*$ between screened and unscreened subsamples as a summary statistic of the impact of the fifth force. Our likelihood formalism on the other hand uses the full information on a galaxy-by-galaxy basis to quantitatively constrain

the theory parameters. In particular, the use of both components of $\vec{r}_*$ on the plane of the sky brings significant additional constraining power.

(2) Our model for the gravitational field includes not only the mass in halos hosting bright objects, but also the mass in fainter halos and the interhalo component of the density field [40]. This eliminates pathological cases in which $\Phi_{ex}$ is apparently 0 due simply to incompleteness in the maps of [24] (e.g., [103], Table 1). It comes further with full point-by-point probability distribution functions for $\Phi$ and $\vec{a}_5$, which enable us to precisely model uncertainties in the predicted signal.

(3) Our sample size is hugely larger than [36]'s: 10 822 vs 245 galaxies. This is partly due to the use of a more complete Alfalfa catalog and partly because our gravitational maps, based on 2M++ rather than SDSS, cover the entire sky. We are not therefore forced to impose a cut on angular position. As shown in Figs. 14 and 15(b), a large sample size is critical for beating down the bound on $\Delta G/G_N$.

Reference [36] constrains the average $r_*$ due to modified gravity to be $\lesssim 1$ kpc. Given our results on the average size of the signal in Sec. VA, this corresponds to $\Delta G/G_N \lesssim 1$ for $\lambda_C = 5$ Mpc. This is significantly weaker than our own constraints due to the extra machinery of point (i) and increased sample size of point (iii), despite the increased screened fractions resulting from the extra mass of point (ii). While our constraints are limited by the observational uncertainty of the HI centroids, those of [36] were limited primarily by the sophistication of their theoretical modeling. Indeed, our analysis outperforms their forecast that a predicted offset of 0.5 kpc would be detectable with 24 500 dwarf galaxies at $\langle\theta\rangle = 24$ arcsec: we are sensitive to a predicted offset of $\sim 0.01$ kpc ($r_* \sim 1$ kpc for $\Delta G/G_N = 1$ and $\Delta G/G_N$ constraints $\sim 0.01$) using 11 000 galaxies—with no selection on density—at $\langle\theta\rangle = 36$ arcsec.

The strongest constraints previously derived on astrophysical scales come from comparing screened and unscreened distance indicators to nearby galaxies. For $|\Phi_c| \lesssim 10^{-6}$ distances derived from screened tip-of-the-red-giant-branch stars would be expected to differ from those of unscreened Cepheids; their consistency implies $\Delta G/G_N \lesssim 0.2 (1\sigma)$ for $|\Phi_c| > 10^{-6}$ and $\Delta G/G_N \lesssim 1$ for $|\Phi_c| = \text{few} \times 10^{-7}$ [104]. Similar results are achievable comparing the gas and stellar rotation curves of isolated dwarf galaxies [103]. Despite a weak signal and large observational uncertainties, our inference is several times stronger due to a combination of huge sample size, great range of environments probed, and use of a vector rather than scalar observable which effectively affords two orthogonal signals in the plane of the sky. No previous analysis has reached the level of sensitivity required to probe the region $\{\lambda_C \approx 1.8 \text{ Mpc}, \Delta G/G_N \approx 0.025\}$ in which we find evidence for a signal.
Additional dynamical galactic signals which have not yet been used to place quantitative constraints on fifth forces include warping of stellar disks, asymmetries in disk kinematics and offsets in the rotation curves of stars and gas [29,36]. Galaxy warps are the subject of [79], and the other effects will be studied in future work.

## 2. On smaller and larger scales

Testing fifth forces on astrophysical scales is a relatively new endeavor: other direct constraints come only from signals at much smaller scales, while cosmological tests are capable of producing independent though indirect and weaker limits. As discussed in Sec. I, unscreened fifth forces with ranges $10^{-14} < \lambda_c/m < 10^{-6}$ are strongly constrained by torsion balance experiments, lunar laser ranging and planetary data within the Solar System. Although our $\Delta G/G_N$ constraints are weaker than these by up to 4 orders of magnitude, reflecting the complex nature of astrophysical systems, they provide the possibility of filling in the parameter space of possible distance-dependent (or potential-, acceleration- or curvature-dependent) modifications to GR, which incorporate a fifth force, at the galaxy scale [105].

The strongest constraint from Solar System tests for $\lambda_c$ as large as that of interest here comes from the Eddington light bending parameter measured by Cassini, which requires $\Delta G/G_N \lesssim 10^{-5}$ [13]. Our constraint in the absence of screening is 1–2 orders of magnitude worse and also requires a differential coupling of the fifth force to stars, gas and dark matter. Nevertheless we use a qualitatively different signal, and, as discussed in Sec. V D, our method is straightforwardly scalable with future data and may be expected to reach the sensitivities of Solar System tests with next-generation surveys.

Without a means of shielding the scalar field from the screening influence of the surrounding dense environment, screened theories cannot be probed within the Milky Way if the screening threshold puts it in the screened regime. This is the case for most viable models. Nevertheless, it is possible to construct laboratory setups, in particular by means of a vacuum chamber with thick walls, in which the chameleon field is decoupled from the field outside [106]. This allows the creation of regions where the field remains unscreened, enabling tests of the mechanism by means of Casimir forces [107], levitated microspheres [108] and atom [109] and neutron [110] interferometry. All of these tests probe short-range fifth forces; see [111,112] for a review.

A number of other tests are possible at cosmological scales (see the review of [99]). In general, a cosmological fifth force causes an enhancement of the growth rate of structure at late times, which affects the CMB through the integrated Sachs-Wolfe and Sunyaev-Zel'dovich effects as well as gravitational lensing [113,114]. Constraints are also possible by means of the galaxy power spectrum [115];

comparison of dynamical and lensing masses of elliptical galaxies [116]; redshift space distortions [117]; and the abundance [118–120], internal structure [121] and relative dynamical and lensing masses [122,123] of clusters. These constraints are typically couched in terms of Hu-Sawicki models of $f(R)$, where $\Delta G/G_N = 1/3$ and $f_{R0}$ is related to $\lambda_c$ and $\Phi_c$ through Eq. (5). The strongest constraints are at the $f_{R0} \sim 10^{-5}$ level, showing cosmological tests to be significantly less useful for these theories than astrophysical ones.

## D. Suggestions for further work

Our work provides a case study of the use of the gravitational maps of [40] to constrain aspects of modified gravity on the scale of galaxies and their environments. This scale is characterized by $1 \text{ kpc} \lesssim r \lesssim 50 \text{ Mpc}$, $10^{-8} \lesssim |\Phi| \lesssim 10^{-4}$, $10^{-17} \lesssim |\vec{a}|/\text{km s}^{-2} \lesssim 10^{-13}$ and curvature $10^{-57} \lesssim K/\text{cm}^{-2} \lesssim 10^{-50}$.

Even continuing to focus solely on fifth force phenomenology, a range of intragalaxy statistics besides the separation of stars and gas may be expected to provide new information: as mentioned in Secs. I and VI C 1, a screened fifth force would induce warps in stellar disks and asymmetries in kinematics, as well as boost the rotation velocity of unscreened mass relative to screened mass [29]. Provided these signals could be suitably quantified, each could be used to derive constraints on $\Delta G$ and $\lambda_c$ by a method directly analogous to our own. As these analyses would utilize independent data and possess sensitivity to different types of galaxies in different environments, their results would be complementary. For the evidence for a chameleon- or symmetron-screened fifth force with $\lambda_c \simeq 2 \text{ Mpc}$ presented here to be convincing, it must be corroborated (or shown to be below the sensitivity level of current data) by each of these signals. We do this for galaxy warps in [79]. It is also possible to make predictions for future or present data sets from our best-fit $\{\lambda_c, \Delta G/G_N\}$ values, either by taking the maximum-likelihood point or by means of the Bayesian posterior predictive distribution. This would help to avoid confirmation bias in the analysis of that data. As discussed in Sec. VI B, the non-fifth-force part of the likelihood function could be refined by examining the signals in cosmological hydrodynamical $\Lambda$CDM simulations and the inferences made fully self-consistent by basing the determination of galaxy and large scale structure formation on simulations in modified gravity.

The constraints derivable from all of these signals will straightforwardly strengthen as future galaxy surveys begin observation. Given the amount of information evidently present in these data (as quantified by the strength of our constraints) and the great range of theories testable in tandem through simple phenomenological parametrizations, we propose that gravitational physics be incorporated as a key science driver of these surveys (see also [14,124,125]). Constraints based on intragalaxy
statistics may in fact scale better with survey quality than conventional tests based on intergalaxy statistics due to lower sensitivity to cosmic variance. At the very least they are complementary and come for little extra effort or expense with ongoing survey programs.

As discussed in [40], our methods may find more general use beyond the search for fifth forces. Any theory predicting deviations from GR that correlate with the local gravitational field may in principle be tested; this includes not only the other screening mechanisms listed in Sec. I (kinetic screening roughly occurs at a critical value of acceleration and Vainshtein at a critical value of curvature [23]), but also generic equivalence principle violations [e.g., in modified Newtonian dynamics (MOND), whose external field effect imposes an acceleration threshold] and models of dynamical dark energy which entail novel behavior near the curvature scale of $\Lambda$. Indeed, it is even possible to perform model-independent tests by simply correlating galaxies' gravitational properties with potential modified gravity signals, although the lack of a prediction would preclude forward modeling and Bayesian inference in this case. This method relies on the assumption that any breakdown of GR would impose a different dependence of dynamical behavior on gravitational environment than would contaminating “galaxy formation” effects, as we believe to be the case for the present signal $\vec{r}_*$. To check this assumption a more precise determination of the impact of baryonic physics would be helpful. This program would build gravity phenomenology from the ground up in galaxies, rather than test specific theories top down.

We made no *a priori* distinction here between galaxies in different environments, but rather fed the entire Alfalfa data set into our inference framework. This approach is ideal when all galaxies provide useful information on the parameters to be inferred: in our case, even massive galaxies or those in dense environments are of interest because their differences with the remainder inform the constraint on the screening threshold $\Phi_c$. The screened subsample for a given $\Phi_c$ calibrates the behavior of the signal in the absence of fifth forces. By contrast, information on some modified gravity theories may only be gleaned from galaxies in very specific environments. For example, testing the external field effect in equivalence-principle-violating theories would require galaxies in a dominant external field. The prototypical model in this class, MOND, requires $a_{ex} > a_0 \approx 1.2 \times 10^{-10} \text{ m/s}^2$ for the external field to dominate, which is at the extreme high end of the $a_{ex}$ distribution of 2M++ galaxies in the local Universe [40]. In such cases it will be preferable to begin by searching for regions of the local Universe satisfying certain constraints on the gravitational field and to analyze only the subset of galaxies situated therein.

To generalize the discussion in Sec. V D, we close by listing generic desiderata for surveys and modeling to be useful for probing fundamental physics with galaxies.

(i) Any test that depends on a galaxy's gravitational environment must make use of a reconstruction of the gravitational field along the lines of [40]. Improving the precision of the $\Phi$, $\vec{a}$ and $K$ maps requires principally an increase in the limiting magnitude out to which galaxies, the primary mass tracers, are visible. This is particularly important for pushing studies to higher redshift $z \gtrsim 0.05$ where current all-sky surveys become substantially incomplete. This will improve statistics and enhance the range of environments able to be probed. We estimate in [40] the improvement in the precision of gravitational parameter inferences afforded by upcoming photometric surveys. This information may be augmented by weak lensing data, spectroscopic information and photometry at other wavelengths to further constrain the overall distribution of mass.

(ii) Robust detections of the fifth force signals described above require low-redshift, high resolution data at multiple wavelengths to distinguish galaxies' different mass components. In general, data quality is more important than sample size, area or redshift range. There is no requirement that data be homogeneous or uniform across the sky, since spatial statistics of the galaxy population as a whole are not in question. This enables smaller surveys and isolated observations to play a larger role than is possible for traditional probes.

(iii) Resolved spectroscopy is necessary to study dynamical modified gravity signals such as differences in the rotation curves of stars and gas, kinematical asymmetries and mass discrepancies. These are complicated in the case of screened theories by the fact that the most common tracer of stellar kinematics, $H\alpha$ emission, derives from the ionized Strömgren spheres surrounding stars which are likely to be at least partly unscreened. The screened kinematics of the stars themselves must therefore be derived from stellar absorption lines, which are considerably more difficult to measure. Spectroscopy will also be useful to tighten constraints on galaxies' central dynamical masses, which provide the restoring forces on stellarcentroids that compensate for their insensitivity to $F_5$. The combination $\vec{a}_5 \Delta G(M(<r)/r^2)^{-1}$ appears in most of the formulas for screening signals [29].

(iv) When observational uncertainties are small, inference is limited by the precision with which the signals may be theoretically modeled on a galaxy-by-galaxy basis. This determines the width of the fifth force part of the likelihood function and is set primarily by $M(<r)$ and the gravitational parameters $\vec{a}_5$ and $\Phi$. Our calculation of the latter depends on the BORG inference of the smooth density field, and the precision of that machinery therefore feeds directly into the strength of our results.
Improvements are ongoing in the form of the BORG likelihood function, method for density field reconstruction and input data used. Our inference of $\bar{a}_5$ and $\Phi$ depends also on the galaxy-halo connection used to assign halos and hence total dynamical mass distributions to galaxies based on their photometry. In the future, additional observables will be brought to bear on constraining halo properties and both the scatter and systematic uncertainty in the galaxy-halo connection will fall.

## VII. CONCLUSION

We use the observed displacements between galaxies' stellar and gas mass centroids in the $D < 100$ Mpc part of the complete Alfalfa catalog to constrain fifth forces of range 0.4–50 Mpc that couple differentially to stars, gas and dark matter. Our primary case study is screened modified gravity, where the effect of a new scalar field is hidden in high-density environments by means of the chameleon or symmetron mechanism. In this case, stars self-screen, and hence decouple from the fifth force, in otherwise unscreened galaxies. We deploy the gravitational maps of [40] to identify screened and unscreened regions of the $z < 0.05$ Universe for a given fifth force range $\lambda_C$ and screening threshold $\Phi_c$, and we forward model the spatial offset between the optical and HI emission of each Alfalfa galaxy as a function of $\lambda_C$, $\Phi_c$ and fifth force strength $\Delta G$. We construct a Bayesian likelihood framework to constrain these parameters from the data, propagating fully the uncertainties in the gravitational maps, the galaxy structure inputs to the signal prediction, and the observations.

Models with $\Delta G/G_N \simeq 1$ and $\{\lambda_C, |\Phi_C|\}$ in the range $\{0.4–50$ Mpc, $2 \times 10^{-8} – 4 \times 10^{-4}\}$ produce HI-optical offsets of typical magnitude $\sim 0.1–1$ kpc for galaxies in the gravitational environments of the local Universe. This is slightly larger when screening is switched off so that all (rather than only unscreened) mass within $\lambda_C$ contributes to the fifth force. We find weak but generally positive correlation between the direction of the observed optical-HI displacement on the plane of the sky and that of an external screened fifth force field and between the magnitude of that displacement and the prediction of screened models. We quantify the statistical significance of these correlations under two main assumptions for the noise in the Alfalfa data. The fiducial model, which sets the noise Gaussian and equal on average to the scatter in measured optical-HI displacement in bins of the signal-to-noise ratio of the detection, yields $6.6\sigma$ evidence for a component of the offsets deriving from a screened fifth force, with best-fit solution $\lambda_C = 1.8$ Mpc, $\Delta G/G_N = 0.025$. Similar fifth force models without screening achieve no increase in likelihood over the general relativistic case $\Delta G = 0$. We validate this detection by a number of means but cannot prove that it may not be caused by other effects such

as galaxy formation physics. We will study such possible degeneracies further in future work.

A more conservative (and lower likelihood) model doubles the magnitude of the noise, enabling us to place the robust constraints $\Delta G/G_N \lesssim 0.1$ ($1\sigma$) for $\lambda_C = 0.5$ Mpc, and $\Delta G/G_N \lesssim$ few $\times 10^{-4}$ for $\lambda_C = 50$ Mpc in the case with screening, and $\Delta G/G_N \lesssim$ few $\times 10^{-3}$ and $\Delta G/G_N \lesssim$ few $\times 10^{-4}$ respectively in the case without. In $f(R)$ gravity (where $\Delta G/G_N = 1/3$) this corresponds to $f_{R0} \lesssim 2 \times 10^{-8}$. Even given this artificial inflation of the uncertainties, these constraints are 2 orders of magnitude stronger than those from cosmology, 1 order of magnitude stronger than those from distance indicators and rotation curve tests on galaxy scales, and approaching the strength of Solar System fifth force tests. This attests to the great constraining power of our data and methodology.

Our study is the first to leverage “big data” from galaxy surveys for gravitational physics by means of an intragalaxy signal. Such data may be expected to provide crucial information on gravity in the fully nonlinear regime in the future, as ever more powerful galaxy surveys are brought to bear. Besides probing qualitatively different regions of gravitational parameter space, intragalaxy tests should scale in strength more steeply with survey performance than conventional probes due to reduced sensitivity to sample variance and maximal utilization of information at low redshift where precision is greatest. Under optimistic assumptions for the noise characteristics of the HI data, applying our methodology to radio detections from SKA-class observatories would provide sensitivity to $\Delta G/G_N \sim \mathcal{O}(10^{-9})$ in screened theories with $\lambda_C \simeq 5$ Mpc. Finally, we discuss prospects for a range of additional intragalactic signals to yield further and complementary information on fundamental physics in this underexplored regime. Crucially, these analyses will confirm or rule out the possible screened fifth force detection achieved here.

## ACKNOWLEDGMENTS

This work was supported by St. John's College, Oxford. P. G. F. acknowledges support from Leverhulme, Science and Technology Facilities Council (STFC), Beecroft Institute for Particle Astrophysics and Cosmology (BIPAC) and the European Research Council (ERC). G. L. acknowledges support by the L'Agence nationale de la recherche (ANR) Grant No. ANR-16-CE23-0002 and from the Labex Institut Lagrange de Paris (ILP) (Reference No. ANR-10-LABX-63) part of the Idex Sorbonne Universités à Paris pour l'Enseignement et la Recherche (SUPER) (ANR-11-IDEX-0004-02). We thank Martha Haynes for sharing the complete Alfalfa catalog before public release, Phil Bull for information on SKA and Referee A of an earlier paper for helpful suggestions on validating the detection. Computations were performed at Oxford and SLAC.
[1] S. Weinberg, Phys. Rev. Lett. **19**, 1264 (1967).

[2] B. Ratra and P. J. E. Peebles, Phys. Rev. D **37**, 3406 (1988).

[3] C. Brans and R. H. Dicke, Phys. Rev. **124**, 925 (1961).

[4] C. Wetterich, Nucl. Phys. **B302**, 668 (1988).

[5] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. **513**, 1 (2012).

[6] E. G. Adelberger, B. R. Heckel, and A. E. Nelson, Annu. Rev. Nucl. Part. Sci. **53**, 77 (2003).

[7] E. G. Adelberger, B. R. Heckel, C. W. Stubbs, and W. F. Rogers, Annu. Rev. Nucl. Part. Sci. **41**, 269 (1991).

[8] C. D. Hoyle, U. Schmidt, B. R. Heckel, E. G. Adelberger, J. H. Gundlach, D. J. Kapner, and H. E. Swanson, Phys. Rev. Lett. **86**, 1418 (2001).

[9] J. C. Long, H. W. Chan, and J. C. Price, Nucl. Phys. **B539**, 23 (1999).

[10] J. Chiaverini, S. J. Smullin, A. A. Geraci, D. M. Weld, and A. Kapitulnik, Phys. Rev. Lett. **90**, 151101 (2003).

[11] E. Samain et al., Astron. Astrophys. Suppl. Ser. **130**, 235 (1998).

[12] J. G. Williams, X. X. Newhall, and J. O. Dickey, Phys. Rev. D **53**, 6730 (1996).

[13] B. Bertotti, L. Iess, and P. Tortora, Nature (London) **425**, 374 (2003).

[14] B. Jain and J. Khoury, Ann. Phys. (Amsterdam) **325**, 1479 (2010).

[15] J. Khoury, arXiv:1011.5909.

[16] J. Khoury and A. Weltman, Phys. Rev. D **69**, 044026 (2004).

[17] E. Babichev, C. Deffayet, and R. Ziour, Int. J. Mod. Phys. D **18**, 2147 (2009).

[18] A. I. Vainshtein, Phys. Lett. **39B**, 393 (1972).

[19] K. Hinterbichler and J. Khoury, Phys. Rev. Lett. **104**, 231301 (2010).

[20] S. M. Carroll, V. Duvvuri, M. Trodden, and M. S. Turner, Phys. Rev. D **70**, 043528 (2004).

[21] P. Brax, C. van de Bruck, A.-C. Davis, and D. Shaw, Phys. Rev. D **82**, 063519 (2010).

[22] P. Brax, A.-C. Davis, B. Li, and H. A. Winther, Phys. Rev. D **86**, 044015 (2012).

[23] J. Khoury, arXiv:1312.2006.

[24] A. Cabré, V. Vikram, G.-B. Zhao, B. Jain, and K. Koyama, J. Cosmol. Astropart. Phys. **7** (2012) 034.

[25] B. Falck, K. Koyama, and G.-B. Zhao, J. Cosmol. Astropart. Phys. **7** (2015) 049.

[26] G.-B. Zhao, B. Li, and K. Koyama, Phys. Rev. Lett. **107**, 071303 (2011).

[27] G.-B. Zhao, B. Li, and K. Koyama, Phys. Rev. D **83**, 044007 (2011).

[28] L. Hui, A. Nicolis, and C. W. Stubbs, Phys. Rev. D **80**, 104002 (2009).

[29] B. Jain and J. VanderPlas, J. Cosmol. Astropart. Phys. **10** (2011) 032.

[30] R. Giovanelli et al., Astron. J. **130**, 2598 (2005).

[31] R. Giovanelli et al., Astron. J. **133**, 2569 (2007).

[32] M. P. Haynes et al., Astron. J. **142**, 170 (2011).

[33] B. R. Kent, R. Giovanelli, M. P. Haynes, A. M. Martin, A. Saintonge, S. Stierwalt, T. J. Balonek, N. Brosch, and R. A. Koopmann, Astron. J. **136**, 713 (2008).

[34] A. M. Martin, R. Giovanelli, M. P. Haynes, A. Saintonge, G. L. Hoffman, B. R. Kent, and S. Stierwalt, Astrophys. J. Suppl. Ser. **183**, 214 (2009).

[35] A. Saintonge, R. Giovanelli, M. P. Haynes, G. L. Hoffman, B. R. Kent, A. M. Martin, S. Stierwalt, and N. Brosch, Astron. J. **135**, 588 (2008).

[36] V. Vikram, A. Cabré, B. Jain, and J. T. VanderPlas, J. Cosmol. Astropart. Phys. **8** (2013) 020.

[37] W. Hu and I. Sawicki, Phys. Rev. D **76**, 064004 (2007).

[38] A. De Felice and S. Tsujikawa, Living Rev. Relativity **13**, 3 (2010).

[39] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. **82**, 451 (2010).

[40] H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, Mon. Not. R. Astron. Soc. **474**, 3152 (2018).

[41] G. Lavaux and M. J. Hudson, Mon. Not. R. Astron. Soc. **416**, 2840 (2011).

[42] C. Conroy, R. H. Wechsler, and A. V. Kravtsov, Astrophys. J. **647**, 201 (2006).

[43] A. V. Kravtsov, A. A. Berlind, R. H. Wechsler, A. A. Klypin, S. Gottlober, B. Allgood, and J. R. Primack, Astrophys. J. **609**, 35 (2004).

[44] R. M. Reddick, R. H. Wechsler, J. L. Tinker, and P. S. Behroozi, Astrophys. J. **771**, 30 (2013).

[45] B. V. Lehmann, Y.-Y. Mao, M. R. Becker, S. W. Skillman, and R. H. Wechsler, Astrophys. J. **834**, 37 (2017).

[46] P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, Astrophys. J. **762**, 109 (2013).

[47] S. W. Skillman, M. S. Warren, M. J. Turk, R. H. Wechsler, D. E. Holz, and P. M. Sutter, arXiv:1407.2600.

[48] B. Diemer and A. V. Kravtsov, Astrophys. J. **799**, 108 (2015).

[49] J. Jasche, F. S. Kitaura, B. D. Wandelt, and T. A. Enßlin, Mon. Not. R. Astron. Soc. **406**, 60 (2010).

[50] J. Jasche and B. D. Wandelt, Mon. Not. R. Astron. Soc. **425**, 1042 (2012).

[51] J. Jasche and B. D. Wandelt, Mon. Not. R. Astron. Soc. **432**, 894 (2013).

[52] J. Jasche, F. Leclercq, and B. D. Wandelt, J. Cosmol. Astropart. Phys. **1** (2015) 036.

[53] G. Lavaux and J. Jasche, Mon. Not. R. Astron. Soc. **455**, 3169 (2016).

[54] J. Jasche and G. Lavaux, arXiv:1806.11117.

[55] B. M. B. Henriques, S. D. M. White, P. A. Thomas, R. Angulo, Q. Guo, G. Lemson, V. Springel, and R. Overzier, Mon. Not. R. Astron. Soc. **451**, 2663 (2015).

[56] A. A. Dutton et al., Mon. Not. R. Astron. Soc. **416**, 322 (2011).

[57] F. Lelli, Galaxies **2**, 292 (2014).

[58] F. Lelli, S. S. McGaugh, J. M. Schombert, and M. S. Pawlowski, Astrophys. J. Lett. **827**, L19 (2016).

[59] M. Milgrom, Phys. Rev. Lett. **117**, 141101 (2016).

[60] B. Famaey and S. S. McGaugh, Living Rev. Relativity **15**, 10 (2012).

[61] J. Janz, M. Cappellari, A. J. Romanowsky, L. Ciotti, A. Alabi, and D. Forbes, Mon. Not. R. Astron. Soc. **461**, 2367 (2016).

[62] S. S. McGaugh, Astrophys. J. **609**, 652 (2004).

[63] R. H. Sanders, Astron. Astrophys. Rev. **2**, 1 (1990).
[64] W. J. G. de Blok, *Adv. Astron.* **2010**, *1* (2010).

[65] D. C. Rodrigues, A. del Popolo, V. Marra, and P. L. C. de Oliveira, Mon. Not. R. Astron. Soc. **470**, 2410 (2017).

[66] M. R. Blanton, E. Kazin, D. Muna, B. A. Weaver, and A. Price-Whelan, Astron. J. **142**, 31 (2011).

[67] M. Bernardi *et al.*, Astron J. **125**, 1817 (2003).

[68] M. R. Blanton and S. Roweis, Astron. J. **133**, 734 (2007).

[69] G. Chabrier, Publ. Astron. Soc. Pac. **115**, 763 (2003).

[70] M. P. Haynes and R. Giovanelli, Astron. J. **89**, 758 (1984).

[71] J. L. Sérsic, *Atlas de Galaxias Australes* (Observatorio Astronomico, Cordoba, Argentina, 1968).

[72] L. Ciotti and G. Bertin, Astron. Astrophys. **352**, 447 (1999).

[73] P. Prugniel and F. Simien, Astron. Astrophys. **321**, 111 (1997).

[74] R. B. Tully and P. Fouque, Astrophys. J. Suppl. Ser. **58**, 67 (1985).

[75] A. Begum, J. N. Chengalur, I. D. Karachentsev, S. S. Kaisin, and M. E. Sharina, Mon. Not. R. Astron. Soc. **365**, 1220 (2006).

[76] M. Geha, M. R. Blanton, M. Masjedi, and A. A. West, Astrophys. J. **653**, 240 (2006).

[77] H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, arXiv:1802.07206.

[78] J. Alsing, E. Berti, C. M. Will, and H. Zaglauer, Phys. Rev. D **85**, 064041 (2012).

[79] H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, arXiv:1807.11742.

[80] M. Santos, D. Alonso, P. Bull, M. B. Silva, and S. Yahya, *Proc. Sci. AASKA14* (2014) 21.

[81] S. Yahya, P. Bull, M. G. Santos, M. Silva, R. Maartens, P. Okouma, and B. Bassett, Mon. Not. R. Astron. Soc. **450**, 2251 (2015).

[82] J. Sakstein, Phys. Rev. D **97**, 064028 (2018).

[83] M. A. W. Verheijen, T. A. Oosterloo, W. A. van Cappellen, L. Bakker, M. V. Ivashina, and J. M. van der Hulst, The Evolution of Galaxies Through the Neutral Hydrogen Window, **1035**, 265 (2008).

[84] D. R. DeBoer et al., *IEEE Proc. 97*, 1507 (2009).

[85] S. Kaviraj et al., Mon. Not. R. Astron. Soc. **467**, 4739 (2017).

[86] N. Khandai, T. Di Matteo, R. Croft, S. Wilkins, Y. Feng, E. Tucker, C. DeGraf, and M. Liu, Mon. Not. R. Astron. Soc. **450**, 1349 (2015).

[87] J. Schaye et al., Mon. Not. R. Astron. Soc. **446**, 521 (2015).

[88] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu, G. Snyder, D. Nelson, and L. Hernquist, Mon. Not. R. Astron. Soc. **444**, 1518 (2014).

[89] T. K. Chan, D. Kereš, J. Oñorbe, P. F. Hopkins, A. L. Muratov, C.-A. Faucher-Giguère, and E. Quataert, Mon. Not. R. Astron. Soc. **454**, 2981 (2015).

[90] P. F. Hopkins, D. Kereš, J. Oñorbe, C.-A. Faucher-Giguère, E. Quataert, N. Murray, and J. S. Bullock, Mon. Not. R. Astron. Soc. **445**, 581 (2014).

[91] T. Sawala et al., Mon. Not. R. Astron. Soc. **457**, 1931 (2016).

[92] G. S. Stinson, C. Brook, A. V. Macciò, J. Wadsley, T. R. Quinn, and H. M. P. Couchman, Mon. Not. R. Astron. Soc. **428**, 129 (2013).

[93] N. E. Chisari et al., arXiv:1801.08559.

[94] W. A. Hellwing, M. Schaller, C. S. Frenk, T. Theuns, J. Schaye, R. G. Bower, and R. A. Crain, Mon. Not. R. Astron. Soc. **461**, L11 (2016).

[95] M. Schaller, C. S. Frenk, R. G. Bower, T. Theuns, A. Jenkins, J. Schaye, R. A. Crain, M. Furlong, C. Dalla Vecchia, and I. G. McCarthy, Mon. Not. R. Astron. Soc. **451**, 1247 (2015).

[96] A. J. Mead, C. Heymans, L. Lombriser, J. A. Peacock, O. I. Steele, and H. A. Winther, Mon. Not. R. Astron. Soc. **459**, 1468 (2016).

[97] E. Puchwein, M. Baldi, and V. Springel, Mon. Not. R. Astron. Soc. **436**, 348 (2013).

[98] F. Fontanot, E. Puchwein, V. Springel, and D. Bianchi, Mon. Not. R. Astron. Soc. **436**, 2672 (2013).

[99] L. Lombriser, Ann. Phys. (Amsterdam) **526**, 259 (2014).

[100] L. Lombriser, B. Li, K. Koyama, and G.-B. Zhao, Phys. Rev. D **87**, 123511 (2013).

[101] D. Shi, B. Li, J. Han, L. Gao, and W. A. Hellwing, Mon. Not. R. Astron. Soc. **452**, 3179 (2015).

[102] A.-C. Davis, E. A. Lim, J. Sakstein, and D. J. Shaw, Phys. Rev. D **85**, 123006 (2012).

[103] V. Vikram, J. Sakstein, C. Davis, and A. Neil, Phys. Rev. D **97**, 104055 (2018).

[104] B. Jain, V. Vikram, and J. Sakstein, Astrophys. J. **779**, 39 (2013).

[105] T. Baker, D. Psaltis, and C. Skordis, Astrophys. J. **802**, 63 (2015).

[106] C. Burrage, E. J. Copeland, and E. A. Hinds, J. Cosmol. Astropart. Phys. **3** (2015) 042.

[107] R. S. Decca, D. López, E. Fischbach, G. L. Klimchitskaya, D. E. Krause, and V. M. Mostepanenko, Phys. Rev. D **75**, 077101 (2007).

[108] A. D. Rider, D. C. Moore, L. M. Blakemore, M. Lu, and G. Gratta, Phys. Rev. Lett. **117**, 101101 (2016).

[109] P. Hamilton, M. Jaffe, P. Haslinger, Q. Simmons, H. Müller, and J. Khoury, Science **349**, 849 (2015).

[110] P. Brax, G. Pignol, and D. Roulier, Phys. Rev. D **88**, 083004 (2013).

[111] C. Burrage and J. Sakstein, J. Cosmol. Astropart. Phys. **11** (2016) 045.

[112] C. Burrage and J. Sakstein, Living Rev. Relativity **21**, 1 (2018).

[113] L. Lombriser, A. Slosar, U. Seljak, and W. Hu, Phys. Rev. D **85**, 124038 (2012).

[114] Y.-S. Song, H. Peiris, and W. Hu, Phys. Rev. D **76**, 063517 (2007).

[115] J. Dossett, B. Hu, and D. Parkinson, J. Cosmol. Astropart. Phys. **3** (2014) 046.

[116] T. L. Smith, arXiv:0907.4829.

[117] K. Yamamoto, G. Nakamura, G. Hütsi, T. Narikawa, and T. Sato, Phys. Rev. D **81**, 103517 (2010).

[118] M. Cataneo, D. Rapetti, F. Schmidt, A. B. Mantz, S. W. Allen, D. E. Applegate, P. L. Kelly, A. von der Linden, and R. G. Morris, Phys. Rev. D **92**, 044009 (2015).
[119] S. Ferraro, F. Schmidt, and W. Hu, Phys. Rev. D **83**, 063503 (2011).

[120] F. Schmidt, A. Vikhlinin, and W. Hu, Phys. Rev. D **80**, 083505 (2009).

[121] L. Lombriser, F. Schmidt, T. Baldauf, R. Mandelbaum, U. Seljak, and R. E. Smith, Phys. Rev. D **85**, 102001 (2012).

[122] A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon, K. Koyama, and R. C. Nichol, J. Cosmol. Astropart. Phys. **4** (2014) 013.

[123] H. Wilcox et al., Mon. Not. R. Astron. Soc. **452**, 1171 (2015).

[124] B. Jain, Phil. Trans. R. Soc. A **369**, 5081 (2011).

[125] B. Jain et al., arXiv:1309.5389.