Reconstructing the gravitational field of the local universe

TL;DR

This paper maps key gravitational variables in the local universe using galaxy surveys and simulations to evaluate their potential in testing gravity theories and constraining models.

Contribution

It introduces a method combining observational data and simulations to map gravitational variables around galaxies up to 200 Mpc, enhancing gravity tests.

Findings

Galaxies significantly extend the scope of gravity tests.

Mapping these variables helps constrain specific gravity models.

Upcoming surveys will improve these constraints.

Abstract

Tests of gravity at the galaxy scale are in their infancy. As a first step to systematically uncovering the gravitational significance of galaxies, we map three fundamental gravitational variables -- the Newtonian potential, acceleration and curvature -- over the galaxy environments of the local universe to a distance of approximately 200 Mpc. Our method combines the contributions from galaxies in an all-sky redshift survey, halos from an N-body simulation hosting low-luminosity objects, and linear and quasi-linear modes of the density field. We use the ranges of these variables to determine the extent to which galaxies expand the scope of generic tests of gravity and are capable of constraining specific classes of model for which they have special significance. Finally, we investigate the improvements afforded by upcoming galaxy surveys.

Tables

Table 1. Table 1: Statistics describing the distributions of correction factors relating proxy values from 2M++ objects to those from all halo mass. “Calibration” refers to our procedure for estimating the correction factors from { d 𝑑 d , N 10 subscript 𝑁 10 N_{10} , Φ vis subscript Φ vis \Phi_{\text{vis}} } (Section  2.2 ), which can be seen to reduce their uncertainties. Note that the first three columns include all well-resolved halos in the simulation, while the last pertains to the correction factors actually assigned at the positions of the 2M++ galaxies (as in Fig.  1 ).

Correction factor	Median	Mean	St. dev.	St. dev. after calibration
$c_{\Phi }$	1.1	1.5	1.4	0.18
$c_a$	1.1	3.9	58	0.53
$c_{\theta }$	0.07	0.29	0.51	0.33
${\log }_{10}\left(c_K\right)$	0.09	0.45	0.83	0.14

Table 2. Table 2: Parameters forming the input to our calculations of Φ Φ \Phi , a → → 𝑎 \vec{a} and K 𝐾 K , and the source of the distributions from which they are drawn. In this work we take approximate maximum-likelihood values for each; in future work full uncertainties will be estimated by Monte Carlo marginalisation. Subscripts ‘t’ and ‘o’ in the final row denote true and observed values respectively.

Input	Origin of parent distribution
$\overrightarrow{c}\equiv \{c_{\Phi },c_a,c_{\theta },c_K\}$	Regions of N-body simulations with
	given distributions of visible mass
Smoothed density field	Galaxy number density and peculiar
	velocity fields (Lavaux & Jasche, 2016)
$P(M_{\text{vir}},c|L;\alpha ,{\sigma }_{\text{AM}})$	Abundance matching mock with
	fixed halo proxy and scatter
$P(\alpha ,{\sigma }_{\text{AM}})$	Galaxy clustering analysed by AM
	(Lehmann et al., 2015)
$P(L_{\text{t}},{\overrightarrow{x}}_{\text{t}}|L_{\text{o}},{\overrightarrow{x}}_{\text{o}})$	Galaxy survey measurements

Equations

$${\lambda }_C\approx 32\sqrt{f_{R0}/10^{-4}}\text{Mpc},$$

$${\Phi }_{\text{crit}}/c^2\approx \frac{3}{2}{f}_{R0}.$$

$$M=m-A_k(l,b)-k(z)+e(z)-D_L(z),$$

$$\rho (r)=\frac{M_{\text{vir}}{c}^2}{4\pi r(R_{\text{vir}}+cr{)}^2(\ln (1+c)-c/(1+c))},$$

$${\Phi }_{\text{vis}}=-\sum _{\text{i}}\frac{G{M}_{\text{vir,i}}}{r_{\text{i}}}\frac{\ln (1+\frac{c_{\text{i}}{r}_{\text{i}}}{R_{\text{vir,i}}})}{\ln (1+c_{\text{i}})-\frac{c_{\text{i}}}{1+c_{\text{i}}}},$$

$${\vec{a}}_{\text{vis}}=-\sum _{\text{i}}\frac{G{M}_{\text{vir,i}}}{r_{\text{i}}}\frac{\frac{1}{r_{\text{i}}}\ln (1+\frac{c_{\text{i}}{r}_{\text{i}}}{R_{\text{vir,i}}})-\frac{c_{\text{i}}/R_{\text{vir,i}}}{1+c_{\text{i}}{r}_{\text{i}}/R_{\text{vir,i}}}}{\ln (1+c_{\text{i}})-\frac{c_{\text{i}}}{1+c_{\text{i}}}}\widehat{r_{\text{i}}},$$

$$K_{\text{vis}}=\sum _{\text{i}}\sqrt{48}\frac{G{M}_{\text{vir,i}}}{r^3}.$$

Full text

††pubyear: 2017††pagerange: Reconstructing the gravitational field of the local universe–LABEL:LastPage

1 Introduction

General Relativity (GR) has been the reigning paradigm of gravity for almost a century, and yet there is no shortage of alternatives. A range of possible reasons to extend or replace GR has been explored in the literature, including the attempt to overcome theoretical or conceptual difficulties, obviate the need for dark energy or alleviate the cosmological constant problem, account more precisely for astrophysical phenomena, or simply investigate the range of gravity theories consonant with basic physical principles or that may have interesting observational consequences (see e.g. Clifton et al. 2012 and references therein). These diverse motivations have led to a large and heterogeneous parameter space of models which cannot be thoroughly probed by traditional tests. Instead, systematic progress will require the synthesis of evidence from across the range of scales accessible to observation and experiment.

To date, most effort has been devoted to probing gravity in three types of system: the laboratory, the Solar System, and the linear and quasi-linear cosmological regimes. In the laboratory, tests of the equivalence principle (EP) on which GR is premised have now reached $\mathcal{O}(10^{-13})$ precision (Schlamminger et al., 2008), severely limiting deviations from the inverse square law. Within the Solar System, the Parameterised Post-Newtonian (PPN) framework has enabled the coefficients of general metric theories to be constrained at leading relativistic order, providing tight limits on the coupling of extra fields to matter and generic deviations from standard weak-field equations of motion (Nordtvedt, 1968; Will, 1993). In cosmology, GR is combined with the hypotheses of dark matter and dark energy to form the $\Lambda$CDM model, which may be probed via its predictions for Big Bang Nucleosynthesis, the Cosmic Microwave Background, the expansion history of the universe and the growth rate of structure (e.g. Bull et al. 2016 and references therein). Analogously to PPN, frameworks have also been devised to test generic deviations from standard cosmological metrics (e.g. Baker et al. 2013). In none of these cases has a convincing deviation from GR been found.

Nevertheless, these tests do not cover the full range of systems in which deviations from GR may appear, and modified gravity theories may be constructed which satisfy all current experimental bounds and yet exhibit divergent behaviour elsewhere. The possibility that GR may break in some systems – or as a function of certain variables – but not others decouples tests in different regimes, and introduces the possibility that novel gravitational signals may be sequestered in regions of parameter space so far unexplored. Indeed, the notion that physical theories break down at critical values of certain variables, and hence in select systems, is not new: Galilean Relativity gives way to Special at $v$ near $c$, and Newtonian gravity to GR at $\Phi$ near $c^{2}$. Other variables, in which present tests span only a limited range of values, may mark the onset of new gravitational regimes.

The aim of this paper is to contribute to the development of tests of GR in a qualitatively different and relatively under-explored regime: the galactic. In particular, we see the galaxy scale as able to extend the scope of gravitational probes in three main ways:

Many scalar-tensor theories of gravity must employ a “screening mechanism” to suppress their fifth force in the Solar System and laboratory. This may be achieved either by the scalar field acquiring a high mass and hence short range when large (chameleon screening; Khoury & Weltman 2004), or by kinetic (kinetic screening; Babichev et al. 2009) or higher order (Vainshtein screening; Vainshtein 1972) terms in the Lagrangian becoming dominant. In each case the degree of screening is governed by a different function of the scalar field, which simulations have shown to correlate with simple Newtonian proxies: the potential $\Phi$ for chameleon screening, acceleration $\vec{a}$ for kinetic, and some function of the Riemann tensor – we will call it the “curvature” $K$ – for Vainshtein (Cabré et al., 2012; Khoury, 2013). For typical scalar-tensor theories, the transitions between screened and unscreened regions occur at values of these proxies characteristic of galaxies and their environments (Hui et al., 2009; Jain & VanderPlas, 2011); GR is recovered at the larger values that describe laboratory and Solar System tests. Further, for a broad class of chameleon-like theories laboratory fifth force constraints imply a range for the scalar field at cosmological densities of $\lesssim 1$ Mpc, rendering its impact on linear perturbations and the universe’s expansion history negligible (Wang et al., 2012). This makes the intermediate galaxy scale the ideal one at which to test such models. 2. 2.

A dependence of the kinematics of any system on the external gravitational field constitutes a violation of the equivalence principle (EP), and a fortiori of GR. The EP has been tested precisely in the lab (at $a\approx 10\>\text{m}\,\text{s}^{-2}$) and inner Solar System ($a\approx 10^{-2}\>\text{m}\,\text{s}^{-2}$). Galaxies, however, have characteristic accelerations at least six orders of magnitude lower, and hence allow for EP tests in a very different region of parameter space. Beside the case of screening described above (in which EP violation is typically effective at the macroscopic level, rather than fundamental at the level of the action), a paradigmatic example of a theory that exploits this gap to violate the EP is Modified Newtonian Dynamics (MOND; Milgrom 1983a, b, c). This phenomenon is known in MOND as the “external field effect,” whereby a large external acceleration can render a system’s dynamics Newtonian or quasi-Newtonian even for internal accelerations much below the threshold value $a_{0}\approx 10^{-10}\>\text{m}\,\text{s}^{-2}$. As explained further in Section 4.1, galaxies’ mass-to-light ratios and the shapes of their rotation and velocity dispersion profiles would be expected to correlate with external field strength when the latter was sufficiently high, enabling the theory to be tested by means of a prediction complementary to the more conventional ones concerning the internal field alone. Indeed, galaxy environments are among the only with accelerations at and around $a_{0}$, and hence are uniquely capable of testing this and related models. 3. 3.

Baker et al. (2015) mapped out the theoretical and experimental gravity parameter space in terms of $\Phi$ and $K$, and discovered that the curvature values of galaxies are not currently probed by any observational test. These values separate the relatively well-studied small-scale regime at high $K$ from the ‘troubled’ cosmological regime at low $K$, leading Baker et al. (2015) to suggest they may be a natural place for novel gravitational physics to emerge. Since galaxies are the only systems to inhabit this region of parameter space, they are the only ones in principle capable of probing it, regardless of measurement precision.

These motivations suggest we focus on three gravitational variables at the galaxy scale: the Newtonian potential $\Phi$, acceleration $\vec{a}$ and curvature $K$. More generally, these are among the most basic descriptors of the gravitational field, and may therefore be expected a priori to characterise the transitions between gravitational regimes. In the weak-field limit for a set of point masses (as describes galaxies), they scale as $M/r$, $M/r^{2}$ and $M/r^{3}$ respectively, and hence provide a range of relative scalings with the variables $M$ and $r$ that determine an object’s gravitational influence.

The first step towards identifying novel gravitational physics associated with one or more of these variables is to map out their values over the local universe. That is the task of this paper. In particular, we will build maps of $\Phi$, $\vec{a}$ and $K$ out to $\sim 200\>\text{Mpc}$ by combining the contributions of galaxies measured in an all-sky survey, halos in an N-body simulation hosting galaxies too faint to be observed, and a smooth density field not captured by the halo model. Given a set of galaxies with measurements for potential modified gravity signals, these maps will allow the signals’ correlations with each gravitational variable to be determined, and hence the theories above to be investigated and gaps in the experimental parameter space filled in. This will be the subject of future work.

The only previous work along these lines is Cabré et al. (2012), who focus on $\Phi$ as an estimator of the degree of screening in chameleon theories. Besides extending to $\vec{a}$ and $K$, our analysis will build on this pioneering study in several ways: we will use a more complete and homogeneous all-sky source catalogue, apply more sophisticated methods for determining source objects’ mass distributions from their magnitudes, and calibrate and supplement the basic maps with N-body simulations and estimates of the quasi-linear density field respectively. Unlike Cabré et al. (2012), we will make no attempt to relate our proxies to features of specific theories through modified gravity simulations, but focus instead on simply mapping the gravitational variables as robustly as possible. This will allow our results to retain full generality for future application to any model in which these variables are significant.

The structure of this paper is as follows. In Section 2 we describe our method for constructing maps of $\Phi$, $\vec{a}$ and $K$, and in Section 3 we present the results. Section 4 applies our findings to modified gravity, describes the main sources of uncertainty, and discusses the improvements in precision achievable by upcoming galaxy surveys. Section 5 concludes.

2 Method

We make two general methodological comments before detailing our 3-step procedure.

While each of $\Phi$, $\vec{a}$ and $K$ has an “internal” and “external” component, the first due to an object’s own mass and the second that of its environment, we will focus solely on the environmental contribution, which is a function purely of spatial position. The internal component depends on objects’ masses, and must therefore be calculated specifically for a given test galaxy sample. As $\Phi$, $\vec{a}$ and $K$ are our estimators for the degrees of freedom in which GR may receive corrections, we refer to them hereafter as “proxies.”

The Newtonian potential $\Phi$ diverges when calculated in a sphere of $r\rightarrow\infty$, making it necessary to impose a cutoff distance $r_{\text{max}}$. Chameleon-screened theories that motivate our investigation of $\Phi$ supply a natural choice for $r_{\text{max}}$, the Compton wavelength $\lambda_{C}$ of the scalar field that determines the effective range of the fifth force. Masses significantly further apart than $\lambda_{C}$ do not interact via the scalar coupling even if the field is otherwise unscreened. In $f(R)$ gravity, $\lambda_{C}$ is set by the background field value $f_{R0}$, which also determines the potential $\Phi_{\text{crit}}$ at which screening becomes operative (Hu & Sawicki, 2007; Cabré et al., 2012):

[TABLE]

[TABLE]

The $f_{R0}$ scale that may be probed by galaxies is therefore set by typical values of $\Phi$, which is a function of $r_{\text{max}}$. We use a fiducial value of $r_{\text{max}}=10\>\text{Mpc}$, corresponding to $f_{R0}\approx\Phi_{\text{crit}}\approx 10^{-5}$ for $\lambda_{C}=r_{\text{max}}$, but show results also for $r_{\text{max}}=3\>\text{Mpc}$. Since our other proxies ($\vec{a}$ and $K$) fall off more steeply with $r$ (and individual contributions to $\vec{a}$ sum vectorially), sources beyond $3\>\text{Mpc}$ typically contribute little, making them fairly insensitive to the choice of $r_{\text{max}}$ as we demonstrate explicitly in Section 3.3. Their values are normally set by a few nearby objects.

2.1 Primary source catalogue

We base our analysis on the 2M++ galaxy catalogue (Lavaux & Hudson, 2011), a synthesis of 2MASS, 6dF and SDSS data. This is an optimum catalogue for our purposes, for three reasons. First, it was designed with the goal of high all-sky completeness out to 200 Mpc. Complete sky coverage will prevent our having to restrict our maps to the footprint of a single survey, and 200 Mpc is around the largest distance that potential modified gravity signals are robustly measurable at present (Vikram et al., 2013). Deeper surveys (e.g. the SDSS main sample) exist only over part of the sky. Second, it has a homogeneous limiting $K$-band magnitude ($m_{K}<12.5$), which will facilitate our modelling of the contributions to the proxies from missing objects by means of an N-body simulation in Section 2.2. Finally, the catalogue has already been used to estimate the smoothed density field in the survey volume (Lavaux & Jasche, 2016), which provides an important contribution to our proxies that is largely independent of those above.

To determine the values of $\Phi$, $\vec{a}$ and $K$ sourced by the 2M++ galaxies, we first estimate the mass distributions of their dark matter halos. We utilise the technique of “halo abundance matching” (AM), which maps galaxies to halos produced in an N-body simulation by postulating a nearly monotonic relationship between luminosity and some function of halo virial mass and concentration (Kravtsov et al., 2004; Conroy et al., 2006; Moster et al., 2010). For a suitable scatter and function of halo properties, this has been shown to yield excellent agreement with galaxy clustering statistics (e.g. Reddick et al. 2013; Lehmann et al. 2015) and moderate agreement with galaxies’ internal dynamics (e.g. Trujillo-Gomez et al. 2011; Desmond & Wechsler 2015, 2016; Desmond 2017). We use the 2M++ $K$-band luminosity function measured in Lavaux & Hudson (2011), and rockstar (Behroozi et al., 2013) halos from the darksky-400 simulation (Skillman et al., 2014), a $(400\>\text{Mpc/h})^{3}$ box with $4096^{3}$ particles run with the 2hot code (Warren, 2013). Our specific AM model will be that of Lehmann et al. (2015), which matches luminosity to $v_{\text{vir}}(v_{\text{max}}/v_{\text{vir}})^{\alpha}$, with best-fitting values $\alpha=0.6$ and uniform Gaussian scatter $\sigma_{\text{AM}}=0.16$ dex. We caution that these parameters were derived using an $r$-band luminosity function from SDSS rather than a $K$-band one from 2M++, although we have verified that a basic counts-in-cells clustering statistic out to $10$ Mpc is consistent between the 2M++ dataset and our AM catalogue. This method enables us to generate an absolute $K$-band magnitude for each halo in the simulation box.111In principle there is further information in the group statistics of the 2M++ catalogue, for two reasons. First, using the group luminosity function rather than that of individual galaxies may improve the estimate of the halo masses derived from AM at the bright end. Second, the velocity dispersions of galaxies in a group provide complementary information to the luminosity on the distribution of dark matter mass. Folding this in would improve the precision of the results, although it is beyond the present scope of our work.

Next, we next calculate the absolute magnitude $M$ of each galaxy in the 2M++ catalogue:

[TABLE]

where $m$ is the apparent magnitude, $A_{k}(l,b)$ describes dust absorption in the direction $(l,b)$, $k(z)$ is the k-correction due to redshifting of the spectrum, $e(z)$ is the correction for stellar population evolution, and $D_{L}$ is the luminosity distance. Following Lavaux & Hudson (2011), we take $k(z)=-2.1z$, $e(z)=0.8z$ and $A_{k}(l,b)=0.35E_{\text{(B-V)}}(l,b)$, and calculate the extinction factor $E_{(B-V)}$ from the dust map of Schlegel et al. (1998). We associate each 2M++ galaxy with the halo in the simulation closest to it in absolute magnitude. We assume this halo to have a Navarro-Frenk-White (NFW) density profile (Navarro et al., 1996),

[TABLE]

with virial mass $M_{\text{vir}}$, virial radius $R_{\text{vir}}$ and concentration $c\equiv R_{\text{vir}}/r_{\text{s}}$ as output by rockstar.

With this mass distribution in place, we are in position to calculate the proxy values sourced by each 2M++ object. For $\Phi$ and $\vec{a}$ we use standard forms for the NFW profile (Cole & Lacey, 1996):

[TABLE]

[TABLE]

for source object $i$ at distance $r_{\text{i}}$ from the test point, where $\hat{r}$ points from the source halo centre to the test point. The subscript ‘vis’ denotes that these contributions to the proxies derive from objects visible to the 2M++ survey; further contributions will be described in Sections 2.2 and 2.3. We measure curvature using the Kretschmann scalar $K\equiv(R^{\alpha\beta\gamma\delta}R_{\alpha\beta\gamma\delta})^{1/2}$, both because it is non-zero in vacuum and to promote compatibility with Baker et al. (2015), who find it to be of use in synthesising laboratory, astrophysical and cosmological constraints on gravity. While in detail the Kretschmann scalar for a set of point masses is nonlinear, we find the error incurred by approximating it as the linear sum of individual contributions to be small. We assume also that $K$ can be adequately calculated by treating each source halo as a point mass at its centre. As we make no attempt to estimate the accuracy of this approximation, $K$ ought not to be considered more than an order-of-magnitude estimator of the true curvature.

[TABLE]

2.2 Restoring missing halos

The 2M++ catalogue clearly does not include all mass within 200 Mpc, and hence the calculation of Section 2.1 underestimates $\Phi$, $\vec{a}$ and $K$. Part of the missing mass resides in halos hosting galaxies too faint to be observed in the 2M++ survey yet sufficiently massive to be well-resolved in the N-body simulation to which we have already had recourse. This mass can be filled in by calibrating with that simulation as follows:

Define a well-resolved halo as one containing at least $1000$ particles within $R_{\text{vir}}$, which for the darksky-400 simulation equates to a mass cut $M_{\text{vir}}>7.63\times 10^{10}\>h^{-1}\>M_{\odot}$. Above this threshold halos have reliable concentration measurements (Diemer & Kravtsov, 2015), allowing Eqs. 5-7 to be evaluated accurately, and the halo mass function is complete (Reed et al., 2013), ensuring that all such halos are identified by the halo finder. 2. 2.

Using the absolute magnitude assigned to each halo from AM, and assuming the observer to be at the centre of the box, determine the apparent magnitude from Eq. 3.222Although the N-body box is not intended to represent the actual local universe, we nevertheless use the angular dependence of the Schlegel et al. (1998) dust map to calculate the extinction at the position of each halo. As this term is small, a different angular dependence would have little effect. Thus split the halos into those that would be visible in the 2M++ catalogue ($m<12.5$) and those that would not ($m>12.5$). 3. 3.

At the position of each well-resolved halo, calculate $\Phi$ from all other well-resolved halos within 10 Mpc ($\Phi_{\text{hal}}$), and from the visible ones alone ($\Phi_{\text{vis}}$). (It is advisable to use the positions of the halos themselves rather than random points in the box because the latter would not reproduce the clustering properties of real galaxies. In particular, low density regions would be much more prominently represented.) Take $c_{\Phi}\equiv\Phi_{\text{hal}}/\Phi_{\text{vis}}$ as the multiplicative “correction factor” required to map the latter potential onto the former.333Note that due to redshift incompleteness in the 2M++ catalogue, as discussed in Lavaux & Hudson (2011), fewer objects are actually observed than pass the magnitude cut, an effect more pronounced in some directions than others. This causes our correction factors to be underestimated, as we assume that all objects with $m<12.5$ contribute to the observable proxies. However, as the survey achieved high completeness in all regions above the galactic disk (Lavaux & Hudson, 2011), this effect is insignificant there relative to the other sources of uncertainty in our analysis. Repeat to obtain the corresponding correction factors for acceleration and curvature, $c_{a}\equiv|\vec{a}_{\text{hal}}|/|\vec{a}_{\text{vis}}|$ (describing the change in the magnitude of $\vec{a}$), $c_{\theta}\equiv\frac{\vec{a}_{\text{hal}}\cdot\vec{a}_{\text{vis}}}{|\vec{a}_{\text{hal}}||\vec{a}_{\text{vis}}|}$ (describing the change in the direction of $\vec{a}$) and $c_{K}\equiv K_{\text{hal}}/K_{\text{vis}}$. Denote the set of correction factors by $\vec{c}\equiv\{c_{\Phi},c_{a},c_{\theta},c_{K}\}$. Record also the distance to the halo in question ($d$), the number of visible objects within 10 Mpc of it ($N_{10}$), and $\Phi_{\text{vis}}$. We find these observables to correlate strongly with $\vec{c}$, and hence will use them to predict it. 4. 4.

Partition the {$c_{\Phi}$, $c_{a}$, $c_{\theta}$, $c_{K}$, $d$, $N_{10}$, $\Phi_{\text{vis}}$} space into discrete cells and populate it with the halos in the box. This estimates their joint frequency distribution. We sample logarithmically in $c_{K}$ and $\Phi_{\text{vis}}$, take $35$ bins in each direction, and choose the upper and lower limits to enclose at least $95$ per cent of the simulation points. 5. 5.

For each test galaxy in the real universe, estimate $\vec{c}$ by drawing random values for each of its components from this distribution at the {$d$, $N_{10}$, $\Phi_{\text{vis}}$} of that galaxy. This measures $\vec{c}$ conditional on {$d$, $N_{10}$, $\Phi_{\text{vis}}$}, and retains the full covariance among its components. 6. 6.

Estimate the total potential and curvature values for the test galaxy by multiplying the estimates from the visible mass alone by the corresponding correction factors, i.e. $\Phi_{\text{hal}}=c_{\phi}\times\Phi_{\text{vis}}$ and $K_{\text{hal}}=c_{K}\times K_{\text{vis}}$. For acceleration, first multiply $|\vec{a}_{\text{vis}}|$ by $c_{a}$, then randomly rotate the resulting vector through $c_{\theta}$.

2.3 Modelling the remaining missing mass

The combination of halos hosting visible galaxies and those well-resolved in our simulation accounts for only $\sim 10$ per cent of the universe’s mass; the remainder is located either in smaller halos or outside any virialised structure (e.g. van Daalen & Schaye 2015). The analysis of Section 2.2 therefore continues to underestimate the proxies. To model the remaining mass we use the results of Lavaux & Jasche (2016), who apply 2LPT in a Bayesian framework to derive the probability distribution for the $z=0$ density field consistent with the number density and peculiar velocities of the 2M++ galaxies. This is done on a 3D grid of spacing $2.3\>h^{-1}\>\text{Mpc}$ and therefore provides an accurate estimate of the power spectrum on scales $k\lesssim 1\>h\>\text{Mpc}^{-1}$.

Although this method reconstructs approximately $100$ per cent of the mass in the 2M++ volume, the smoothing scale of $\sim 3.4$ Mpc means that the screening proxies cannot be estimated from this mass alone; their scaling with $r^{-n}$ makes them sensitive to smaller-scale overdensities such as we have already modelled in Sections 2.1 and 2.2. We estimate the proxy values due to all mass by summing the contributions from the resolved halos and 2M++ galaxies, derived in Section 2.2, to those from the mass field derived by smoothly interpolating the density grid of Lavaux & Jasche (2016) at a representative point in their MCMC chain.444The results are not significantly altered by choosing a different high-likelihood point from the chain, although see Section 4.2 for further discussion of this source of uncertainty. As this method double-counts a small fraction of the mass, it systematically overestimates $\Phi$, $a$ and $K$;555This provides further motivation for limiting the contribution in Section 2.2 to halos with $>1000$ particles: including smaller halos (or using the particles themselves) would double-count mass more egregiously. our results for the minimum degree of screening – and hence the strongest limits on scalar-tensor theory screening thresholds afforded in principle by galaxy observations – are therefore conservative. The two contributions do however provide power on largely complementary scales: the halos may be seen as filling in the power from the 1-halo term that is washed out by the $3.4$ Mpc smoothing. We provide further discussion of this difficulty, including means by which it may be overcome, in Section 4.2.

Our code for performing all the operations of this section is publicly available at https://www2.physics.ox.ac.uk/contacts/people/desmond.

3 Results

In a real application, our method would be used to evaluate $\Phi$, $\vec{a}$ and $K$ at the positions of galaxies with measured values for potential modified gravity signals. Since we do not yet possess such a sample, we instead show results at the positions of galaxies in our source catalogue, which at least sample realistic galaxy environments. In this section we will investigate the impact of unseen halo mass and the reliability of our method of estimating it (Section 3.1), the relative contributions of the three mass components (Section 3.2), and the overall proxy distributions (Section 3.3).

3.1 Modelling unseen mass

We begin with our determination of the correction factors $c_{\Phi}$, $c_{a}$, $c_{\theta}$ and $c_{K}$: histograms of these quantities are shown in Figure 1 for galaxies within 200 and 50 Mpc. The averages and standard deviations of these distributions are reported in the second and third columns of Table 1, while the fourth contains the average dispersions at fixed values of $d$, $N_{10}$ and $\Phi_{\text{vis}}$. These measure the uncertainties on our estimate of $\vec{c}$ using the method of Section 2.2. The most significant features of these results are the following:

As the correction factors are not negligible, they are important in any calculation of $\Phi$, $\vec{a}$ or $K$. Tails towards high values indicate environments dominated by unseen mass. 2. 2.

On average, $c_{K}>c_{a}>c_{\Phi}$, and the variances of the distributions increase also in that order (this is due to contributions below the relative frequency cutoff in Fig. 1 and hence not visible there). This is the order of decreasing sensitivity of the proxy to the separation $r$ of source and test point: $|\vec{a}|$ and $K$ are more dependent than $\Phi$ on nearby low-mass objects which are less likely to be included in 2M++. This makes our method most reliable for determining $\Phi$ and least for $K$, as the latter is more likely to be dominated by a single faint object. Indeed, at some test points $K_{\text{hal}}$ exceeds $K_{\text{vis}}$ by more than an order of magnitude (note the logarithm in the final row of Table 1), making a method based on a magnitude-limited source galaxy sample unreliable. 3. 3.

The widths of the $\vec{c}$ distributions are significantly smaller after calibrating with {$d$, $N_{10}$, $\Phi_{\text{vis}}$} than before, showing that this method improves the precision of our proxy calculations. The remaining variance contributes to the final uncertainty in our maps, as we discuss further in Section 4.2. 4. 4.

The correction factors are larger further away since a smaller fraction of objects passes the 2M++ magnitude cut.

3.2 Relative contributions to the proxies

Next we show in Figure 2 the relative contributions to the proxies from the observed 2M++ galaxies, the halos from the N-body box and the smoothed density field. That they are roughly the same order of magnitude shows that none are negligible. In the case of $\Phi$, the fractional contribution of the 2M++ galaxies is roughly uniform up to 80 per cent, indicating that in few cases is a fraction of $\Phi$ greater than this sourced by those galaxies alone. The invisible halos in the simulation typically contribute $\sim 20\%$, and rarely more than $50\%$. Some test galaxies have no or very halos or 2M++ objects within 10 Mpc, making both $\Phi_{\text{vis}}$ and $\Phi_{\text{hal}}$ very small; the potential for these objects effectively derives purely from the smoothed density field.

Interpretation of the relative contributions of the three components to $\vec{a}$ (Fig. 2) is complicated by the fact that they do not sum directly at a given test point, but rather vectorially. Depending on the phases, this means that the magnitudes of the acceleration due to any one of the components may be larger than the total acceleration, making the fractional contribution greater than 1. The 2M++ part peaks at around $75\%$, and most often the remaining $25\%$ (with roughly aligned directions) comes from unseen halos. Since mass in the low-frequency modes of the smoothed density field is by definition fairly homogeneous (and that in the zero-frequency mode entirely so), it tends to have little effect on the acceleration as contributions from opposite sides of the test galaxy roughly cancel out. That said, the small fraction of test points with little 2M++ or unseen halo mass within 10 Mpc have accelerations dominated by the low-frequency modes, as evidenced by the peak in the green curve at $\sim 1$.

We find the curvature typically to be dominated by a single component, most often the 2M++ galaxies (Fig. 2). Since $K$ falls more steeply with $r$ than $\Phi$ or $\vec{a}$, and hence samples on average a smaller volume around a given test point, it is less affected by mass in the smoothed density field.

3.3 Distributions of gravitational variables

Figure 3 shows the distributions of $\Phi$, $\vec{a}$ and $K$ across our test objects, including source objects up to $r_{\text{max}}=10\>\text{Mpc}$ (blue) or 3 Mpc (red). As mentioned in Section 1, the choice of $r_{\text{max}}$ is more important for $\Phi$ than $a$ or $K$ due to its weaker fall off with the separation of test and source object: most of the mass contributing to $a$ and $K$ is located within 3 Mpc. We find the proxy distributions not to depend significantly on the masses of the test galaxies, and thus they provide an indication of the range of values probed by any galaxy sample without explicit selection on environment.

Finally, in Figure 4 we display the proxy values across a $300\>\text{Mpc}\times 300\>\text{Mpc}$ slice of the local universe, again using $r_{\text{max}}=10$ Mpc. The Milky Way is located at $x=y=0$, and the sampling resolution is 1.5 Mpc. The increasingly detailed level of structure in these maps when moving from $\Phi$ to $a$ to $K$ reflects again the greater sensitivity of the proxy to the separation $r$ of source and test points. Their approximate uniformity over the entire area attests to the accuracy of our method for restoring missing halos, of which there are more further away.

These maps reveal promising regions in which to search for modified gravity; for example, regions of deep blue in Fig. 4(a) have very low $\Phi/c^{2}$ (down to few $\times 10^{-7}$), allowing any galaxies within them to remain unscreened in chameleon theories with small background scalar field values. Conversely, regions of deep red in Fig. 4(b) have particularly high accelerations, making them most likely to harbour external field-dominated galaxies in MOND.

4 Discussion

In this section we discuss the implications of our results for galaxy-scale tests of modified gravity (Section 4.1), and then describe the major uncertainties in our maps and the potential for future improvement (Section 4.2).

4.1 Applications to modified gravity

We begin with a few general inferences about the usefulness of galaxy-scale tests of gravity from the distributions of our proxy values (Fig. 3).

In chameleon scalar-tensor theories, the value of $\Phi$ marking the transition between the screened and unscreened regimes is set by the background value of the scalar field (Hu & Sawicki, 2007). A test of screening involves a comparison between putatively screened and unscreened samples, and hence the lowest value of $\Phi$ determines the strength of the constraint that may in principle be placed on the theory. Figure 3 suggests that a significant number of galaxies are likely to be environmentally unscreened for $\Phi_{\text{crit}}/c^{2}=10^{-7}$ and $r_{\text{max}}=\lambda_{C}=1\>\text{Mpc}$, giving galaxy-scale tests the potential to probe $f_{R0}$ to at least the $10^{-7}$ level (see Eqs. 1 and 2). Although our results differ in detail from those of Cabré et al. (2012), due to our more sophisticated treatment of missing mass, we share this optimistic conclusion. 2. 2.

Typical external accelerations are approximately a few $\times 10^{-15}\>\text{km}\,\text{s}^{-2}$, with a range of $\sim 10^{-16}-{\rm few}\times 10^{-14}$ (Fig. 3). Given the relatively small difference between the red and blue curves in that figure, especially at the high end, this is not likely to change much if $r_{\text{max}}$ was extended. In MOND-type models, a galaxy enters the external field-dominated regime when the external acceleration $a_{\text{ex}}$, due to surrounding mass, exceeds the acceleration $a_{\text{in}}$ generated by the galaxy itself, provided $a_{\text{in}}$ is less than the characteristic acceleration $a_{0}\approx 1.2\times 10^{-13}\>\text{km}\,\text{s}^{-2}$. Thus our distribution of $a_{\text{ex}}$ indicates the maximum internal accelerations – and hence surface densities – that galaxies can have for the external field to be important. A typical value $a_{\text{ex}}\sim{\rm few}\times 10^{-15}\>\text{km}\,\text{s}^{-2}$ would correspond to surface density $\Sigma\approx 20\>M_{\odot}\>\text{pc}^{-2}$, which is very small even for a low surface brightness galaxy. Thus only unusually high $a_{\text{ex}}$ values would appreciably impact dynamics. Galaxies with $a_{\text{ex}}>a_{0}$ (of which Fig. 3 suggests there are very few) should exhibit fully Newtonian behaviour. We caution however that we have measured the proxies only at the positions of the 2M++ galaxies; more promising galaxies for tests of the external field effect are dwarfs very close to a single massive object which provides the bulk of the external acceleration (such as the dSphs of the Milky Way and M31; McGaugh & Milgrom 2013a, b). It is unlikely that these would be included among our test galaxies due to high apparent magnitudes at any distance beyond the local group. At the positions of such galaxies located by other means, however, our map would be expected to report $a_{\text{ex}}$ accurately provided the galaxies sourcing the field were included in the 2M++ catalogue.

At the other end of the scale, by spherical symmetry (and unlike for $\Phi$ and $K$), the smoothed density field does not provide an irreducible background value of $\vec{a}$, and hence there is no fundamental reason why regions of arbitrarily low $\vec{a}$ could not exist. This enables screening mechanisms characterised by acceleration (e.g. kinetic) to be tested in principle to arbitrary precision in galaxies, provided a sample in sufficiently remote environments could be compiled. 3. 3.

The external curvature is typically around $K=\text{few}\>\times 10^{-54}\>\text{cm}^{-2}$, with a range of roughly $10^{-56}-10^{-51}$. This sets the precision with which Vainshtein screening may be tested at the galaxy scale. These values are at the low end of the “curvature desert” identified in Baker et al. (2015), enabling galaxies potentially to inhabit any part of this region when their internal curvature contributions are added.666Baker et al. (2015) consider only the internal contributions. For models in which the important variable is the total proxy value, our results should therefore be considered summed to the galaxy lines in their fig. 1. As $K$ falls off rapidly with $r$, these contributions will typically greatly exceed the environmental ones considered here. Note also that the lower limit of our distribution, $K\approx 10^{-56}\>\text{cm}^{-2}$, coincides roughly with the “Lambda” line of Baker et al. (2015)’s fig. 1, which is approximately the curvature of the cosmological background.

More detailed or model-specific conclusions concerning modified gravity will require a test galaxy sample for which measurements of potential modified gravity signals have independently been made. Correlations of these signals with the proxies determined from our maps (plus the internal contributions from the galaxies themselves) may then be used to search for a transition to a new gravitational regime or test specific theories for which these proxies are important.

Such signals will depend on the theory in question. In scalar-tensor theories, the stars in unscreened galaxies may self-screen, leading to differences between their kinematics and those of the gas and dark matter. In particular, the stars would fall more slowly in an external field, leading to displacements between the centroids of stellar and Hi disks preferentially aligned with the external acceleration, and to an enhancement of Hi over H$\alpha$ rotation curves. The offset between halo and disk centres may then cause stellar disks to warp, and both the photometry and kinematics of the disk to develop asymmetries (Jain & VanderPlas, 2011; Vikram et al., 2013).

In MOND, the value of $|\vec{a}|$ sets the ratio of a system’s dynamical to baryonic mass, and also the shape of the rotation or velocity dispersion profile. In galaxies with internal acceleration $a_{\text{in}}\ll a_{0}$, rotation curves would be flat for $a_{\text{ex}}<a_{\text{in}}$ (the standard deep-MOND regime), but Keplerian for $a_{\text{ex}}>a_{\text{in}}$ (the external field-dominated regime; Famaey & McGaugh 2012). The galaxy’s total mass discrepancy, measured by $\mathcal{D}\equiv\frac{RV^{2}}{M_{\text{bar}}G}$, would then be $\sim a_{0}/a_{\text{ex}}$ for $a_{\text{ex}}<a_{0}$ (the quasi-Newtonian regime) and $1$ for $a_{\text{ex}}>a_{0}$ (the Newtonian limit). The correlation of $\mathcal{D}$ with $a_{\text{ex}}$ – the external-field analogue of the mass discrepancy–acceleration relation (McGaugh, 1999; Lelli et al., 2017) – would have a characteristic shape and a feature at $a_{\text{ex}}=a_{0}$.

Finally, we note that our maps may also find application in the study of galaxy formation. While correlations between $\Phi$, $a$ or $K$ and galaxy signals may be evidence for novel gravitational physics, they would more likely attest to galaxy growth or evolution effects in $\Lambda$CDM. They may therefore be useful in studies of quenching, conformity, satellite–host interactions and the environment dependence of the galaxy–halo connection. These phenomena will also constitute important systematics in tests of gravity.

4.2 Error analysis and future prospects

4.2.1 Sources of uncertainty

Here we describe the most significant uncertainties in our maps, in roughly decreasing order of importance. This will give an indication of the areas on which future analyses should focus to best improve precision, and the ways in which upcoming galaxy surveys will be able to contribute. The parameters that form the input to our analysis, and the distributions from which they are drawn, are summarised in Table 2.

We estimated the correction factors required to convert the proxies sourced by visible mass to values due to all resolved halos on a point-by-point basis using their correlations with distance, the number of 2M++ objects within 10 Mpc, and the potential from these objects. However, as shown in Table 1, significant uncertainties remain, especially for $|\vec{a}|$ ($\sim 200\%$) and $K$ ($\sim 0.7$ dex). These uncertainties increase with distance. Since this uncertainty derives from the fraction of objects included in the galaxy survey, it is best reduced by employing a deeper survey. Over a third of the sky the magnitude limit may already be improved by using the SDSS main sample, with a Petrosian $r$-band magnitude limit of $17.7$, and the Taipan survey will soon provide similar data in the South. Within the next decade, the Dark Energy Spectroscopic Instrument (DESI) will push around two magnitudes deeper, directly locating an even greater proportion of the local universe’s mass. Indeed, we find that for a DESI-like catalogue with a limiting $r$-band magnitude of $19.7$ and average $r-K$ colour of $2.7$ (e.g. Bell et al. 2003), essentially all halos with $>1000$ particles in our simulation would be expected to host visible galaxies out to 200 Mpc, all but eliminating the need for correction factors. We refer the reader to Jain (2011), Jain & Khoury (2010) and Jain et al. (2013) for further discussion of survey requirements and optimisation for the purposes of gravitational physics. 2. 2.

The analysis of Lavaux & Jasche (2016) derives a probability distribution for the smoothed density field of the local universe from the positions and velocities of the 2M++ galaxies. The uncertainty that the width of this distribution propagates into our maps may be estimated by averaging over the results from all field configurations in the posterior, effectively marginalising over the possibilities. Although we leave this to future work, we have verified by direct substitution of several high-likelihood configurations that at typical locations the resulting uncertainties in $\Phi$, $a$ and $K$ are subdominant to those described above. We caution however that this uncertainty is most important in void regions containing few massive halos, which are the best places to search for screening. 3. 3.

Galaxies of given luminosity may reside in halos with a range of properties, so that the observed distribution of light does not uniquely determine the distribution of total mass. Thus the 2M++ galaxy sample corresponds to a range of possible $\Phi_{\text{vis}}$, $\vec{a}_{\text{vis}}$ and $K_{\text{vis}}$ fields, of which we show simply example realisations. Within the AM framework, this derives from the combination of a statistical uncertainty due to the range of possible halo $M_{\text{vir}}$ and $c$ for a galaxy of given luminosity, quantified by $\sigma_{\text{AM}}$, and a potential systematic error due to uncertainties in the model parameters $\alpha$ and $\sigma_{\text{AM}}$ themselves. Approximately, $\sigma_{\text{AM}}=0.16$ (Lehmann et al., 2015) produces a $\sim 0.16$ dex uncertainty in the mass of the halo of a galaxy of given luminosity, and hence in the values of $\Phi$, $a$ and $K$ sourced by that halo. This uncertainty is most significant for test points surrounded by many 2M++ galaxies. In addition, baryonic effects may cause halo profiles to differ from the NFW form that we assume (e.g. Blumenthal et al. 1986; Governato et al. 2010; Pontzen & Governato 2012). 4. 4.

Uncertainties in the measured magnitude, position and distance of each source object propagate into uncertainties in our maps.

These uncertainties may all be accounted for in tandem by generating many Monte Carlo realisations of our maps with each input parameter drawn randomly from its set of allowed values, thus producing a posterior distribution of $\Phi$, $a$ and $K$ at each test point. These may be further propagated into constraints on modified gravity parameters given a full inference framework and a test galaxy sample. This will be the subject of future work.

4.2.2 Overestimation of total mass

As discussed in Section 2.3, including halos in addition to the density field of Lavaux & Jasche (2016) overestimates the total mass in the universe by $\sim 10$ per cent, and hence $\Phi$, $a$ and $K$ by a comparable factor. Here we describe two methods by which this discrepancy may be eliminated.

First, estimate for each test point the total mass within the surrounding volume of interest (in our fiducial case a sphere of radius $10$ Mpc) using a smoothed field akin to Lavaux & Jasche (2016). This is reasonable provided the smoothing scale is less than the size of the volume of interest. Then estimate the true configuration of that mass as follows. A part is associated with visible galaxies (e.g. in 2M++), and may therefore be situated more precisely in halos around those galaxies specified by inverse AM, as in Section 2.1. Another part is in mass elements captured by the reference N-body simulation. As in Section 2.2, the distribution of this mass may be estimated by examining regions of the simulation with similar distributions of visible mass, although in this case all the particles in the simulation may be used rather than just well-resolved halos because there is no danger of double-counting mass. This provides a plausible small-scale configuration of the smoothed mass (and hence the corresponding $\Phi$, $a$ and $K$ at the test point), and the possible configurations may then be marginalised over to estimate the uncertainty in the proxy values. Any remaining mass must be located in structures below the resolution limit of the simulation, and the maximum-entropy assumption for this mass is that it is smoothly distributed, as in the original field. This procedure has the further advantage of reconstructing modes of the power spectrum below the $2.4\>h^{-1}$ Mpc limit of Lavaux & Jasche (2016) field alone (see their figure 2), enhancing the utility of that framework. 2. 2.

It is necessary to add by hand the small-scale power due to halos because the 2M++ density field cannot be resolved below $2.4\>h^{-1}$ Mpc. These halos come from a simulation box bearing no relation to the local universe, necessitating the correlation with observables described in Section 2.2. However, the analysis of Lavaux & Jasche (2016) also infers the initial conditions of structure formation in the 2M++ volume. This opens the possibility of running simulations constrained to match the smoothed density field at $z=0$. Such simulations would automatically place the correct mass in the simulation volume, as well as producing realisations of the density field resolved to the level of the particle mass. The distributions of $\Phi$, $a$ and $K$ maps would then be obtained by running multiple simulations consistent with the 2M++ number density and peculiar velocity fields, thereby marginalising over both the smoothed $z=0$ density field and the mass distribution on smaller scales. This method has the potential to obviate the need for all the inputs of Table 2 beside the second.

4.2.3 Further considerations for tests of gravity

A range of possible $\Phi$, $\vec{a}$ and $K$ fields corresponds to an uncertainty in the degree of screening of a given galaxy (under scalar–tensor gravity) or in its effective force law (under MOND). This will introduce marginal cases when dividing a test galaxy sample into different gravitational regimes, weakening null tests of GR and constraints on modified gravity parameters. The magnitude of the total uncertainty therefore determines the combination of minimum signal strength and sample size that would be statistically significant.

A full comparison with modified gravity expectations would have further uncertainties. To make predictions for a specific theory, it would likely be invalid in detail to assume $\Lambda$CDM to calculate the halo mass function, the galaxy–halo connection, the relation between absolute and apparent magnitude, or the quasi-linear modes of the 2M++ field. We justify our doing so by arguing that the types of modified gravity theory that motivate our work are likely to generate larger deviations from GR on galactic than cosmological scales. The impact on cosmology – and hence those parts of our analysis that rely on $\Lambda$CDM – would be a second order effect in the prediction of the final signals.

This argument is a priori least plausible in the case of MOND, which is unique among the models motivating our work in that it seeks to eliminate dark matter. MOND is not connected continuously to GR and $\Lambda$CDM as scalar-tensor theories are, but rather postulates a discrete deviation at the level of the Newtonian limit. The level of systematic error between our calculated $\vec{a}$ and that in MOND likely depends on the specific formulation of the theory. In QUMOND (Milgrom, 2010), the MOND acceleration is simply an algebraic function of the Newtonian acceleration due to baryons alone, and the difference may therefore be represented by a distribution of “phantom dark matter.” Although there is no a priori reason why this distribution should be correctly given by AM (for one there is no longer a halo mass function), the fact that observational estimates of the relation between visible and dynamical mass roughly agree with the AM prediction (Behroozi et al., 2010) suggests that this technique would give a sensible estimate for the amount of phantom dark matter as a function of luminosity in the mass range of interest. We would however expect $c_{a}$ to be nearer 1 in this case, as there would be no contribution to $\vec{a}$ from mass unassociated with light. As this tends to increase $|\vec{a}|$, our result would constitute an upper limit on the external acceleration field in MOND, making our conclusion that there are very few objects in the fully external field-dominated regime robust.

Other formulations, however, include a curl field in the relation between Newtonian and MONDian acceleration, which introduces further nonlinearities into the dependence of the total acceleration on the (baryonic) mass distribution, in addition to a potential offset between the direction of the true acceleration and the one we have calculated (e.g. Llinares et al. 2008; Famaey & McGaugh 2012). In such cases our method may give a poor approximation to the acceleration governing the MOND external field effect; we leave further investigation of this effect to future work.

5 Conclusion

We have combined an all-sky galaxy catalogue, a high-resolution N-body box and an estimate of the smoothed density field of the local universe to map the gravitational potential $\Phi$, acceleration $\vec{a}$ and curvature $K$ to a distance of 200 Mpc. As the fundamental variables of the gravitational field, these may be expected on general grounds to mark the transitions between gravitational regimes, and at the galaxy scale they are of specific importance in screened scalar-tensor theories and EP-violating models such as MOND. We use our maps to determine the sensitivity that galaxy-scale tests of modified gravity may be expected to achieve, and identify promising regions of the local universe. Our method may be easily adapted to any source and test galaxy sample, and we make public the code for doing so.

From the overall distributions of the gravitational variables we draw four general conclusions.

• •

When calculated in spheres of radius 10 (3) Mpc (typical Compton wavelengths of a light scalar field), the distribution of $\Phi/c^{2}$ peaks at around $10^{-5}$ ($10^{-6}$), and has a minimum of $10^{-6}$ (few $\times 10^{-8}$) set largely by the smooth component of the density field. Since the background scalar field in chameleon-screened scalar-tensor theories (e.g. $f_{R0}$ in $f(R)$ gravity) is of order the threshold in $\Phi/c^{2}$ at which screening becomes operative, astrophysical measurements are capable of constraining this to at least the $10^{-7}$ level by investigating galaxies with low internal $\Phi$ values. This makes galaxy-scale tests of screening significantly stronger than cosmological probes (e.g. Lombriser et al. 2012).

• •

The external acceleration field is typically around a few $\times 10^{-15}\>\text{km}\,\text{s}^{-2}$ and rarely exceeds the characteristic acceleration $a_{0}\approx 10^{-13}\>\text{km}\,\text{s}^{-2}$ at which visible and dynamical masses diverge. Thus, in the context of MOND-type models, few objects will generically be in the fully external field-dominated regime in which dynamics is Newtonian even for arbitrarily small internal acceleration. Locating such objects to test the external field effect would therefore require targeted searches around known mass concentrations.

• •

Curvature values are centred around $10^{-55}\>\text{cm}^{-2}$, at the low end of the “curvature desert” that Baker et al. (2015) identifies as a potentially interesting region for modified gravity.

• •

The precision with which $\Phi$, $\vec{a}$ and $K$ may be mapped is limited by our ability to detect faint objects at large distance, and galaxy catalogues from future surveys such as DESI and LSST will therefore greatly reduce uncertainties. Theoretical advances such as the creation of simulations constrained to match the local volume may also be used to improve precision.

In future work we will improve our total mass estimation, model fully the uncertainties in our maps by Monte Carlo marginalisation, and apply our framework to specific tests of gravity in galaxies.

Acknowledgements

We thank Julien Devriendt, Bhuvnesh Jain, Vinu Vikram and Risa Wechsler for helpful discussions during the course of this work, and Tessa Baker, Gary Mamon, Stacy McGaugh and an anonymous referee for comments on the manuscript. HD received support from the Balzan foundation via New College, Oxford, the U.S. Department of Energy under contract number DE-AC02-76SF00515, and St John’s College, Oxford. PGF is supported by ERC H2020 693024 GravityLS project, the Beecroft Trust and STFC. GL acknowledges financial support from “Programme National de Cosmologie and Galaxies” (PNCG) of CNRS/INSU, France. GL is supported by ANR grant ANR-16-CE23-0002.

This work made use of the Dark Sky simulations, which were produced using an INCITE 2014 allocation on the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory. We thank the Dark Sky collaboration for creating and providing access to this simulation, and Sam Skillman and Yao-Yuan Mao for running rockstar and consistent trees on it, respectively. Additional computation was performed at SLAC National Accelerator Laboratory.

This research was supported by the DFG cluster of excellence “Origin and Structure of the Universe” (www.universe-cluster.de). This work was partly made in the ILP LABEX (under reference ANR-10-LABX-63) supported by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02.