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

HAL
open science

# Fifth force constraints from galaxy warps

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 galaxy warps. Physical Review D, 2018, 98 (8), pp.083010. 10.1103/PhysRevD.98.083010. hal-01867485

HAL Id: hal-01867485

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

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 galaxy warps

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 bd Arago, 75014 Paris, France*
*and Sorbonne Universités, Institut Lagrange de Paris (ILP), 98 bis bd 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 31 July 2018; published 15 October 2018)

Intragalaxy signals contain a wealth of information on fundamental physics, both the dark sector and the nature of gravity. While so far largely unexplored, such probes are set to rise dramatically in importance as upcoming surveys provide data of unprecedented quantity and quality on galaxy structure and dynamics. In this paper, we use warping of stellar disks to test the chameleon- or symmetron-screened fifth forces which generically arise when new fields couple to matter. We take r-band images of mostly late-type galaxies from the *Nasa Sloan Atlas* and develop an automated algorithm to quantify the degree of U-shaped warping they exhibit. We then forward model the warp signal as a function of fifth-force strength, $\Delta G/G_N$, and range, $\lambda_C$, and the gravitational environments and internal properties of the galaxies, including full propagation of the non-Gaussian uncertainties. Convolving this fifth-force likelihood function with a Gaussian describing astrophysical and observational noise and then constraining $\Delta G/G_N$ and $\lambda_C$ by Markov chain Monte Carlo, we find the overall likelihood to be significantly increased ($\Delta \log(L) \approx 20$) by adding a screened fifth force with $\lambda_C \approx 2$ Mpc and $\Delta G/G_N \approx 0.01$. The variation of $\Delta \log(L)$ with $\lambda_C$ is quantitatively as expected from the correlation of the magnitude of the fifth-force field with the force's range, and a similar model without screening achieves no increase in likelihood over the General Relativistic case $\Delta G = 0$. Although these results are in good agreement with a previous analysis of the same model using offsets between galaxies' stellar and gas mass centroids [H. Desmond et al., Phys. Rev. D **98**, 064015 (2018).], we caution that the effects of confounding baryonic and dark matter physics must be thoroughly investigated for the results of the inference to be unambiguous.

DOI: 10.1103/PhysRevD.98.083010

## I. INTRODUCTION

Our current theory of fundamental physics is both successful and puzzling. We are able to explain almost all phenomena at high energies and on cosmological scales in terms of a subset of gauge theories of particles and fields combined with the theory of General Relativity. Yet, while the agreement with data is remarkable, a number of conceptual and practical problems have emerged: why is there a hierarchy of scales, e.g., and what is the nature of the Universe's dark sector?

At the heart of almost all attempts to solve these problems are extra degrees of freedom, or new fundamental fields. These emerge when effective field theories break down and when symmetries or conservation laws are broken. The Higgs field is a notable example within the standard model, as is the attempt outside of it to understand inertia by endowing the gravitational constant with dynamics [1]. They also arise naturally in effective descriptions of cosmology from UV-complete theories. The ubiquity of extra fields in attempts to extend the standard model makes them the smoking gun for new physics.

New fields generically result in new (fifth) forces. Furthermore, through quantum consistency at low energy, these fields will couple nonminimally to the gravitational

<footnote>*harry.desmond@physics.ox.ac.uk</footnote>
sector, making the forces gravitational [2]. A notable example is a scalar field $\varphi$ with a potential $V(\varphi) \propto m^2\varphi^2 + \lambda\varphi^4$; quantum corrections will inevitably lead to a nonminimal coupling of the form $\xi\varphi^2 R$ (where $R$ is the Ricci scalar) with $\xi \propto \lambda$ [3].

Gravitational fifth forces can be interpreted as arising from corrections to the Newtonian potential $\Phi$. Near a localized body of mass $M$,

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

where $G_N$ is Newton's constant and $\Delta G$ and $m$ parametrize the relative strength and inverse range of the fifth force. If $m$ is large, the force is short range (i.e., small Compton wavelength $\lambda_C$), but if $m \approx 0$, the fifth force may compete with gravity on all scales. Stringent constraints have been placed on $\Delta G/G_N$ and $m$ in a variety of regimes, to the point that universally coupled fifth forces are effectively irrelevant on both laboratory and astrophysical scales [4].

Over the last decade, however, a better understanding of fifth forces in a cosmological setting has emerged. It has been shown that many theories are equipped with a gravitational screening mechanism, such that either the range of the fifth force is vastly reduced or its strength becomes negligible in particular environments. For example, in the case of chameleon [5] or symmetron [6] screening, the fifth forces switch off in regions of high density, as quantified by a deep Newtonian potential; in the case of Vainshtein screening [7], fifth forces are suppressed close to sufficiently massive bodies. Fifth forces may therefore be undetectable in the Solar System and laboratory while playing an important role cosmologically or in regions of space with weak gravitational field.

Screening necessitates that one searches for astrophysical fifth forces in special places, i.e., in particular cosmic environments. The first step, then, is to construct a reliable gravitational map of a volume of interest. By taking into account the various contributions sourcing the gravitational field, it is possible to demarcate regions of the local Universe by Newtonian potential, acceleration, or curvature. In Ref. [8], we presented state-of-the-art maps of this type out to $z = 0.05$, built from the 2M++ all-sky redshift survey [9] combined with the halo population of an N-body simulation [10] and constrained realizations of the long wavelength modes of the density field [11].

With the ability to identify “interesting” regions, one can now look for specific signatures of fifth forces. A particularly fertile, and mostly unexplored, regime is within galaxies, where large and growing data sets exist unhampered by cosmic variance. The trade-off is that potential degeneracies arise from nongravitational (dissipative) galactic physics as well as uncertainties in the location of the dark matter that governs galaxies’ dynamics and drives their evolution. Nevertheless, conservative models of

all possible contaminants combined with statistical samples of galaxies in a range of gravitational environments should allow tests of screened fifth forces with existing data more powerful than are possible either within the Solar System or in cosmology.

One such attempt looked for offsets between the center of emission of the stars (which typically screen themselves) and HI gas of individual galaxies (Refs. [12,13], hereafter D18). We demonstrated great sensitivity to $\Delta G/G_N$ and $m$: with strongly conservative measurement errors, we set $\Delta G/G_N < \text{few} \times 10^{-4}$ for fifth-force ranges $\sim 50$ Mpc, while a more realistic noise model led to a $6.6\sigma$ detection of $\Delta G/G_N \approx 0.025$ at a range $\lambda_C \approx 1.8$ Mpc. While this result is striking and somewhat intriguing, there may be unaccounted-for systematics that need to be better understood; in particular, the differential impact of baryonic physics on stars and gas may be significant, and a proper assessment of this effect awaits study in high-resolution cosmological hydrodynamical simulations of galaxy formation.

We focus here on another signal of fifth forces, the warping of stellar galactic disks. Most late-type galaxies appear to be warped to some degree, with the most prominent examples observed in HI. There is weaker evidence for warping in stellar disks, although statistically significant samples have recently begun to be compiled. The origin of warps is not fully understood; they are likely due largely to environmental factors such as tidal interactions or intergalactic magnetic fields, although the observation of warps in isolated galaxies requires that internal processes play some role. For a review of the observational and theoretical status of warps, see Ref. [14].

Another source of warping is a screened fifth force [15]. If the (screened) stellar disk lags behind the halo center in an overall unscreened galaxy, the differential force due to the dark matter across the disk establishes a bowl-shaped equilibrium profile for the stellar mass. This effect is more pronounced the stronger the fifth force is and the better its alignment with the disk normal is, and is more visible on the plane of the sky the nearer the galaxy is to being viewed edge on.

In this paper, we develop in detail the formalism for warping under a fifth force and create a likelihood framework with which the properties of the force may thereby be inferred. We then compile and analyze a catalog of warps in mostly late-type galaxies to constrain $\Delta G/G_N$ and $\lambda_C$. Our analysis parallels in many respects that of D18, although—and crucially—its different observational and theoretical inputs give it a largely disjoint set of systematics. Our inference is therefore complementary to D18 and, as we shall see, reinforces its conclusions.

The structure of the paper is as follows. In Sec. II, we model galaxy structure under a fifth force and develop a summary statistic to quantify the overall strength of a U-shaped warp. In Sec. III, we describe our observational
sample. In Sec. IV, we detail our method, including the measurement of warp strength in the observations, calculation of the predicted fifth-force signal, and comparison of the two with a Bayesian likelihood model. In Sec. V, we present constraints on $\Delta G/G_N$ and $\lambda_C$ and document a range of checks aimed at validation. Section VI discusses our results, with particular emphasis on potentially confounding systematics, complementarity to D18, and necessary further work. Section VII concludes.

## II. GALAXY WARPS UNDER A FIFTH FORCE

### A. Heuristic overview

The basic requirement for the type of warping we investigate is that a galaxy's stellar disk and dark matter halo experience different accelerations due to their common environment. This may be caused by an additional force besides standard gravity acting on one of the two components, either a long-range interaction with surrounding mass or an additional interaction between the particles comprising one of the two subsystems themselves. Effects of the latter type arise from dynamical friction and dark matter self-interactions; in this paper, we focus on the former.

Consider a thin stellar disk at the center of a dark matter halo, with both mass components accelerating together in an external field. Now, suppose that the halo is subject to an additional force which the stars do not feel, so that it experiences an extra acceleration $\vec{a}_5$. The halo center will move ahead of the stellar centroid until the acceleration of the latter due to its offset $r_*$ from the former exactly compensates for its not feeling the extra force. This sets up a potential gradient across the disk that causes it to warp into a U-shape; heuristically, as the total halo acceleration is lower at larger galactocentric radius, it must point at a smaller angle to $\vec{a}_5$ to equal the fixed acceleration difference between the two components (Fig. 1). We calculate

<figure><img src="image_4.png" /><figcaption>FIG. 1. Schematic illustration of the formation of a U-shaped warp from a separation of stellar disk and halo centers. $\vec{a}_5$ is the acceleration resulting from the fifth force, which in the case of thin-shell screening (e.g., chameleon and symmetron) is typically felt by the dark matter but not the stars.</figcaption></figure>

the expected warp curve precisely in Sec. II B, and in Sec. II C, we integrate over it with a suitable kernel to create a single measure of the overall U-shaped warp strength. We then calculate this quantity analytically for power-law halo density profiles.

The case in which we are particularly interested is a chameleon- or symmetron-screened fifth force [5,6], characterized by a strength relative to gravity, $\Delta G/G_N$, and a range $\lambda_C$ inversely related to the scalar field's mass. The self-screening parameter $\chi = \phi_0/(2\beta M_{pl})$—where $\phi_0$ is the background value of the scalar field, $\beta$ is its coupling coefficient to normal matter, and $M_{pl}$ is the Planck mass—sets the threshold $|\Phi_c|$ in the Newtonian potential above which an object is screened [16–18]. For values of $|\Phi_c|$ below $\sim 10^{-6}$ ($c \equiv 1$), main-sequence stars screen themselves even if the galaxy as a whole is unscreened [19]. Thus, the halo feels a fifth force due to surrounding unscreened mass, while the stellar disk does not, causing an effective equivalence principle violation of precisely the type necessary to generate the warping described above.

### B. Warp curve expected from a screened fifth force

In this section, we calculate the shape of the stellar disk resulting from a screened fifth force and investigate its dependence on both fifth-force and galaxy structural parameters. We begin by defining a cylindrical coordinate system with the $z$ axis along the disk normal and origin coincident with the halo center and consider a star moving in a circular orbit around the $z$ axis (Fig. 1). D18 (and previously [15]) considered the equilibrium separation $r_* \equiv z(x=0)$ of the stellar centroid from the halo center,

$$ \frac{M(<r_*)}{r_*^2} \hat{r}_* = \vec{a}_5 \frac{\Delta G}{G_N^2}, \qquad (2) $$

where $M(<r)$ is the enclosed dark matter plus gas mass. We now generalize this calculation to require equilibrium at an arbitrary point along the disk. The condition is

$$ a_{5,z} \frac{\Delta G}{G_N} = a_{h,z} = -a_h \frac{z}{r}, \qquad (3) $$

where $r = \sqrt{z^2 + x^2}$, $a_h$ is the magnitude of the acceleration due to the halo, and subscript “, z” denotes $z$ projection. Assuming the halo is spherically symmetric [$a_h = G_N M(<r)/r^2$], the disk must satisfy

$$ z(x) = -a_{5,z} \frac{\Delta G}{G_N^2} \frac{r^3}{M(<r)}. \qquad (4) $$

Thus, negative $a_{5,z}$, as in Fig. 1, implies positive $z$. The amplitude of this effect depends on the same function of the external fifth-force field, scalar coupling, and total density profile as the offset $\vec{r}_*$ in Eq. (2), viz. $a_5 \Delta G/M(<r)$, but is now a function of position along the disk. For realistic halo
density profiles which fall with $r$, $z$ is an increasing function of $x$ so that the disk acquires a U shape bending away from $\vec{a}_{5,z}$.

We now assume that the gravitational restoring force is dominated by the dark matter, so that $M \simeq M_{\text{halo}}$, and that the warp is small, so that $z(x) \ll x$. The first assumption is justified in the Appendix, and the second is justified postfacto from the warp strengths resulting from reasonable fifth-force models (Sec. VA). This simplifies the above equation to<fnref n="1" />

$$z(x) = -a_{5,z} \frac{\Delta G}{G_N^2} \frac{|x|^3}{M_{\text{halo}}(<x)}. \quad (5)$$

For given $a_5$, the strength of the warp is maximized when the disk falls face on in the external fifth-force field. Conversely, if the disk falls edge on, no warping in the direction of the disk normal is expected, although asymmetries will develop in the plane of the disk and in its kinematics [15].

### C. Summarizing the warp curve

To simplify our inference and minimize the impact of noise, we compress the warp curve into a scalar summary statistic. Our choice is motivated by the desire to maximize sensitivity to U-shaped warps, which are predicted by the fifth-force model and less commonly observed with large magnitude than S-shaped ones. Following Ref. [20], we choose

$$w_1 \equiv \frac{1}{L_x^3} \int_{-L_x}^{L_x} |x| z(x) dx, \quad (6)$$

where $L_x$ is the distance along the major axis from the center of the galaxy out to which the warp curve is calculated. The choice of $L_x$ is a trade-off between maximizing the information included from the outer regions of the disk where the warp is most pronounced, on the one hand, and reducing the contamination from low signal-to-noise regions where sky fluctuations dominate, on the other. Our fiducial choice is $L_x = 3R_{\text{eff}}$ (with $R_{\text{eff}}$ the half-light radius of the disk), although we check that qualitatively similar results follow from 2 or 4 $R_{\text{eff}}$. $w_1$ is unitless and zero for an S-shaped warp. The determination of $w_1$ in the observational data is described in Sec. IVA; here, we continue to focus on the model prediction.<fnref n="2" />

To perform Bayesian inference, we must forward model the signal, which requires solving Eqs. (5) and (6) to find the expected $w_1$ for given fifth-force and galaxy parameters. As galaxies occupy the central regions of halos (typically within the scale radius), we begin by approximating the dark matter density over the extent of the galaxy by a single power law in $r$.<fnref n="3" /> To allow for more general profiles than NFW (e.g., to take account of the effect of baryonic feedback), we keep the power-law index free and assign a specific value only when a numerical answer is required:

$$\rho(r) \simeq \rho(r_s)(r/r_s)^{-n}. \quad (7)$$

We pivot the power law at the halo scale radius $r_s$ because this is typically well measured in NFW fits to the halos produced in N-body simulations. These will be employed below to estimate halo properties, in conjunction with empirical models for the galaxy-halo connection. This gives an enclosed mass,

$$M(<R) = \frac{4\pi\rho_{rs}r_s^nR^{3-n}}{3-n}, \quad (8)$$

where $\rho_{rs} \equiv \rho(r_s)$. Assuming that the warp is small ($z \ll x$, justified in Sec. VA),

$$R = (x^2 + z^2)^{1/2} \approx x + O(z^2). \quad (9)$$

Plugging Eqs. (8) and (9) into Eq. (5), we find the warp curve

$$z(x) \approx K|x|^n, \quad (10)$$

where

$$K = - \frac{a_{5,z} \Delta G}{G_N^2} \frac{3-n}{4\pi\rho_{rs}r_s^n}. \quad (11)$$

We show example warp curves for three different values of $n$ in Fig. 2.

To facilitate comparison with the observations, we now define a new variable $z'$ with a mean of 0 across the disk out to $L_x \equiv 3R_{\text{eff}}$,

$$\langle z \rangle = \frac{1}{L_x} \int_0^{L_x} z(x) dx = \frac{K}{n+1} L_x^n, \quad (12)$$

so that

$$z'(x) \equiv z(x) - \langle z \rangle = K \left( |x|^n - \frac{L_x^n}{n+1} \right). \quad (13)$$

<footnote n="1">¹Note that in this approximation $z = 0$ at $x = 0$, so that Eq. (2) is not reproduced. An overall offset is not of interest here but only the variation of $z$ with $x$; in practice we will redefine the $z$ axis anyway so that the mean position of the warped disk is at $z' = 0$. (This also increases independence from D18.)</footnote>

<footnote n="2">²The compression of the warp curve into $w_1$ necessarily throws away information and may therefore be expected to weaken our inference. It is not clear, however, whether finer-grained features of the warp curve possess significant additional constraining power over the noise; we leave investigation of this issue to future work.</footnote>

<footnote n="3">³We have checked that $r_s > 3R_{\text{eff}}$ for the great majority of galaxies in most abundance matching realizations (see Sec. IV B).</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Expected warp curves deriving from a fifth force according to Eq. (10), as a function of the power-law slope n of the halo density profile within the stellar disk. To evaluate K, we take typical illustrative values $r_s = 10$ kpc, $\rho_{rs} = 10^7 M_\odot$ kpc$^{-3}$, $a_{5,z} = 10^{-16}$ km/s, and $\Delta G/G_N = 1$. The disk’s major axis lies along x, and z is the orthogonal direction on the plane of the sky. Fifth-force warping is an $\mathcal{O}$(10 pc) effect; sensitivity to it requires a statistical sample of galaxies in a range of gravitational environments.</figcaption></figure>

Finally, we evaluate the integral for $w_1$ in these transformed coordinates,

$$w_1 = \frac{2}{L_x^3} \int_0^{L_x} z'(x) x dx = \frac{2K}{L_x^3} \left( \frac{L_x^{n+2}}{n+2} - \frac{L_x^{n+2}}{2(n+1)} \right) \\ = K \frac{n}{(n+1)(n+2)} L_x^{n-1}. \qquad (14)$$

where we have used the symmetry of the warp about the midpoint of the disk to halve the range of the integral. Substituting for $L_x$ and K yields

$$w_1 = - \frac{n(3-n)}{(n+1)(n+2)} a_{5,z} \frac{\Delta G}{G_N^2} \frac{1}{4\pi\rho_{rs}} \frac{(3R_{\text{eff}}/r_s)^n}{3R_{\text{eff}}}. \qquad (15)$$

### III. OBSERVATIONAL DATA

We use the *Nasa Sloan Atlas* (NSA),<fnref n="4" /> a compilation of UV, optical, and near-IR data for galaxies within ~250 Mpc. To parallel D18, we will be interested in fifth-force ranges up to 50 Mpc, which, given that our screening maps extend only to ~200 Mpc [8], forces us to restrict our observational data set to $D \lesssim 100$ Mpc. This is because we need to know the screening properties not only of the test galaxies but also of the objects sourcing the fifth forces felt by those galaxies, as screened objects do not contribute to $\bar{a}_5$. We also cut at $M_* > 10^9 M_\odot$, necessary for precise characterization of the halos of the test galaxies as required for evaluating Eq. (15). We use the r-band images to measure the warp (Sec. IVA) and record also the galaxy stellar mass, $M_*$, the apparent minor-to-major axis

ratio, $b/a$, the angular half-light radius, SERSIC_TH50, and the angle of the disk’s major axis relative to north, SERSIC PHI, which we will use in modeling the warps (Sec. IV B). We show three example raw images in Fig. 3.

In Sec. II, we modeled test galaxies as thin disks oriented exactly edge on. Although it is possible to restrict our observational sample to such disks ($b/a = 0.15$), this would drastically reduce our statistics. A puffy disk that is oriented edge on may continue to be analyzed with the method of Sec. II B, provided that $z(x)$ is taken to be the center of mass across the slice of the disk at x and that tides across the galaxy are not significant. As most of our galaxies are isolated, this should be a reasonable approximation. For a galaxy not oriented fully edge on, it is more difficult to predict exactly how the warp should appear. As a compromise between statistics and possible systematic error, we cut at $b/a = 0.5$. Although we cannot gauge the uncertainty in the expected $w_1$ ensuing from a fifth force for galaxies of more complex geometry, we can to first order model the change in contribution from other physics by introducing a dependence on $b/a$ of the width of the non-fifth-force part of the likelihood function. This is described in Sec. IV C. With these cuts, our final sample contains 4,206 galaxies at a mean distance of 74 Mpc.

### IV. METHOD

Our overall method is similar to D18 in which further details may be found. We forward model $w_1$ as a function of the global parameters $\Delta G/G_N$ and $\lambda_C$ and the galaxy-specific parameters $R_{\text{eff}}$, $r_s$, $\rho_{rs}$, $\Phi$, and $\bar{a}_5$. We specify probability distributions for the galaxy-specific parameters to create a fifth-force likelihood function $\mathcal{L}_5(w_1|\Delta G/G_N, \lambda_C)$ for each NSA galaxy, convolve it with a Gaussian that models the effect of non-fifth-force physics, and constrain $\{\lambda_C, \Delta G/G_N\}$ by Markov chain Monte Carlo (MCMC). We begin with our measurements of $w_1$ from the NSA sample before detailing our forward model.

#### A. Measuring $w_1$

We begin by rotating each image through SERSIC PHI so that the galaxy’s major axis is along the direction we designate as x. Taking the perpendicular axis on the plane of the sky to be z, we set the origin of our coordinate system at the centroid of the galaxy’s r-band emission and use units of pixel size. We scan in the x direction through all pixels between $\pm 3R_{\text{eff}}$ and for each vertical slice calculate the intensity-weighted mean value of z, $\bar{z}(x)$, between $z = \pm L_z \equiv 3(b/a)R_{\text{eff}}$,

$$\bar{z}(x) = \frac{\sum_{z=-L_z}^{z=L_z} zI(x, z)}{\sum_{z=-L_z}^{z=L_z} I(x, z)}, \qquad (16)$$

where the sum is over pixels and $I(x, z)$ is the measured flux in pixel $(x, z)$. For the observational data, $\bar{z}$ replaces z

<footnote n="4">⁴http://www.nsatlas.org/.</footnote>
<figure><img src="image_6.png" /></figure>

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

<figure><img src="image_8.png" /><figcaption>FIG. 3. Example r-band images of galaxies from the *Nasa Sloan Atlas*, labeled by International Astronomical Union (IAU) designation and apparent axis ratio *b/a*.</figcaption></figure>

in Eqs. (6) and (15) as a measure of the average position of the disk at given *x*. The choice of $L_z$ is intended to roughly cover the extent of the disk along its minor axis; we have checked that our conclusions are qualitatively unaffected by changing $L_z$ within 2–4 (b/a)$R_{\text{eff}}$.

The warp curves $\bar{z}(x)$ of the galaxies in Fig. 3 are shown in Fig. 4. These show that there is some offset between the origin of the coordinate system and the center of the disk; to remove the bias this would induce in $w_1$, we redefine the *z* axis so that the center of the warp curve is at $z' = 0$:

$$ \langle z \rangle = \frac{1}{6R_{\text{eff}}} \sum_{x=-3R_{\text{eff}}}^{x=3R_{\text{eff}}} \bar{z}(x), \quad \bar{z}'(x) \equiv \bar{z}(x) - \langle z \rangle. \quad (17) $$

Finally, we evaluate Eq. (6) as

$$ w_1 = \frac{1}{(3R_{\text{eff}})^3} \sum_{x=-3R_{\text{eff}}}^{x=3R_{\text{eff}}} |x| \bar{z}'(x). \quad (18) $$

## B. Forward modeling $w_1$

As $w_1$ vanishes in overall-screened galaxies and is otherwise proportional to the fifth-force field $\vec{a}_5$, the first step to forward modeling it on a galaxy-by-galaxy basis is to precisely characterize galaxies’ gravitational environments. To do so, we utilize the maps of Ref. [8], where both $\Phi_{\text{ex}}$ and $\vec{a}_5$ are calculated as a function of fifth-force range $\lambda_C$ at every point in space to a distance ~200 Mpc.

Reference [8] combines the contributions from mass associated with galaxies in the 2M++ all-sky redshift survey [9] (linked to halos by abundance matching), resolved halos in an N-body simulation hosting galaxies too faint to be observed (the effect of which is correlated with observables to model point by point in the real Universe), and mass distributed in linear and quasilinear modes of the density field not captured by the halo model. The latter derives from the BORG algorithm [11,21–25], which propagates information from the 2M++ galaxy number density field to the current total density field and its initial conditions by means of a Bayesian likelihood framework, assuming best-fit Planck Λ cold dark matter (CDM) cosmology and marginalizing over a bias model. The parameters forming the inputs to this calculation are drawn from probability distributions describing their possible values, so that a probability distribution is derived for $\Phi_{\text{ex}}$ and $\vec{a}_5$ at the position of each test galaxy rather than a single value. Further details on our deployment of these maps may be found in D18, Secs. IV. A–IV. C, and the distributions of the input parameters themselves are summarized in Table 1 (rows 1–4) of that paper. We project $\vec{a}_5$ onto the normal of each disk (on the plane of the sky) using

$$ a_{5,z} = a_{5,\delta} \sin(\theta) - a_{5,\alpha} \cos(\theta), \quad (19) $$

where $a_{5,\alpha}$ and $a_{5,\delta}$ are the right ascension (RA) and declination (DEC) components of $\vec{a}_5$ at the test point

<figure><img src="image_9.png" /><figcaption>FIG. 4. Warp curves out to $\pm 3R_{\text{eff}}$ derived from the images in Fig. 3, with corresponding $w_1$ calculated from Eq. (18).</figcaption></figure>
and $\theta \equiv \text{SERSIC\_PHI}$ is the angle of the disk’s major axis east of north on the sky.

The total screening proxy $\Phi$ is the sum of $\Phi_{ex}$ and an internal contribution $\Phi_{in}$ due to the galaxy’s own mass. We take $\Phi_{in} = V_{disp}^2$, with $V_{disp}$ as recorded in the NSA, which we assign a 10% uncertainty as in D18 (Table 1, row 5). We draw $N_{MC} = 1000$ Monte Carlo realizations of the model, each independently sampling the $\Phi_{in}$, $\Phi_{ex}$, and $\vec{\alpha}_5$ distributions separately for each test galaxy. For a given $\lambda_C$, this gives a probability $f \equiv N(|\Phi| < |\Phi_C|)/N_{MC}$ for the galaxy to be unscreened, where $|\Phi_C| = 1.5 \times 10^{-4}(\lambda_C/32Mpc)^2$ is the screening threshold in the Newtonian potential for Hu-Sawicki $f(R)$ [26].

Next, we estimate $r_s$, $\rho_{rs}$, and $R_{eff}$ for each test galaxy. We begin by performing abundance matching (AM); see Refs [27–29]) to assign a stellar mass to each halo in the DARKSKY-400 simulation [10], postprocessed with the ROCKSTAR halo finder [30]. We use the stellar mass function of Ref. [31], which has been shown to match the stellar mass determinations of the NSA from KCORRECT, and the specific AM parametrization of Ref. [32], which is known to accurately reproduce galaxy clustering as a function of both stellar mass and r-band luminosity. We use the best-fit parameters $\alpha = 0.6$ and $\sigma_{AM} = 0.16$ for the AM proxy and universal Gaussian scatter in stellar mass at fixed halo mass, respectively. We then associate each test galaxy with the object in the abundance-matched catalog closest to it in $M_*$, taking $r_s$ from the ROCKSTAR output and calculating $\rho_{rs}$ assuming a NFW profile. As AM is inherently stochastic when $\sigma_{AM} > 0$, we generate 200 mock catalogs in this way and randomly select one for each Monte Carlo realization of our model. $R_{eff}$ is derived by multiplying the angular major-axis half-light radius of a two-dimensional (2D) Sérsic fit to the light profile (SERSIC\_TH50 in the NSA) by the angular diameter distance to the galaxy (assuming a $\Lambda$CDM cosmology $h=0.7$, $\Omega_{\Lambda}=0.7$, $\Omega_m = 0.3$, $\Omega_k = 0$) and circularizing with a further factor of $(b/a)^{1/2}$.

We now possess all the information necessary to forward model $w_1$ by Eq. (15) for each NSA galaxy as a function of $\Delta G/G_N$ and $\lambda_C$. We take a fiducial value for $n$ [Eq. (7)] of 0.5, corresponding to a central dark matter density slope that has been flattened somewhat from a primordial NFW profile by baryonic feedback. This is reasonable in light of both observations and simulations (e.g., Refs. [33,34]), although we check that varying $n$ in the range 0.1–1.5 has at most a factor of 4 effect on maximum-likelihood $\Delta G/G_N$ values.

### C. Likelihood model and inference method

Each Monte Carlo realization of our model produces a possible value of $w_1$ for given $\{\lambda_C, \Delta G/G_N\}$, so that the full set of realizations specifies the likelihood. Rather than assume a specific form for this, we approximate it empirically as a histogram by dividing the $N_{MC}$ $w_1$ values into $N_{bins} = 10$ equal-width bins between the smallest and largest values $w_{1(0)}$ and $w_{1(1)}$,

$$ \begin{aligned} & \mathcal{L}_{5,i}(w_1 | \Delta G, \lambda_C) \\ & = (1 - f)\delta(w_1) + f \sum_{j=0}^{N_{\text{bins}}-1} P_j \delta(w_1 - \Delta G w_j), \end{aligned} \quad (20) $$

where $i$ is the galaxy index, $\Delta w_1 \equiv (w_{1(1)} - w_{1(0)})/N_{\text{bins}}$ is the bin width, $w_j \equiv w_{1(0)} + (j + 1/2)\Delta w_1$ is the position of the $j$th bin center, and $P_j \equiv N_j/N_{\text{MC}}$ is the fraction of model realizations falling in bin $j$. The first term gives the screened contribution to $\mathcal{L}_5$, and the second term gives the unscreened contribution.

The full likelihood function, however, ought to contain an additional component that permits $w_1$ to be nonzero even in the absence of a fifth force; this describes measurement error in $w_1$ as well as a contribution from any other type of physics. We will see in Sec. VA that by Eq. (15) our best-fit fifth-force models produce $w_1$ values on average 2 orders of magnitude smaller than those observed, so that this additional likelihood component is in fact required to provide the bulk of the signal. For want of additional information, we take this component to be Gaussian and choose its width such that it accounts on average for the entire signal. This conservative choice guarantees that to first order the entire data set can be accounted for without fifth forces.

The simplest choice would be to take this width to be the same for all galaxies and equal to the overall standard deviation in observed $w_1$, $\sigma(w_{1,\text{obs}})$. However, this would bias our inference if the dispersion correlates significantly with galaxy properties. We find $\sigma(w_{1,\text{obs}})$ to correlate most strongly with the apparent minor-to-major axis ratio of the disk $b/a$ (Fig. 5), describing the variation in measured warp strength when the disk is inclined from edge on or has significant thickness. Neglecting this correlation would

<figure><img src="image_10.png" /><figcaption>FIG. 5. Dispersion in observed $w_1$ plotted against apparent galaxy axis ratio $b/a$ in the NSA sample. We take 30 bins between $b/a = 0.15$ and 0.5. The strong trend and large scatter indicate that this correlation must be taken into account in the non-fifth-force part of the likelihood model.</figcaption></figure>
cause $\sigma(w_{1,\text{obs}})$ to be overestimated at small $b/a$ and underestimated at large $b/a$, skewing the part of the measured $w_1$ for which the fifth-force component of the likelihood is required to account. We therefore make the width of the Gaussian noise component of the likelihood a function of $b/a$ by setting it equal to $\sigma(w_{1,\text{obs}})$ separately in 30 bins of $b/a$ between the minimum and maximum values 0.15 and 0.5. This ensures that for any $b/a$ the data can be fully accounted for to good accuracy regardless of the behavior of the fifth force, provided only that this force provides a subdominant contribution to $w_1$. To leading order, then, $w_1$ is purely noise.

Convolving $\mathcal{L}_{5,i}$ in Eq. (20) with this Gaussian (of width $\sigma_i$ for galaxy $i$) yields the total likelihood function

$$ \mathcal{L}_i(w_1|\Delta G, \lambda_C) = (1-f) \frac{\exp\{-w_1^2/2\sigma_i^2\}}{\sqrt{2\pi}\sigma_i} + f \sum_{j=0}^{N_{\text{bins}}} P_j \frac{\exp\{-(w_1 - \Delta G w_j)^2/2\sigma_i^2\}}{\sqrt{2\pi}\sigma_i}. \quad (21) $$

We find practically identical constraints on $\{\Delta G/G_N, \lambda_C\}$ if $\sigma_i$ is parametrized as a function of $b/a$ and fitted along with the fifth-force parameters. Our noise model is formally identical to that of D18 (Secs. IV.E and IV.F), except that it depends on $b/a$ rather than signal-to-noise ratio $s$.

We show three separate principal analyses. For the first two, we perform one-dimensional inference for $\Delta G/G_N$, separately for 20 logarithmically uniformly spaced values of $\lambda_C$ in the range 0.4–50 Mpc. In the first, we consider a screened fifth-force model as described above. In the second, we switch screening off so that the unscreened fraction $f$ for each test galaxy is 1 and all mass within $\lambda_C$ contributes to $\bar{a}_5$. (For this to lead to warping, the fifth force must couple to the dark matter but not the stars, which is likely contrived; we show this merely as a foil for the screening case.) For the third analysis, we again include screening and perform a full 2D inference of both $\Delta G/G_N$ and $\lambda_C$.

We check that our analysis is converged with number of Monte Carlo realizations of the model, that the AM galaxy-halo connection and smooth density field from BORG are thoroughly sampled, that our MCMC is converged with the number of steps, and that our principal results are insensitive to reasonable variations in $M(<r)$ and the assumed uncertainties in galaxy and halo properties. Further validation is documented in Sec. V C.

# V. RESULTS

## A. Comparison of observed and predicted signal

In Fig. 6, we compare the measured $|w_1|$ values and their scatter to those predicted by a model with $\lambda_C = 5$ Mpc,

$\Delta G/G_N = 1$, both with and without screening. The expectation $\langle|w_1|\rangle$ is the mode of the $N_{\text{MC}}$ realizations of the model (unscreened component only), which is the location of the bin with maximum $P_j$ [Eq. (21)]. $\sigma(|w_1|)$ is the minimal width enclosing 68% of the realizations, and $\sigma(|w_{1,\text{obs}}|)$ is the $b/a$-dependent scatter described in Sec. IV C. We see that the values predicted by this model are typically a factor of a few larger than the observed ones, indicating that if the predicted and observed warps for individual galaxies were aligned then a fifth force of this strength could account for the entire signal. By contrast, for the maximum-likelihood $\Delta G/G_N$ values of D18 (~0.02, which we will also find to be of most interest here), $\langle|w_1|\rangle$ is over an order of magnitude smaller than $|w_{1,\text{obs}}|$ on average, indicating that the better part of each warp derives from non-fifth-force effects. In this case, our approximation of Sec. IV C is valid, and the fifth-force contribution to the likelihood can be considered a small perturbation to the noise. That the $\langle|w_1|\rangle$ and $\sigma(w_1)$ distributions are similar between the screening and no screening runs at this $\lambda_C$ indicates that most of the mass sourcing $\bar{a}_5$ is unscreened in either case.

In Fig. 7, we show the correlations of the predicted and observed $|w_1|$ values with the Newtonian potential $|\Phi|$ and fifth-force acceleration $a_5$. As in D18, the predicted signal rises with $|\Phi|$ and $a_5$, cuts off in the screened case at $|\Phi| > |\Phi_c|$ when the test galaxies become screened (red), and is otherwise slightly reduced by the screening of the source masses (blue). Spearman's, $\rho_s$, and Pearson's, $\rho_p$, coefficients for the blue points are 0.60 and 0.57, respectively, for the correlation with $\Phi$ and 0.74 for the correlation with $a_5$. The green points are for the model in which screening is entirely switched off. The data, on the other hand, show no such correlations ($\rho_s, \rho_p = 0.0$; bottom panels), indicating again that a fifth force can account for at most a small part of the signal. This requires $\Delta G/G_N \ll 1$. Figure 8 plots the predicted and observed $w_1$ directly against one another.

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

We now employ the likelihood formalism of Sec. IV to constrain $\lambda_C$ and $\Delta G/G_N$. To begin, we fix $\lambda_C$ to one of its 20 values between 400 kpc and 50 Mpc and show in Figs. 9 and 10 a range of posteriors for $\Delta G/G_N$ for models both with (red) and without (green) screening. In cases in which the posteriors are peaked at $\Delta G/G_N = 0$, the model without screening is more tightly constrained because screened galaxies are predicted to have $w_1 = 0$ and hence do not contribute to the constraint.

In Fig. 11(a), we show the increase in maximum log-likelihood $\Delta \log(\mathcal{L})$ over the general relativistic (GR) case $\Delta G=0$ achieved for any $\Delta G/G_N$ at each $\lambda_C$ (solid lines; the dashed lines will be explained in Sec. V C). We see a clear peak at $\lambda_C$ similar to D18, 1.3–2.3 Mpc. Further, the best-fit $\Delta G/G_N$ values are in
<figure><img src="image_11.png" /><figcaption>FIG. 6. Histograms of the expectation (upper) and uncertainty (lower) of the fifth-force predicted (left) and observed (right) signal $w_1$ over all galaxies in our data set. $\langle |w_1| \rangle$ is the modal average of 1000 Monte Carlo realizations of the model with $\Delta G/G_N = 1$, $\lambda_C = 5$ Mpc, given ten bins of $w_1$ between its maximum and minimum values. The uncertainty is the minimal width enclosing 68% of the model realizations, which is the width of the fifth-force part of the likelihood function. For the observed signal, the uncertainty is the standard deviation of $w_1$ conditioned on apparent axis ratio $b/a$ (see Sec. IV C). This is the width of the non-fifth-force part of the likelihood function.</figcaption></figure>

good agreement; the three points with highest likelihood are $\{\lambda_C/\text{Mpc}, \Delta G/G_N, \Delta\log(\mathcal{L})\} = \{1.37, 0.019, 20\}$, $\{1.80, 0.011, 17\}$, and $\{2.32, 0.0085, 16\}$. This corresponds to a maximum AIC difference [35] of 36 for a difference of 2 degrees of freedom and hence a probability ratio of the baseline GR model to the screened model of $e^{-18}$. The $\Delta G/G_N$ posteriors are ~7σ discrepant with $\Delta G = 0$. These results are highly statistically significant even with a trial factor ~1/50 to account for the look-elsewhere effect.

The models without screening achieve only very small $\Delta\log(\mathcal{L})$ for any $\lambda_C$, indicating that the improvement in the goodness of fit depends on the partition of both test and source galaxies into screened and unscreened subsets, which is a function of both their internal mass distribution and their environment. In conjunction with the mock data tests presented in Sec. V C, we take this to be valuable

evidence concerning the validity of the fifth-force interpretation. Other effects would not be expected to correlate $w_1$ with $\Phi$ and $\tilde{a}_5$ as screening does.

The trend in best-fit $\Delta G/G_N$ with $\lambda_C$ is shown in Fig. 11(b). Models with larger $\lambda_C$ tend to give larger $\tilde{a}_5$ because they incorporate the contribution from mass out to a larger distance from a given test galaxy and hence prefer smaller $\Delta G/G_N$ for fixed $w_{1,\text{obs}}$. At larger $\lambda_C$ values than shown, $\Delta G/G_N \sim 0$.

Finally, Fig. 12 shows the constraints on $\Delta G/G_N$ and $\lambda_C$, and their covariance, for the full 2D inference in which both are free. Most of the probability density is located at the maximum-likelihood solution $\lambda_C = 1.37$ Mpc—with $0.01 \lesssim \Delta G/G_N \lesssim 0.06$—although $\lambda_C$ up to 2.3 Mpc is statistically marginally allowed. The trimodality in Fig. 12(a) likely reflects the fact that in this case most of the model realizations fell into three nonadjacent bins of
<figure><img src="image_12.png" /><figcaption>FIG. 7. Correlations of predicted (upper) and observed (lower) signals with Newtonian potential $\Phi$ (left) and fifth-force acceleration $a_5$ (right). The calculation of the prediction is as Fig. 6 for the fiducial case $\lambda_c = 5$ Mpc, $\Delta G/G_N = 1$, with $\langle\Phi\rangle$ and $\langle a_5 \rangle$ the modal averages over the 1000 model realizations. In the top row, the green points show the case in which screening is switched off, the red points are for when it is included, and the blue points are the same as the red except without the $w_1$ values of screened galaxies set to 0. The black dashed line shows the screening threshold $|\Phi_c|$ for this model. The median sizes of the error bars in the four directions (minimal 68% bounds for the prediction, Gaussian $1\sigma$ for the observations) are shown by the red crosses, although we suppress the vertical uncertainty in the bottom panels as it is simply the scatter in the points.</figcaption></figure>

$w_1$ for many galaxies, each one preferring a slightly different value for $\Delta G/G_N$ to achieve the best fit to the data.

### C. Validation

#### 1. Mock data with best-fit $\Delta G/G_N$

We turn now to methods for validating our results. We begin by generating mock data with an imposed $\Delta G/G_N$ corresponding to one of the three maximum-likelihood solutions of Fig. 11(a). We do this by first calculating $\langle w_1 \rangle$ for each galaxy—the mode of the model realizations if the galaxy is likely unscreened ($f > 0.5$) and 0 otherwise—and then scattering it by $\sigma_i$. This gives the

mock data noise characteristics identical to the real data; were the fifth-force model correct, the two data sets would also have identical correlations with the gravitational environment. We then refit the model, first recalculating $\sigma_i$ by binning the mock $w_1$ values in $b/a$, as for the real data, and then repeating the MCMC to constrain $\Delta G/G_N$ at each $\lambda_c$. In particular, we record the maximum-likelihood $\Delta G/G_N$ and the increase in log-likelihood $\Delta \log(\mathcal{L})$ over $\Delta G = 0$ that this achieves. We average these values over ten mock data sets to sample the stochasticity induced by the noise and show the results as the dashed lines in Fig. 11: the magenta line uses mock data generated with $\{\lambda_c/\text{Mpc}, \Delta G/G_N\} = \{1.37, 0.019\}$, the red line $\{1.8, 0.011\}$, the cyan line
<figure><img src="image_13.png" /><figcaption>FIG. 8. Correlation of predicted and observed $w_1$ for the model of Figs. 6 and 7.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 9. Posterior of $\Delta G/G_N$ for $\lambda_C = 1.8$ Mpc, both with and without screening. As explained in Sec. V C, we take the screened model of this Compton wavelength to provide the best overall solution to our inference.</figcaption></figure>

$\{2.32, 0.0085\}$, and the green line is as the red except refit by the model without screening.

Mock data with a signal at $\lambda_C = 1.37$ or 2.32 Mpc cannot produce $\Delta \log(\mathcal{L})$ values as high as measured from the data, indicating that in these cases the noise provides too large a contribution to the overall likelihood relative to the signal. However, this is not true for data with a signal at $\lambda_C = 1.8$ Mpc, which gives an almost identical peak in $\Delta \log(\mathcal{L})$ at $1.3 \lesssim \lambda_C$/Mpc $\lesssim 2.3$ as the real data. Further, it roughly reproduces the smaller peaks seen at $\lambda_C \approx 0.65$ and 14 Mpc. This suggests that these peaks may not be noise in our experiment [which appears from the green lines in Fig. 11(a) to be at the $\Delta \log(\mathcal{L}) \pm 1$ level] but rather arise from correlations at the test galaxy positions between the $\vec{a}_5$ and $\Phi$ fields resulting from the “real” fifth-force range $\sim 1.8$ Mpc and those which would be generated from those other ranges. As these correlations must be dependent on the locations at which the fields are sampled, we would expect them to have diminished as the number of galaxies in the sample (probing different gravitational environments) was increased. The sample size in D18 2.5 times larger than here may explain why subsidiary peaks in $\Delta \log(\mathcal{L})$ were not seen in Fig. 8(a) of that work. The use of a vector rather than scalar observable may also make it rarer for correlations between $\vec{a}_5(\lambda_C)$ and $\Phi(\lambda_C)$ at different $\lambda_C$ to propagate significantly into $\Delta \log(\mathcal{L})$.

As for the real data, fitting the mock data with a model without screening (green) achieves no noticeable $\Delta \log(\mathcal{L})$ for any $\lambda_C$, indicating either that the $\vec{a}_5$ fields of the two models are poorly correlated even at the same $\lambda_C$ or, more likely, that the inference depends crucially on the test galaxies’ unscreened likelihood fractions $f$. In the model without screening, $f \equiv 1$.

These results suggest that $\lambda_C = 1.8$ Mpc, $\Delta G/G_N = 0.011$ provides the best solution to our inference problem, and we take it as our fiducial model from now on. This is precisely the maximum-likelihood $\lambda_C$ of D18, with best-fit $\Delta G/G$ differing by a factor of only 2.3.

In Fig. 13, we generate and refit 250 mock data sets with the fiducial model to examine the scatter between them.

<figure><img src="image_15.png" /><figcaption>FIG. 10. Posteriors on $\Delta G/G_N$ for three values of $\lambda_C$, as indicated in the captions, with and without screening. In (a), the posterior for the model without screening is at very small $\Delta G$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 11. (a) shows the maximum log-likelihood increase $\Delta \log(\mathcal{L})$ over the GR case $\Delta G = 0$, for any $\Delta G/G_N$, separately for 20 $\lambda_C$ values in the range 0.4–50 Mpc (solid lines). We show results both with and without screening. The model without screening (green) never achieves $\Delta \log(\mathcal{L}) > 2$, indicating both that this model cannot improve upon GR and that the noise in our experiment is under control. The model with screening, on the other hand, (red) achieves $\Delta \log(\mathcal{L}) > 15$ for $\lambda_C = 1.3–2.3$ Mpc. The dashed lines show the average results when the models are applied to ten mock data sets with signals injected by hand at $\lambda_C = 1.37, 1.80$ or $2.32$ Mpc, as indicated in the legend, and the corresponding maximum-likelihood $\Delta G/G_N$ value 0.019, 0.011, or 0.0085. The magenta, red, and cyan dashed lines are reconstructed by the model with screening (‘s’), and the green dashed line is reconstructed by the model without screening (‘s’). That the red dashed and solid lines are similar over the entire range of $\lambda_C$ indicates that mock data with $\lambda_C = 1.8$ Mpc and $\Delta G/G_N = 0.011$ are practically indistinguishable from the real data as seen by our inference framework. The smaller peaks in $\Delta \log(\mathcal{L})$ away from the preferred $\lambda_C$, also approximately reproduced by the mock data, may indicate correlations between the $\tilde{a}_5$ fields resulting from different fifth-force ranges. (b) shows the best-fit $\Delta G/G_N$ values corresponding to the points of largest $\Delta \log(\mathcal{L})$ at each $\lambda_C$. At larger $\lambda_C$, $\Delta G/G_N \approx 0$. The dashed lines show the results from mock data with $\lambda_C = 1.8$ Mpc and $\Delta G/G_N = 0.011$, which gives similar results.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 12. Posteriors of $\Delta G/G_N$ and $\lambda_C$ from a 2D inference. (a) and (b) show the marginalized constraints, while (c) gives the full 2D posterior. The great majority of points in the chain lie in the $0.1 < \log_{10}(\lambda_C/\text{Mpc}) < 0.2$ bin.</figcaption></figure>

Figure 13(a) shows the distribution of maximum-likelihood $\Delta G/G_N$; that the reconstructed values cluster around the input value indicates that we are able to reconstruct a known truth without bias. Figure 13(c) shows the corresponding distribution of $\Delta \log(\mathcal{L})$, demonstrating the signal to be as strong relative to the noise in the mock data as the real data and hence indicating reliability of our noise model. Finally, Fig. 13(b) correlates $\Delta G/G_N$ and $\Delta \log(\mathcal{L})$, showing that mock data sets with a larger best-fit $\Delta G/G_N$ also achieve a larger $\Delta \log(\mathcal{L})$. The real data lie at the intersection of the red lines. Through the lens

of our likelihood, mock data with $\lambda_C = 1.8$ Mpc and $\Delta G/G_N = 0.011$ appear identical to the observations.

Finally, we check the $\Delta \log(\mathcal{L})$ values derived by refitting any model to mock data generated with $\Delta G = 0$. We find these to be $\lesssim 1$, showing that the peaks in Fig. 11(a) could not follow from data of this form.

## 2. Bootstrap and jackknife resampling

Next, we generate 250 mock data sets by bootstrap or jackknife resampling the real data (in the jackknife case
<figure><img src="image_18.png" /><figcaption>FIG. 13. Results of refitting a screened model with $\lambda_C = 1.8$ Mpc to 250 mock data generated by that model at $\Delta G/G_N = 0.011$, differing only in noise. (a) plots the distribution of best-fit $\Delta G/G_N$ values, and (c) plots the corresponding distribution of $\Delta \log(\mathcal{L})$. The vertical red lines show the values in the real data (0.011 and 17.4, respectively); that these are near the centers of the mock data distributions shows those data to behave very much like the observations. (b) illustrates the correlation between $\Delta G/G_N$ and $\Delta \log(\mathcal{L})$ among the mock data sets, with the real data located at the intersection of the red lines.</figcaption></figure>

retaining a random 70% of the sample). We fit each data set with the screened $\lambda_C = 1.8$ Mpc model and calculate both the maximum-likelihood $\Delta G/G_N$ and the discrepancy from $\Delta G = 0$ exhibited by the posterior, defined as the ratio of the median $\Delta G/G_N$ to the standard deviation. We show the resulting histograms in Fig. 14, in which the red lines show

the corresponding values in the real data. Most bootstraps have a very similar significance of $\Delta G > 0$ to the real data, and most jackknives have a slightly lower significance due to the reduced statistics. There is, however, a sizeable fraction of data sets in each case with much lower significance, $<3\sigma$. Thus, in some samples similar to ours,

<figure><img src="image_19.png" /><figcaption>FIG. 14. (a) and (d) show the distributions of maximum-likelihood $\Delta G/G_N$ values (for the screened model with $\lambda_C = 1.8$ Mpc) inferred from 250 bootstrap (upper) or jackknife (lower) resamples of the NSA data. Each jackknife data set contains 70% of the full sample. The vertical red line shows the result in the full sample (0.011), around which the resample values cluster. (c) and (f) show the corresponding distributions of deviations from $\Delta G = 0$, measured by dividing the median of the $\Delta G/G_N$ posterior by its standard deviation. Although both sets produce $\Delta G > 0$ with lower significance than the full data set on the whole, that the majority of resamples achieve a $\gtrsim 3\sigma$ result indicates that the preference for a fifth force is not a peculiar property of the NSA sample. Finally, (b) and (e) show the correlation between best-fit $\Delta G/G_N$ and discrepancy from $\Delta G = 0$. The real data lie at the intersection of the red lines.</figcaption></figure>
a significant deviation from $\Delta G = 0$ would not be found. The maximum-likelihood $\Delta G/G_N$ values tell a similar story: while resamples in which a significant improvement over GR is possible favor $\Delta G/G_N$ values close to that in the total data and the remainder favor $\Delta G \simeq 0$. Of all the checks we have performed, this one casts the most doubt on the fifth-force interpretation and indicates that more data may be required to validate the model with high statistical confidence. That our results here are more ambiguous than the corresponding Fig. 12 of D18 likely follows from our smaller sample size.

### 3. Inversion of the predicted signal

As a final check, we repeat the inference at $\lambda_C = 1.8$ Mpc with the sign of all predicted $w_1$ values reversed. The posterior of $\Delta G/G_N$ is shown in Fig. 15 (cf. Fig. 9): that the detection goes away shows that the directions of the measured and predicted warps are on the whole aligned on a galaxy-by-galaxy basis when $w_1$ has the correct sign, as required for a component of each warp to derive from $\vec{a}_5$.

Although the present analysis uses fewer galaxies than D18 (4,206 vs 10,822), the smaller uncertainty in both the measurements and theoretical prediction (depending on the halo density out to $\sim r_s$, which is better known than in the central $\sim 100$ pc) allows it to achieve a slightly stronger result. This is quantified by both the maximum $\Delta \log(\mathcal{L})$ (20.2 vs 15.6) and the significance of the $\Delta G/G_N$ posterior from 0 at the maximum-likelihood $\lambda_C$ (7.0$\sigma$ vs 6.6$\sigma$). Further, if the scatter in $\Delta \log(\mathcal{L})$ across $\lambda_C$ for the unscreened model is indicative of its uncertainty, $\Delta \log(\mathcal{L})$ is known here to $\pm 1$ [Fig. 11(a)], while in D18, it is known only to $\pm 4$ (their Fig. 8(a)). Finally, the width of the primary peak in $\Delta \log(\mathcal{L})$ is narrower here than in D18, providing a tighter constraint on $\lambda_C$. The present inference seems somewhat cleaner than D18's.

<figure><img src="image_20.png" /><figcaption>FIG. 15. Posterior of $\Delta G/G_N$ for the screened model with $\lambda_C = 1.8$ Mpc, but with the sign of the predicted $w_1$ reversed (cf. Fig. 9).</figcaption></figure>

## VI. DISCUSSION

### A. Systematics

Our principal systematics are as follows. Many are directly analogous to those of D18, Sec. VI B, in which further details may be found.

We assume that the component of the warps due to anything other than a fifth force can be described by a Gaussian with the width dependent only on apparent axis ratio $b/a$. This describes the appearance of warps on the sky rather than their intrinsic properties. Although the astrophysical processes governing warp formation are likely to correlate warp strength with a number of internal and environmental galaxy properties, the exact nature of these correlations is unknown due to small sample sizes in previous warp studies and uncertainty in the underlying theory. Unless these correlations are similar to those induced by fifth forces, however, they would not be expected to significantly bias our inference. This seems unlikely due to the specificity of the dependence of predicted $w_1$ on halo properties $n, \rho_{rs}$, and $r_s$ and galaxy size $R_{\text{eff}}$ in Eq. (15) as well as the surrounding matter that sets $\Phi_{\text{ex}}$ and $\vec{a}_5$. Further, the U-shaped warps expected from fifth forces and to which we tailor our summary statistic are observationally less common than S shapes and are therefore less likely to receive large contributions from competing astrophysical effects.

We see two ways in which the impact of non-fifth-force physics may be assessed. The first is to measure or derive a set of internal and external properties for a statistical sample of galaxies that may be expected to correlate with warp strength. By calculating the covariances between $w_1$ (or any other summary statistic) and each of these properties, the dependencies of galaxy warping may be empirically determined. This full covariance matrix may then replace the simple $b/a$-dependent variance in the non-fifth-force part of the likelihood function, allowing one to check that a significant correlation with $\vec{a}_5$ remains even when the noise model incorporates finer-grained correlations. We discuss this further in Sec. VI C. The second is to measure warp strengths in a population of galaxies produced by hydrodynamical simulations, which in principle incorporate all relevant baryonic physics in $\Lambda$CDM. It is this population that the non-fifth-force part of the likelihood function ought to describe. Unfortunately, this requires $\mathcal{O}(10 \text{ pc})$ spatial resolution in the central regions of halos, which is currently possible only in zoom simulations with poor statistics.

Our determination of the $\Phi$ and $\vec{a}_5$ fields and the galaxy-halo connection assume a concordance $\Lambda$CDM cosmology, which is inconsistent in detail with that implied by our fifth-force models (e.g., Refs. [36–39]). Nevertheless, the weakness of the fifth force in our preferred scenarios [$\mathcal{O}(1\%)$ that of gravity, with a range $\mathcal{O}(1 \text{ Mpc})$] makes changes to the growth rate of structure and abundance and profiles of halos small. Ultimately, it will be preferable to repeat
our analysis self-consistently with modified gravity
simulations [40,41].

Following Ref. [15], we neglected the self-gravity of the
disk in calculating the warp curve given an offset between
the halo and stellar mass centers. The action of neighboring
stellar rings provides a restoring force on the disk, reducing
the magnitude of the warp in equilibrium. Incorporating
this properly requires the sophisticated machinery of
potential theory [42], but we can provide an order-of-
magnitude estimate of its effect by simply summing the
forces due to the disk elements in the warp configuration of
Sec. IIC and comparing the resultant to $\vec{a}_5$. We show this
calculation in the Appendix; the upshot is that the disk self-
force is $O(10^{-2})$ times that from the halo, leading to an
$O(1\%)$ reduction in predicted $w_1$. This is almost certainly
too small to impact our inference.

Finally, we have assumed that the warp curve (10) would
be observed at all times, i.e., that warps due to fifth forces
are quasistable phenomena. This is likely accurate provided
that the fifth-force field driving the warp varies on time-
scales longer than the dynamical time required for the
equilibrium of the disk to become established. However,
the excitation modes of the disk required for the main-
tenance of a warp are not fully understood, and it may be
possible for the observed warp curve to vary significantly
with time even for fixed $\vec{a}_5$. This is likely to be more
important when the warp is forced by one or a few
individual massive objects, which will induce stronger
tidal effects and perturbations to the galaxy's velocity than
the case in which $\vec{a}_5$ receives small contribution from many
masses within $\lambda_C$. Fortunately, unscreened objects—which
are the only ones to experience a fifth-force warp at all—
fall mostly in the latter class.

B. Comparison with D18

Figure 11 is very similar to D18’s Fig. 8. Both analyses
find 6–7σ evidence for a chameleon or symmetron-
screened fifth force with range 1.3 < λ<sub>C</sub>/Mpc < 2.3 and
strength 0.01 ≲ ΔG/G<sub>N</sub> ≲ 0.03. There is a small offset
between the maximum-likelihood λ<sub>C</sub> values of the two
analyses—1.4 Mpc here vs 1.8 Mpc in D18—although it is
within the uncertainties. Our analysis of mock data with an
injected ΔG/G<sub>N</sub> (Sec. V C) suggests that 1.8 Mpc is in fact
the preferred solution here also. The best-fit ΔG/G<sub>N</sub> values
are somewhat farther apart—0.011 here vs 0.025 in D18 at
λ<sub>C</sub> = 1.8 Mpc—and formally discrepant at just over 3σ.

These differences, however, seem trifling compared to
the general agreement that we find. This is all the more
remarkable given that the two analyses use largely inde-
pendent samples of galaxies and signals with qualitatively
different dependences on galaxy and halo structural param-
eters. Although the environmental screening (Φ<sub>ex</sub>) and
fifth-force acceleration (ā<sub>5</sub>) fields are common to both
analyses, they are evaluated at different points for different
samples, making their values largely independent. As they

are sourced by masses within $\lambda_C \approx 1.8$ Mpc, regions of
the field farther apart than this are uncorrelated. The total
degree of screening depends also on $\Phi_{\text{int}}$, which is
determined differently in the two analyses. Finally, the
relevant projections of $\vec{a}_5$ also differ: in D18, $\vec{a}_5$ was
projected onto the plane of the sky to yield independent RA
and DEC projections of the vector signal $\vec{r}_*$, while here it is
projected onto the disk normals, to determine warping, and
then onto the plane of sky to set the observed signal $w_1$. We
therefore do not believe there to be significant covariance
between these two analyses; their results ought to be largely
independent.

We caution, however, that systematic uncertainties may
impact the best-fit values of ΔG/G<sub>N</sub> and λ<sub>C</sub>. The largest in
either analysis is likely the central halo density ρ<sub>0</sub> that
appears in D18’s Eq. (17) and especially its scatter between
galaxies. Figure 16 of that paper shows that varying this
scatter between 0.4 and 1.4 dex causes a variation in best-fit
ΔG/G<sub>N</sub> of 3 orders of magnitude, completely dwarfing the
factor of 2.3 difference with the result here. The present
analysis is much less prone to uncertainties of this type
because w<sub>1</sub> is determined by the entire extent of the disk
(few kpc), rather than the very central regions that define r<sub>*</sub>
(few × 10<sup>-2</sup> kpc). This is manifest in the appearance of ρ<sub>rs</sub>,
the dark matter density at the halo scale radius, in the
formula for expected w<sub>1</sub> where ρ<sub>0</sub> appears in the formula
for r<sub>*</sub>. Thus, we are able to determine the scatter in the
restoring force due to the dark matter self-consistently by
modeling the full galaxy-halo connection with AM, rather
than having to impose it by hand. The approximate
agreement in the reconstructed ΔG/G<sub>N</sub> between the two
analyses may therefore indicate that the scatter in ρ<sub>0</sub> among
galaxies of the same central surface baryonic mass density
is indeed ~1 dex, as we assumed in D18. Intriguingly,
under a unified model for both warps and gas-star offsets—
which may or may not invoke a screened fifth force—the
consistency of the signals could provide a new handle on
halo properties, potentially enabling them to be determined
with much greater precision than is currently possible.

While the inferred ΔG/G<sub>N</sub> is sensitive to the average
magnitude of the predicted signal at fixed ΔG/G<sub>N</sub>, the
inferred λ<sub>C</sub> is sensitive only to the variation of Φ and ã<sub>5</sub>
with the radius from a test point out to which the
contributions from masses are summed. In particular, λ<sub>C</sub>
is almost fully insensitive to the detailed distribution of
dark matter mass within halos. We therefore believe our
constraints on λ<sub>C</sub> to be significantly more robust than
those on ΔG/G<sub>N</sub> and consider similarity between the
λ<sub>C</sub> − Δlog(L) relations here and in D18 a more important
indicator of the general consistency of the analyses than
similarity between the maximum-likelihood ΔG/G<sub>N</sub> val-
ues. Although we cannot conclude that a screened fifth
force must be responsible for the signals, we can be
confident that whatever is responsible operates on a scale
of ~1.8 Mpc and determines the signals through a feature
common to both, most likely a separation of galaxy mass components.

### C. Previous work and the path ahead

The astronomical study of galaxy warps has been ongoing for several decades (see Refs. [14,43] for reviews). The warp in 21 cm emission of the Milky Way has been known since the 1950s [44,45], and warps of other nearby massive galaxies have been known since the 1970s [46,47]. By the 1990s, sample sizes were large enough to conclude that the majority of galaxies' HI disks were probably warped [48]. Evidence for warping of stellar disks (such as we investigate here), however, has proven harder to accumulate as optical emission peters out at a smaller galactocentric radius than HI. Statistical samples of warps in optical disks have been compiled only more recently [49]. The degree of warping has typically been investigated qualitatively by means of visual inspection, although automated algorithms along the lines of ours have previously been employed [50].

Despite empirical evidence for the prevalence of warps, the processes which drive warp formation are not well understood. In particular, there remains considerable uncertainty on whether warps are generated primarily by factors internal to galaxies or rather by environmental influences [49,51]. Those in the latter class include tidal interactions between galaxies, intergalactic magnetic fields, and accretion of gas or dark matter, while those in the former rely on interaction between the disk and halo, perhaps due to a triaxial halo misaligned with the disk. Although most warped disks have neighbors massive enough to induce strong tidal fields, warps have also been observed in apparently isolated systems [52]. The lifetime of warps is also unclear; while galaxies unlikely have discrete warping modes, they may have long-lived collective excitations resembling warps that are maintained by infall or other interaction with the environment [53]. We can provide only a small piece of information on these puzzles here: that our models without screening achieve no increase in $\log(\mathcal{L})$ over the case $\Delta G = 0$ for any $\{\lambda_C, \Delta G/G_N\}$ indicates that warps are not significantly correlated with the local acceleration field sourced by matter within any distance from 0.4 to 50 Mpc.

The general mechanisms underlying warp formation have ramifications for our inference: the “true” (effective) model of astrophysical warps should determine the non-fifth-force part of our likelihood function. Exactly analogously to D18, the Gaussian component that we convolve with $\mathcal{L}_5$ describes both these additional physical effects and observational uncertainty in the determination of $w_1$. Our assumption of a Gaussian (with width dependent only on $b/a$) will be more accurate the greater the relative importance of the observational uncertainty, which it is better suited to describe. We cannot estimate these uncertainties

here, requiring us to extract the information directly from the data.

While most astrophysical processes are expected to generate S- (or integral-) rather than U-shaped warps (and indeed these are more common observationally) [14,54], they almost certainly generate a U-shaped component with strength that correlates significantly with other galaxy properties. In principle, it is possible to use our determinations of warp and galaxy properties to investigate this issue and account in part for additional correlations. For example, simply correlating $w_1$ with a number of internal and environmental galaxy properties such as we have already estimated would provide a handle on which are most important. The internal properties are the observed galaxy size, mass and rotation velocity, and the halo mass and concentration from AM; the external properties are the Newtonian potential, acceleration, curvature and ambient density from our gravitational maps. Lack of significant correlation with these properties would indicate either that variations in the warp strength are driven mainly be observational error, or that the astrophysical processes driving warp formation are largely stochastic. In either case, our likelihood model for the noise would be reliable. Implementing the remaining correlations in the noise model would give a more accurate baseline on which to add the fifth-force signal.

In the context of modified gravity, warps have been investigated previously in Ref. [55], using 495 galaxies from SDSS and ALFALFA. Its authors were unable to place significant constraints due to the relatively small sample size and inability to forward model the warp statistic, but roughly reach sensitivity to $\lambda_C$/Mpc ~ few, $\Delta G/G_N$ ~ few. They forecast that ~8000 isolated dwarf galaxies would be required to test $f_{R0}$ to the $2 \times 10^{-7}$ level, i.e., $\Delta G/G_N = 1/3$, $\lambda_C = 1.4$ Mpc. Due to advances in the screening maps, likelihood function, and inference machinery, we achieve over an order of magnitude greater sensitivity to $\Delta G/G_N$ at that $\lambda_C$ with only ~4000 galaxies, not all of which are isolated dwarves. Our inferences greatly surpass in strength all previous fifth-force searches on either astrophysical [55–58] or cosmological [37,59–69] scales, and the signal they favor has not been constrained by any previous test.

Besides investigating our systematics further—including possible contamination from galaxy formation physics—we plan to extend our framework to the remaining intra-galaxy signals identified in Refs. [15,19]: offsets and asymmetries in the rotation curves of stars and gas [55,58]. It is possible to make predictions for these from the best-fit models here and in D18 by means of either the maximum-likelihood estimator or Bayesian posterior predictive distribution; verification of this would provide stronger evidence than the reconstructions presented thus far. Ideally, one would infer the model parameters from all of these data sets jointly. We will also improve the
resolution of our gravitational maps with constrained N-body simulations using the initial conditions inferred by BORG. Finally, repeating the analyses with independent data of greater quality and quantity in the future—e.g., from DES, WFIRST, Euclid, or LSST—would provide a further check on the results as well as increase the precision of the inference.

In addition to galaxy formation and dynamics, and fifth forces, warps contain information on possible dark matter interactions which have been proposed to resolve small-scale problems in ΛCDM [70,71]. For example, a dark matter self-interaction (SIDM) generates a drag force on a halo as it falls through another halo or the ambient medium, causing a separation between halo and galaxy centers directly analogously to a screened fifth force. This has so far been studied analytically [72,73] and in simulations [20] and implies that dark matter self-interactions constitute another systematic for our fifth-force search: the two types of new physics are degenerate in their effect on $w_1$. This is mitigated, however, by the fact that SIDM would be expected to induce warps pointing predominantly in directions opposite to those from fifth forces, viz. parallel to a galaxy’s peculiar velocity rather than opposite the external fifth-force field. That $w_1$ is on the whole anti-correlated with $\vec{a}_5$ as expected from a screened fifth force therefore implies net antialignment with the SIDM prediction and hence the potential for strong constraints on the interaction cross section $\sigma/m_{\text{dm}}$. We place these constraints quantitatively in upcoming work.

Regardless of the fate of the detection here and in D18, we foresee galaxy-scale tests of fundamental physics becoming a major frontier of astrophysics in the next decade.

## VII. CONCLUSION

With the advent of massive galaxy surveys, intragalaxy probes are set to become increasingly important in our search for new physics. These are complementary to traditional cosmological probes of ΛCDM and modified gravity, come at little if any extra experimental cost, and reach sensitivities to well-motivated screening mechanisms far surpassing those of any other test. In this paper, we search for a chameleon- or symmetron-screened fifth force—as well as an unscreened fifth force with differential coupling to stars and dark matter—by investigating the warping of stellar disks that it predicts. This follows from the spatial offset that develops between the stellar disk and halo center when only one of the two feels the fifth force; in the case of screening, this arises when stars self-screen in otherwise unscreened galaxies.

To probe this effect, we reduce 4,206 r-band images of late-type galaxies in the local 100 Mpc from the *Nasa Sloan Atlas* and develop an automated algorithm for quantifying the degree of U-shaped warping on the plane of the sky. We then calculate the expectation for the warp curve of each

galaxy as a function of the fifth-force range $\lambda_C$ and strength relative to gravity $\Delta G/G_N$, the screening field $\Phi$ and external fifth-force field $\vec{a}_5$, and the internal structural properties of the galaxies and their halos. By propagating the full non-Gaussian uncertainties in each of these inputs by Monte Carlo sampling, we build the fifth-force likelihood function for the observed data $d$, $\mathcal{L}_5(d|\lambda_C, \Delta G/G_N)$, which we convolve with a Gaussian noise component describing both measurement uncertainty and the contribution from non-fifth-force physics. As the expected warp for $\Delta G/G_N$ values of interest, $O(0.01)$, is only $O(1\%–10\%)$ of that measured, the non-fifth-force part of the model accounts for the great majority of the signal. This enables us to set its width equal to the standard deviation of the measured warp in bins of observed axis ratio $b/a$, ensuring that the measurements can be approximately fully accounted for without novel physics. This should make our analysis conservative.

Feeding the data into our likelihood framework and applying Bayes’s theorem, we constrain $\lambda_C$ and $\Delta G/G_N$ by Markov chain Monte Carlo. The overall likelihood is significantly increased ($\Delta \log(\mathcal{L}) \approx 20$) over the General Relativistic case $\Delta G = 0$ by adding a screened fifth force with $\lambda_C \approx 1.8$ Mpc, $\Delta G/G_N \approx 0.01$. At the maximum-likelihood $\lambda_C$, the $\Delta G/G_N$ posterior is inconsistent with 0 at $7\sigma$. Neighboring $\lambda_C$ values exhibit a coherent trend in $\Delta \log(\mathcal{L})$, precisely as would be expected in a fifth-force scenario due to the correlation of $\vec{a}_5$ and $\Phi$ with the fifth-force range. Conversely, models of any fifth-force range but without screening achieve no significant increase in likelihood over GR, again as expected from the correlations inherent to the screened model. We check that mock data without a fifth force could not generate our $\{\lambda_C, \Delta G/G_N\}$ posteriors, that the $\Delta \log(\mathcal{L})$ from the real data is compatible at all $\lambda_C$ with that from mock data generated using the best-fit model, that bootstrap and jackknife resamples yield results similar to the full data set, and that the preference for $\Delta G > 0$ is eliminated by rotating the predicted warp through $180^\circ$.

These results are very similar to those obtained with largely independent data and a completely different signal in Ref. [13]. We cannot be sure that these results could not arise from “galaxy formation” physics in ΛCDM, but we can be confident that some physical process, with a preferred scale and specific correlation with the gravitational environment, is separating galaxies’ stellar, gas, and dark matter mass components. At the very least, our analyses demonstrate the immense and largely undiscovered power of intragalaxy signals for probing fundamental physics.

## ACKNOWLEDGMENTS

This work was supported by St John’s College, Oxford. P. G. F. acknowledges support from Leverhulme, Science and Technology Facilities Council, Beecroft Institute for
Particle Astrophysics and Cosmology, and the European Research Council. G. L. acknowledges support by L'Agence Nationale de la Recherche Grant No. ANR-16-CE23-0002 and from the Labex Institut Lagrange de Paris (reference ANR-10-LABX-63) part of the Idex Sorbonne Universités à Paris pour l'Enseignement et la Recherche (Grant No. ANR-11-IDEX-0004-02). We thank James Binney and Kris Pardo for helpful discussions on warps. Computations were performed at Oxford and SLAC.

## APPENDIX: THE EFFECT OF DISK SELF-GRAVITY

In deriving $w_1$ expected from a fifth force in Sec. II C [Eq. (15)], we assumed that the only restoring force on the stellar disk came from the halo. In reality, however, each stellar ring interacts with its neighbors to provide an additional force that damps the warp. In this Appendix, we estimate the fractional contribution to the overall restoring force from the disk in the warp configuration of Eq. (10). If this is small—as we shall find it to be—then the perturbation to $w_1$ due to the disk is not important.

Take a point on the disk at $(x_0, y_0, z_0)$, and calculate the force in the $z$ direction due to the rest of the disk. The acceleration from an area element $dA$ at $(x, y, z)$ is

$$da_z = \frac{G_N \Sigma(x, y, z) dA(z - z_0)}{((x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2)^{3/2}}, \quad (A1)$$

where $\Sigma$ is the surface stellar mass density. As there is azimuthal symmetry around the $z$ axis, we can, without loss of generality, take $y_0 = 0$ and specify a general point on the disk with cylindrical polar coordinates $(r, \phi)$:

$$x = r \cos(\phi), \quad y = r \sin(\phi), \quad z = Kr^n. \quad (A2)$$

We also have $dA = rd\phi dl$, where $dl$ is an infinitesimal displacement along the disk:

$$dl = (dr^2 + dz^2)^{1/2} = dr(1 + K^2 n^2 r^{2(n-1)})^{1/2}. \quad (A3)$$

Using that $\Sigma$ is independent of $\phi$ by symmetry, the total disk self-force in the $z$ direction is then

$$a_z = G_N K \int_0^\infty dr \int_0^{2\pi} d\phi \times \frac{r\Sigma(r)(r^n - r_0^n)(1 + K^2 n^2 r^{2(n-1)})^{1/2}}{((r\cos(\phi) - x_0)^2 + r^2\sin(\phi)^2 + K^{2n}(r^n - r_0^n)^2)^{3/2}}. \quad (A4)$$

We now make two simplifications that maintain full generality. First, we take the ratio, $f_a$, of $a_z$ with the fifth-force acceleration $a_{5,z}$ that generates the warp. Second, we express $K$ in terms of $w_1$ using Eq. (14):

$$KG_N/a_{5,z} = \frac{3-n}{4\pi\rho_{rs}r_s^n} = w_1 \frac{(n+1)(n+2)}{n} a_{5,z}^{-1} (3R_{\text{eff}})^{1-n}. \quad (A5)$$

Thus,

$$f_a = w_1 G_N \frac{(n+1)(n+2)}{n} a_{5,z}^{-1} (3R_{\text{eff}})^{1-n} I(n, r_0), \quad (A6)$$

where $I$ is the integral over $r$ and $\phi$ in Eq. (A4).

From Eq. (A6), the order of magnitude of $f_a$ is already apparent on dimensional grounds. As $I$ has units of $ML^{n-3}$, and the only dimensions in the integral are $M_*$ and $R_{\text{eff}}$ (which sets $\Sigma$), we must have $I = AM_* R_{\text{eff}}^{n-3}$, where $A$ is a constant of proportionality. Thus, dimensionally,

$$f_a \sim w_1 \frac{M_* G_N}{R_{\text{eff}}^2 a_{5,z}}. \quad (A7)$$

Plugging in some characteristic values for our galaxies—$w_1 = 10^{-4}$, $M_* = 10^9 M_\odot$, $R_{\text{eff}} = 3$ kpc, $a_{5,z} = 10^{-16}$ km s$^{-2}$, and $G = 4.5 \times 10^{-39}$ kpc$^3$ s$^{-2}$ M$^{-1}$—we find $f \sim 10^{-2}$. This is the approximate factor by which the warp curve $z(r)$ is flattened by the self-gravity of the disk; the change is likely insignificant.

A more precise solution is obtained by evaluating Eq. (A4). Given that the warp is small, we can expand the integral in powers of $z/r$; the leading-order contribution is when all $K$-dependent terms are dropped. Assuming an exponential disk

$$\Sigma(r) = \frac{M_*}{2\pi R_d^2} \exp(-r/R_d), \quad (A8)$$

we find

$$I(n, r_0) = \frac{M_*}{2\pi R_d^2} \int_0^\infty dr \int_0^{2\pi} d\phi r \exp(-r/R_d)(r^n - r_0^n) \times ((r \cos(\phi) - r_0)^2 + r^2 \sin^2(\phi))^{-3/2}. \quad (A9)$$

The largest value will arise for $r_0 = 0$ (the center of the disk), since in that case all of the disk is at larger $z$ than the point under consideration. This will therefore give the largest value of $f_a$ and the most stringent test of the approximation that the disk self-gravity is negligible. In this case, the $\phi$ integral trivially evaluates to $2\pi$, and we are left with

$$I(n, 0) = \frac{M_*}{R_d^2} \int_0^\infty dr \exp(-r/R_d) r^{n-2}. \quad (A10)$$

e.g., for $n = 1.5$, using that $R_{\text{eff}} = 1.68 R_d$, we find $I = 3.87 R_{\text{eff}}^{-1.5}$. Thus,

$$f_a = 13 w_1 \frac{M_* G_N}{R_{\text{eff}}^2 a_{5,z}}. \quad (A11)$$
[1] C. Brans and R. Dicke, Phys. Rev. **124**, 925 (1961).

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

[3] T. Markkanen and A. Tranberg, J. Cosmol. Astropart. Phys. 08 (2013) 045.

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

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

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

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

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

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

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

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

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

[13] H. Desmond, P. G. Ferreira, G. Lavaux, and J. Jasche, Phys. Rev. D **98**, 064015 (2018).

[14] J. Binney, Annu. Rev. Astron. Astrophys. **30**, 51 (1992).

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

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

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

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

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

[20] L. F. Secco, A. Farah, B. Jain, S. Adhikari, A. Banerjee, and N. Dalal, Astrophys. J. **860**, 32 (2018).

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

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

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

[24] J. Jasche, F. Leclercq, and B. D. Wandelt, J. Cosmol. Astropart. Phys. 01 (2015) 036.

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

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

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

[28] 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).

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

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

[31] M. Bernardi, A. Meert, R. K. Sheth, V. Vikram, M. Huertas-Company, F. Mei, and F. Shankar, Mon. Not. R. Astron. Soc. **436**, 697 (2013).

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

[33] W. J. G. de Blok, Adv. Astron. **2010**, 1 (2010).

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

[35] H. Akaike, IEEE Trans. Autom. Control **19**, 716 (1974).

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

[37] L. Lombriser, Ann. Phys. (Berlin) **526**, 259 (2014).

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

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

[40] 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).

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

[42] C. Hunter and A. Toomre, Astrophys. J. **155**, 747 (1969).

[43] K. Kuijken and I. Garcia-Ruiz, in *Disks and Disk Galaxies, Proceeding of a conference held in Rome, 2000, ASP Conference Series* Vol. 230, (Astronomical Society of the Pacific, San Francisco, 2001), pp. 401–408.

[44] B. F. Burke, Astron. J. **62**, 90 (1957).

[45] F. J. Kerr, Astron. J. **62**, 93 (1957).

[46] K. Newton and D. T. Emerson, Mon. Not. R. Astron. Soc. **181**, 573 (1977).

[47] D. H. Rogstad, M. C. H. Wright, and I. A. Lockhart, Astrophys. J. **193**, 309 (1974).

[48] A. Bosma, Ph.D. thesis, Groningen University, 1978.

[49] V. Reshetnikov, E. Battaner, F. Combes, and J. Jiménez-Vicente, Astron. Astrophys. **382**, 513 (2002).

[50] J. Jiménez-Vicente, C. Porcel, M. L. Sánchez-Saavedra, and E. Battaner, Astrophys. Space Sci. **253**, 225 (1997).

[51] N. Castro-Rodríguez, M. López-Corredoira, M. L. Sánchez-Saavedra, and E. Battaner, Astron. Astrophys. **391**, 519 (2002).

[52] R. Sancisi, Astron. Astrophys. **53**, 159 (1976).

[53] L. S. Sparke and S. Casertano, Mon. Not. R. Astron. Soc. **234**, 873 (1988).

[54] I. García-Ruiz, R. Sancisi, and K. Kuijken, Astron. Astrophys. **394**, 769 (2002).

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

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

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

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

[59] 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).

[60] J. Dossett, B. Hu, and D. Parkinson, J. Cosmol. Astropart. Phys. 03 (2014) 046.

[61] S. Ferraro, F. Schmidt, and W. Hu, Phys. Rev. D **83**, 063503 (2011).

[62] L. Lombriser, F. Schmidt, T. Baldauf, R. Mandelbaum, U. Seljak, and R. E. Smith, Phys. Rev. D **85**, 102001 (2012).
[63] L. Lombriser, A. Slosar, U. Seljak, and W. Hu, Phys. Rev. D **85**, 124038 (2012).

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

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

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

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

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

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

[70] N. E. Mavromatos, C. R. Argüelles, R. Ruffini, and J. A. Rueda, Int. J. Mod. Phys. D **26**, 1730007 (2017).

[71] D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. **84**, 3760 (2000).

[72] F. Kahlhoefer, K. Schmidt-Hoberg, M. T. Frandsen, and S. Sarkar, Mon. Not. R. Astron. Soc. **437**, 2865 (2014).

[73] J. Kummer, F. Kahlhoefer, and K. Schmidt-Hoberg, Mon. Not. R. Astron. Soc. **474**, 388 (2018).