Screening maps of the local Universe I -- Methodology

TL;DR

This paper presents a new methodology for creating detailed screening maps of the local Universe to identify regions where deviations from General Relativity are suppressed, using Bayesian mass reconstruction and applying it to $f(R)$ gravity models.

Contribution

It introduces a novel approach combining Bayesian mass reconstruction with modified gravity equations to improve predictions of screening effects in the local Universe.

Findings

First qualitative screening results for $f(R)$ gravity.

Demonstration of the methodology's potential for future observational tests.

Framework applicable to various modified gravity models.

Abstract

We introduce the LOCal Universe Screening Test Suite (LOCUSTS) project, an effort to create 'screening maps' in the nearby Universe to identify regions in our neighbourhood which are screened, i.e., regions where deviations from General Relativity (GR) are suppressed, in various modified gravity (MG) models. In these models, deviations from the GR force law are often stronger for smaller astrophysical objects, making them ideal test beds of gravity in the local Universe. However, the actual behaviour of the modified gravity force also depends on the environment of the objects, and to make accurate predictions one has to take the latter into account. This can be done approximately using luminous objects in the local Universe as tracers of the underlying dark matter field. Here, we propose a new approach that takes advantage of state-of-the-art Bayesian reconstruction of the mass distribution in the Universe, which allows us to solve the modified gravity equations and predict the screening effect more accurately. This is the first of a series of works, in which we present our methodology and some qualitative results of screening for a specific MG model, $f(R)$ gravity. Applications to test models using observations and extensions to other classes of models will be studied in future works.

Equations

$${\Phi }_{\mathrm{env}}=\sum \frac{G{M}_i}{r_i}.$$

$${\mathcal{S}}_{\mathrm{EH}}=\int {\mathrm{d}}^{4}x\sqrt{-g}\frac{1}{16\pi G}\left[R+f(R)\right],$$

$$G_{\mu \nu }+f_R{R}_{\mu \nu }-g_{\mu \nu }\left[\frac{1}{2}f(R)-\square f_R\right]-{\nabla }_{\mu }{\nabla }_{\nu }{f}_R=8\pi G{T}_{\mu \nu }^{\mathrm{m}},$$

$$f_R\equiv \frac{\mathrm{d}f(R)}{\mathrm{d}R},$$

$$\square f_R=\frac{1}{3}\left[R-f_{R}R+2f(R)+8\pi G{\rho }_{\mathrm{m}}\right],$$

$$\mathrm{d}{s}^2=(1+2\Psi )\mathrm{d}{t}^2-a^2(t)(1-2\Phi ){\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j,$$

$${\nabla }^2\Phi$$

$${\nabla }^2{f}_R$$

$${\nabla }^2\Phi =4\pi G{a}^2\delta {\rho }_{\mathrm{m}}-\frac{1}{2}{\nabla }^2{f}_R.$$

$$G_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }f(R)\approx 8\pi G{T}_{\mu \nu }.$$

$${\nabla }^2{f}_R+\frac{\partial V_{\mathrm{eff}}\left(f_R\right)}{\partial f_R}=0,$$

$$\frac{\partial V_{\mathrm{eff}}\left(f_R\right)}{\partial f_R}=\frac{1}{3}{a}^2[\delta R\left(f_R\right)+8\pi G\delta {\rho }_{\mathrm{m}}].$$

$${\nabla }^2\Phi \approx \frac{16}{3}\pi G\delta {\rho }_{\mathrm{m}}=\frac{4}{3}{\nabla }^2{\Phi }_{\mathrm{GR}},$$

$$f(R)=-m^2\frac{c_1}{c_2}\frac{(-R/m^2{)}^n}{(-R/m^2{)}^n+1},$$

$$f_R=-\frac{c_1}{c_2^2}\frac{n(-R/m^2{)}^{n-1}}{[(-R/m^2{)}^n+1{]}^2}.$$

$$f_R\approx -n\frac{c_1}{c_2^2}{\left(\frac{m^2}{-R}\right)}^{n+1}.$$

$$f(R)\approx -m^2\frac{c_1}{c_2}\approx -2\Lambda \Rightarrow \frac{c_1}{c_2}=6\frac{{\Omega }_{\Lambda }}{{\Omega }_{\mathrm{m}}}.$$

$$\frac{c_1}{c_2^2}=-\frac{1}{n}{\left[3\left(1+4\frac{{\Omega }_{\Lambda }}{{\Omega }_{\mathrm{m}}}\right)\right]}^{n+1}{f}_{R0}.$$

$$\rho (r)=\frac{{\rho }_0}{r/R_s{\left(1+r/R_s\right)}^2},$$

$$c_{200}\equiv \frac{R_{200}}{R_s}=\mathrm{2.3.}$$

Full text

††pubyear: 2018††pagerange: Screening maps of the local Universe I – Methodology–References

1 Introduction

In recent years, modified gravity (MG) theories (Clifton et al., 2012; Joyce et al., 2015; Koyama, 2016, 2018) have been an active field of research in theoretical, observational and computational cosmology. One of the primary motivations for studying such models is to find alternative models to explain the accelerated cosmic expansion (Riess et al., 1998; Perlmutter et al., 1999), that avoid the theoretical difficulties in the standard $\Lambda$-cold-dark-matter ($\Lambda$CDM) paradigm. Other motivations for MG theories include attempts to find a more complete theory of gravity than General Relativity (GR) and to develop new ways to test the accuracy of GR; the latter is of particular interest since cosmological observations have entered the precision era, and started to allow accurate tests of gravity on length and energy scales vastly different from where GR has been conventionally validated (e.g., Will, 2014).

Being a long-range force, gravity acts on all length scales from sub-atomic to cosmological. Therefore, a deviation from GR’s prescription can in principle be measured on all these scales. Hence, although many of the MG models are originally proposed to tackle a cosmological problem, they can be tested in a huge array of environments or regimes, from laboratory experiments (see Brax et al., 2018, for a recent review), to Solar system and astrophysical objects (see Sakstein, 2018, for a recent review), and to observations at cosmological distances (see Koyama, 2016; Heymans & Zhao, 2018; Cataneo & Rapetti, 2018; Cai, 2018, for some recent reviews).

The requirement that any new theory of gravity must preserve the success of GR on small length scales has important implications on both theories and observations. Theoretically, one is confined to viable MG models, i.e., those that behave sufficiently closely to GR in environments such as the Solar System. One way to achieve this is through a screening mechanism (e.g., Khoury, 2010), by which modifications to the GR force law are suppressed in places of deep gravitational potential or in regions characterised by large gradients and/or by large Laplacians of the potential (like in the Solar system). Observationally, this implies that viable MG models must pass local tests of gravity by design, and thus we may need to turn to astrophysical and cosmological probes for complementary and potentially more stringent tests. The latter has been possible because cosmology concerns typically environments with shallow gravitational potentials or small values of its derivatives, where order unity deviations from GR can occur. MG theories are characterised by a variety of screening mechanisms, which means that a given probe could have very different constraining power for different models. Therefore, it is sensible to explore a wide range of potential cosmological and astrophysical probes. For example, for the popular $f(R)$ gravity model, in which the deviation from GR is controlled by a model parameter $f_{R0}$ (see more details below), the strongest constraints on $f_{R0}$ from cosmology suggest $|f_{R0}|\lesssim 10^{-6}$ (e.g., He et al., 2018; Leo et al., 2019), while the astrophysical constraints are claimed to be stronger (e.g., Jain et al., 2013b; Sakstein et al., 2014).

Even if one is interested in astrophysical constraints, it is often not sufficient to focus only on individual astrophysical objects. This is because, as we have mentioned above, the deviation from GR in many MG models is dependent on not just the astrophysical objects themselves but also the properties of their environments. A dwarf galaxy, for example, can be unscreened (i.e., it experiences a modified gravitational force) if placed in a low-density environment for a specific $f(R)$ model, but the same galaxy may well be screened (i.e., the deviation from GR is efficiently suppressed) if moved to dense environments such as close to a large galaxy cluster. In other words, screening is a nonlinear phenomenon, and the behaviour of (modified) gravity on small scales can not be cleanly disentangled from its behaviour on much larger scales. As a result, the precise knowledge of the total matter distribution in a large region (the environment) is necessary to accurately predict how a modified gravity model would affect the observational properties of an astrophysical object. Not knowing the former could introduce a uncontrolled systematical uncertainty to astrophysical tests of gravity.

Fortunately, observations of the local Universe have now become good enough for us to ‘reconstruct’ the relevant environmental properties needed to understand the screening. The first attempt of making use of such vital information was by Cabre et al. (2012), who estimated at the position of each observed galaxy the Newtonian potential, $\Phi_{\rm env}$, – which determines the screening efficiency for $f(R)$ gravity – from all other neighbouring galaxies:

[TABLE]

A similar but more sophisticated approach was taken by Desmond et al. (2018b), who considered also $\nabla\Phi_{\rm env}$ and $\nabla^{2}\Phi_{\rm env}$, which are quantities controlling the efficiency of other screening mechanisms than the chameleon mechanism exploited by $f(R)$ gravity. More effort was also devoted to obtaining the underlying mass distribution. In Cabre et al. (2012) only the galaxies detected by the Sloan Digital Sky Survey (sdss) were utilised to reconstruct $\Phi_{\rm env}$, while Desmond et al. (2018b) also included the contributions from (i) invisible dark matter haloes – haloes which do not host a galaxy – by using a simulation calibration, and (ii) the underlying total matter field (not necessarily in resolved haloes) at $z=0$, as obtained by a Bayesian density reconstruction technique (Lavaux & Jasche, 2016). The results of the works are 3D maps of the local Universe, which contain values of $\Phi_{\rm env}$, $\nabla\Phi_{\rm env}$ and $\nabla^{2}\Phi_{\rm env}$: these are called screening maps as these quantities determine the screening properties of the leading MG models as a function of location.

In this paper we introduce a new approach to obtain screening maps. Our approach also makes use of the reconstructed total matter field from the observed galaxy catalogues in the local Universe. However, instead of using this density field to calculate quantities such as $\Phi$ and $\nabla\Phi$, we directly use that to solve for the dynamical fields which are responsible for the modification of gravity (and for screening). The main motivation is that, while the above quantities qualitatively determine the efficiency of screening, the quantitative calculation is much more involving: as an example, in Vainshtein-type models it is not $\nabla^{2}\Phi$, but $\nabla^{2}\phi$ and $\nabla^{i}\nabla^{j}\phi\nabla_{i}\nabla_{j}\phi$, where $\phi$ is a scalar field propagating the modified gravity force, that determines the screening, and this is further complicated by the complex cosmic web. This approach, dubbed local universe screening Test Suite, or locusts, solves $\phi$ using the reconstructed density field by employing routines of the MG numerical simulation code ecosmog (Li et al., 2012, 2013a; Li et al., 2013b). This therefore requires the MG model to be clearly specified, and the study will be on a model-by-model basis. On the other hand, because there is only one observed local Universe, the underlying model of gravity – whichever it is – must reproduce the observationally-inferred matter density field. In particular, simulations of different gravity models should produce this same matter density field at $z\sim 0$, perhaps starting from different initial conditions. As a result, we only need to run one single $\Lambda$CDM simulation and output the matter field at various snapshots, and then the modified gravity routine in ecosmog can be used to calculate the screening properties of the model in these snapshots. This is much faster than full MG cosmological simulations, so that we can easily repeat the calculation for hundreds or even thousands (therefore the name locusts) of MG models that densely sample the model and parameter space. Another possibility enabled by this approach is the study of the time evolution of the screening map, which can be obtained by running the MG solver in ecosmog on several close output snapshots and then doing a finite difference.

In this paper we describe the methodology of the locusts simulations, and show the screening maps and some other physical quantities to demonstrate how it works. We do these using a specific MG model – chameleon $f(R)$ gravity – as an example, leaving the application of the method in astrophysical tests and extensions of it to include larger coverage of the local Universe and of more MG models into future works.

The layout of this paper is as follows. Section 2 briefly reviews the chameleon $f(R)$ gravity theory and the simulations used in this work. Section 3 presents our results, including visualisations of the simulated haloes and scalar field compared with the observed distributions of galaxies and galaxy groups, some simple statistics of the behaviour of the fifth force, and detailed properties of the COMA cluster. Finally, we conclude with a short summary and discussion in Section 4.

2 Methods and simulations

2.1 Constrained simulations of the local Universe

We make use of a constrained-realisation N-body simulation (labelled as CS) performed as part of the elucid project (Wang et al., 2014, 2016). The goal of the project is to reproduce the evolution history of our Local Universe by using the reconstructed initial density field from the observed galaxy catalogue. Here, we briefly summarise the reconstruction method as follows. First, a halo-based group catalogue is constructed from the SDSS DR7 galaxy catalogue with their positions and velocities having been corrected to “real space”. Then, a present-day density field is built according to the obtained halo catalogue. Finally, using the Hamiltonian Markov Chain Monte Carlos (HMC) algorithm with particle mesh dynamics, the initial condition is reconstructed from the present-day density field. For a more detailed description, we refer the reader to Wang et al. (2016). The method can effectively trace the $z=0$ massive haloes (\lower 3.22916pt\hbox{\sim}\hbox to0.0pt{\hss\raise 1.1625pt\hbox{>}}10^{13.5}{~{}M_{\odot}}) back to their initial condition, such that the reconstructed initial condition of our Local Universe can be used to study the evolution history of individual galaxy clusters and other cosmic web environments (see Fig. 3).

The initial condition reconstructed above, which is used in this work, features a periodic cubic box with a side length $500h^{-1}$Mpc and $1024^{3}$ dark matter particles. The mass of each simulation particle is $8.3\times 10^{9}h^{-1}{~{}M_{\odot}}$. The cosmological parameters are adopted from the best-fit WMAP5 cosmology (Dunkley et al., 2009): $\Omega_{\rm m}=0.258,\Omega_{\Lambda}=0.742,h=0.72,\sigma_{8}=0.8$ and $n_{\rm s}=0.96$.

2.2 The theoretical model

In this work we focus on a particular class of modified gravity models, $f(R)$ gravity, which is an extension to standard GR by replacing the Ricci scalar $R$ in the Einstein-Hilbert action of gravity with an algebraic function of $R$:

[TABLE]

where $G$ is Newton’s constant and $g$ is the determinant of the metric $g_{\mu\nu}$, with $\mu,\nu=0,1,2,3$.

The modified Einstein equation can be obtained by varying the action, Eq. (2), with respect to the metric $g_{\mu\nu}$, to obtain

[TABLE]

in which $G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ denotes the usual Einstein tensor, $\nabla_{\mu}$ is the covariant derivative compatible with $g_{\mu\nu}$, $\Box\equiv\nabla^{\mu}\nabla_{\mu}$ is the d’Alambertian, and $T^{\rm m}_{\mu\nu}$ is the energy-momentum tensor for matter. The quantity $f_{R}$ in this equation is an extra degree of freedom (a scalar field) of this model, defined by

[TABLE]

whose equation of motion can be obtained by taking the trace of Eqn. (3):

[TABLE]

where $\rho_{\rm m}$ is the density of non-relativistic matter. Therefore, the scalar field $f_{R}$ satisfies a second-order field equation of motion; this means that the modified Einstein equation, (3), which contains fourth-order derivatives of $g_{\mu\nu}$, can be rewritten as a standard second-order Einstein equation with a scalar field.

To investigate the evolution of cosmic structures in the Newtonian regime, we derive the perturbation equations in the Newtonian gauge on a flat Friedmann-Robertson-Walker (FRW) background:

[TABLE]

in which $\Phi=\Phi({\bf x},t)$ and $\Psi=\Phi({\bf x},t)$ are the gravitational potentials, which are functions of the physical time $t$ and the comoving coordinates ${\bf}=\{x^{i}\}$; $\delta_{ij}$ is the 3D spatial metric, and $a(t)$ is the scale factor, which is normalised to $a(t_{0})=a_{0}=1$ at the present day (a subscript 0 denotes the current value of a quantity throughout this paper, unless otherwise stated). In the quasi-static and weak field limits, the system of equations, (3) and (5), can be simplified respectively to:

[TABLE]

in which $\nabla^{2}$ denotes the 3D Laplacian operator, and the density and curvature perturbations are defined respectively as $\delta\rho_{\rm m}\equiv\rho_{\rm m}-\bar{\rho}_{\rm m}$ and $\delta R\equiv R(f_{R})-\bar{R}$; an overbar is used to denote the background value of a quantity. Eq. (7) can be recast in a new form:

[TABLE]

It can be seen clearly that the second term of the right-hand side of Eq. (9) represents a modification to the standard Poisson equation, and we can define $\Phi\equiv\Phi_{\rm GR}-\frac{1}{2}f_{R}$, where $\Phi_{\rm GR}$ is the Newtonian potential in GR, and $-\frac{1}{2}f_{R}$ can be identified as the potential of an additional force – the so-called fifth force, which is propagated by the scalar field $f_{R}$ – between matter particles. The fifth force is not detected in solar system or laboratory tests of gravity (Will, 2014), and these experimental tests place strong constraints on models like this.

To close Eqs. (7,8), one needs the relationship between $f_{R}$ and $R$ such that $\delta R$ can be expressed as a function of the scalar field as $f_{R}$: $\delta R\left(f_{R}\right)$. This can be done by specifying the functional form of $f(R)$, which satisfies the requirement that the resulting $f_{R}$ is a monotonic function of $R$. If $f(R)$ is a slowly-varying function of $R$, i.e., $|f_{R}|\ll 1$, the model has two desirable features:

• •

the terms involving $f_{R}$ in Eq. (3) can be neglected to a good approximation, reducing the Einstein equation to

[TABLE]

If one further approximates $f(R)\approx-2\Lambda$ (recall that $f(R)$ is taken to be nearly constant), with $\Lambda$ being the cosmological constant, then the background expansion history of this model can be made close to that of $\Lambda$CDM. In fact, with suitable choices of $f(R)$ the background expansion histories in the two models can be made exactly identical (He & Wang, 2014).

• •

if $|f_{R}|\ll 1$, one can have $\nabla^{2}f_{R}\sim 0$ and consequently from Eq. (9) we can see that the standard Poisson equation in GR is recovered. If this happens at least in high-density regions, it implies that the fifth force is suppressed in such regions, which can make the model compatible with current local tests of GR.

The suppression of the fifth force in the limit $|f_{R}|\ll 1$ is the result of a suitable choice of $f(R)$; it is a dynamical effect called the screening mechanism. $f(R)$ gravity is a representative example of a wider class of models, called the chameleon model (Khoury & Weltman, 2004), in which the suppression (or screening) of the fifth force works as following: the scalar field $f_{R}$, which propagates the fifth force between matter particles, satisfies Eq. (8), which can be rewritten as

[TABLE]

where $V_{\rm eff}\left(f_{R}\right)$ is an effective potential of the scalar field, given by

[TABLE]

The potential $V_{\rm eff}$ characterises the interactions of the scalar field with itself (the first term on the right side of Eq. (12)) and matter (the second term). For a choice of $f(R)$ such that $V_{\rm eff}\left(f_{R}\right)$ has a global minimum at $f_{R}=f_{R,\rm min}$ and $f_{R,{\rm min}}\rightarrow 0$ as $\delta\rho_{\rm m}\rightarrow\infty$, the fifth force can be suppressed in high-density regions as desired, therefore evading the stringent local constraints on it. Because the behaviour of the fifth force is dependent on the environmental density, the screening mechanism is called the chameleon mechanism. In regions where $|\delta\rho_{\rm m}|\ll 1$, on the other hand, the curvature perturbation $|\delta R|\ll 1$ and so from Eqs. (7,8) one can derive that

[TABLE]

which means that gravity is enhanced by a factor of $1/3$ – an effect that is potentially testable using cosmological observations.

The actual behaviour of the fifth force in $f(R)$ gravity is more complicated that the above intuitive picture, and an accurate solution has to be made by numerically solving Eq. (8) given a matter configuration. In this context, to solve for $f_{R}$ at a given position we need its solution in the neighbourhood as the boundary condition – in other words, to know for certain whether a given cosmological object, such as a star or galaxy, is screened, we need to solve Eq. (8) in a large region encompassing this object, and the solution in that region in turn depends on further nearby regions, and so on. In this picture, screening of the fifth force for an object can be achieved in two ways:

• •

self screening: if the objective is massive enough, it alone can make $|f_{R}|$ small inside and/or nearby, therefore screening the fifth force it feels;

• •

environmental screening: if the object is not massive enough to self screen, but lives near some much larger objects, then $|f_{R}|\ll 1$ can still be satisfied inside and/or near it, causing a suppression of the fifth force it experiences.

To use astrophysical objects in the local Universe to test the fifth force, then, we cannot reliably treat those objects as isolated bodies living on the cosmological background, but have to take into account their larger-scale environments. For this reason a constrained realisation simulation as described in Section 2.1, where the matter distribution mimics that in the real observed local Universe, is ideal as it offers a way to more realistically model the effect of environments in the chameleon screening.

In this work we shall set up the general strategy to carry out constrained realisation simulations in modified gravity models, and present some first results to show how it works. We leave detailed analyses of these simulations that lead to constraints on model parameters to future works. For concreteness, we use the $f(R)$ model proposed by Hu & Sawicki (HS; 2007) as example. This model is given by specifying

[TABLE]

where $m^{2}\equiv 8\pi G\bar{\rho}_{{\rm m}0}/3=H^{2}_{0}\Omega_{\rm m}$ is a parameter of mass dimension 2, $\Omega_{\rm m}$ the density parameter for non-relativistic matter, $H_{0}$ the present-day value of the Hubble expansion rate, and $n$, $c_{1}$ and $c_{2}$ are dimensionless model parameters. The scalar field, Eq. (4), takes the following form:

[TABLE]

To see whether this model can have a background expansion history close to that of standard $\Lambda$CDM, let us consider a $\Lambda$CDM model with $\Omega_{\rm m}\approx 0.3$ and $\Omega_{\Lambda}=1-\Omega_{\rm m}\approx 0.7$, for which we find $|\bar{R}|\approx 40m^{2}\gg m^{2}$, and therefore

[TABLE]

For $n\sim 1$ and $c_{1}/c_{2}^{2}\lesssim 1$, we then have $|f_{R}\left(\bar{R}\right)|\ll 1$, which is the condition by which the background expansion history is close to $\Lambda$CDM, and

[TABLE]

Therefore, once we have specified an (approximate) $\Lambda$CDM background history (by which $c_{1}/c_{2}$ is fixed), the HS $f(R)$ model then has two free parameters, $n$ and $c_{1}/c^{2}_{2}$. The latter is related to the present-day value of the background scalaron, $f_{R0}$,

[TABLE]

The choice of $f_{R0}$ and $n$ fully determines the model.

2.3 The locusts simulations

In this subsection we introduce the locusts simulation suite and briefly describe the simulation technique used.

The locusts simulations are a suite of simulations of various modified gravity models, all starting from an identical initial condition, which itself is obtained as described in Section 2.1. Therefore, they are the first attempt to realistically simulate our local Universe in the context of modified gravity. In particular, one of the primary objectives of locusts is to obtain screening maps, namely a map to show the screening properties at different spatial locations in the local Universe. As stated in the introduction, such screening maps can provide vital information for both cosmological and astrophysical tests of gravity.

While the basic idea is general, in this work we focus on the chameleon $f(R)$ gravity model described in Section 2.2 as explicit example. In particular, we shall specialise to the case of $n=1$, and run simulations for 20 different values of $|f_{R0}|$, ranging from $10^{-7}$ to $10^{-6}$. This parameter range is still compatible with the currently most stringent constraints on $f_{R0}$ from cosmological observations (see, e.g., Cataneo et al., 2015; Liu et al., 2016; Peirone et al., 2017).

The chameleon $f(R)$ simulations used in this work have been done using the ecosmog (Li et al., 2012) code, which is a modified version of the publicly available $N$-body and hydrodynamical simulation code ramses (Teyssier, 2002). This is a particle-mesh code employing the adaptive-mesh refinement technique to achieve high force resolution in dense regions, and parallelised using message passing interface. ecosmog extends ramses by solving the nonlinear field equations which arise from various modified gravity models numerically by the multigrid relaxation method. For details about the implementation in different classes of models, see Li et al. (2012, 2013a); Li et al. (2013b) and references therein. We use an optimised version of ecosmog for the Hu-Sawicki $f(R)$ model, as described in Bose et al. (2017), which is based on a more efficient algorithm to solve the $f(R)$ field equation.

Even with the algorithm optimisation from Bose et al. (2017), running a suite of $>\mathcal{O}(20)$ simulations with different $f(R)$ gravity parameters is still computationally expensive for the resolution and particle number used in locusts. Fortunately, as explained in the introduction, the idea behind locusts does not require us to run full simulations of modified gravity, but only needs one simulation to $z=0$ which provides a mock universe with a underlying matter density field. This underlying density field must be as close to the observationally-inferred density field in the local Universe as possible, and any gravity model should reproduce this same underlying density field. This might be achieved by tuning the initial conditions of the simulations in different models, but the details are not our concern here. Apparently, the simplest way to achieve this is to only run the full simulation (from $z_{\rm ini}=80$ to $z=0$) in the $\Lambda$CDM model, while for the $f(R)$ models we simply run the ecosmog for a few steps – respectively on the particle snapshots of the $\Lambda$CDM simulation at various redshifts – to calculate the behaviour of the scalar field and the fifth force at the redshifts of interest to us. Put in other words, the particle evolution in the locusts simulations is done using Newtonian gravity, while the evaluation of the screening map is done using the complete modified gravity solver.

As mentioned above, the evolution of particle positions and the calculation of the scalar field and screening properties are both performed using ecosmog, which is based on the ramses code. As a rough estimate of the level to which we can trust the simulation density field (i.e., the typical difference between different simulation codes at our resolution), in Fig. 1 we have compared the matter power spectra at $z=0$ predicted by ecosmog and the gadget-2 code (Springel, 2005). We can see there is good agreement – within $1\%$ for $k<3~{}h$Mpc*-1* and $4\%$ for $k<6~{}h$Mpc*-1*. The difference at small scales ($\lesssim\mathcal{O}(1)h$Mpc*-1*) is expected to be much smaller than the typical uncertainty in the density reconstruction.

As another sanity check, in Fig. 2 we plot the halo mass function from the $\Lambda$CDM simulation at $z=0$ (squares) compared with the Tinker et al. (2008) fitting formula (solid line). The halo catalogues in this and other figures of this paper are identified using the phase-space friends-of-friends halo finder rockstar (Behroozi et al., 2013), and the halo mass $M_{200}$ denotes the mass within $R_{200}$, the radius within which the average density is 200 times that of the critical density of the Universe at the halo redshift, $\rho_{\rm crit}(z)$. The simulation output agrees well with Tinker et al. (2008) apart from the high-mass end, and at $M_{200}\lesssim 10^{12}h^{-1}{~{}M_{\odot}}$ (which correspond to haloes with $\lesssim 100$ particles, for which the mass function becomes incomplete due to the low resolution).

3 Results

This section contains the main results of this work. We start with some visualisation and general properties of the fifth force throughout the simulation box, then move on to study statistical properties of the screening maps and the screening around prominent structures in the local Universe, such as the Coma cluster and the SDSS Great Wall.

3.1 Visualisation

Fig. 3 is the visual comparison of a slice taken from the SDSS group catalogue (left panel, in which groups are shown as black dots) with an extraction of the simulation box that is supposed to represent the same region (middle and right panes); the middle panel shows the dark matter density field in the region, while the right panel shows the corresponding scalar field configuration for $f_{R}$. Both simulation results are at $z=0$, and for the right panel a particular $f(R)$ model with $|f_{R0}|=10^{-6}$ is shown for illustration purpose.

We see from Fig 3 that the constrained simulation has successfully reproduced the large-scale structures observed from the SDSS catalogue, noticeably the filamentary patterns on scales of tens of Megaparsecs and above. In particular, from the dark matter distribution we can see clearly the SDSS Great Wall found at $X=-230~{}h^{-1}\rm Mpc{}$ and extending in the vertical coordinate from $-100$ to $50~{}h^{-1}\rm Mpc$ (Gott et al., 2005). The scalar field $f_{R}$, as shown in the right panel, behaves as expected from the chameleon screening mechanism: its value is closer to [math] near clusters and filaments, while approaching the background value $f_{R0}$ further away from these structures. In particular, we note that deep inside void regions the scalar field is nearly uniform, suggesting that the fifth force, which is the gradient of the scalar field, is weak there222However, the fifth force can still be significant near (often small) matter clumps inside these voids, and we shall return to this point later below..

To better compare our constrained simulation with the observational data, we zoom-in on a small region centred on the Great Wall. The results are shown in Fig. 4, where the upper left panel is an enlarged view of the matter density field from the central panel of Fig. 3 using the same colour bar.

Dark matter haloes identified in the constrained simulation are shown as black open circles in the lower left panel of Fig. 4, where the radius of each circle is proportional to the mass of the halo it represents. Overplotted on top are the SDSS galaxy groups which are shown as red filled circles. We find a very good agreement between the positions of simulated haloes and those of SDSS groups, although some outliers do exist. This comparison represents a beautiful illustration of how well the constrained simulation reproduces the large-scale distribution of galaxies.

In the two panels on the right-hand side of Fig. 4 we show the SDSS groups and galaxies overplotted on screening maps for the same zoomed-in region. The coloured map in the upper right panel shows again the scalar field, where we can see more clearly that the scalar field closely traces matter distribution and is nearly homogeneous in low-density regions. The blue and red dots represent SDSS blue and red galaxies respectively in this panel, and the latter also trace well the simulation matter distribution. This suggests that we can use the simulated screening map to predict the scalar field value and fifth force ratio at the positions of the observed objects. In the lower-right panel we show this for groups (filled circles whose sizes indicate the masses of the groups they represent) and red galaxies (dots) – here the colour is used to illustrate the fifth force ratio at the positions of the groups and galaxies, and we can see that the objects are more screened in dense regions than in underdense regions.

3.2 Generic behaviours of the fifth force

Before quantifying the fifth force effects, let us present some results of the general behaviour of the fifth-force-to-standard-gravity ratio across the whole simulation volume. In Fig. 5, we have shown this force ratio at the positions of $10^{5}$ particles randomly selected from the simulation box, where each dot represents the measured value at a simulation particle. The different panels are for different $|f_{R0}|$ values, starting from the least screened case with $10^{-6}$ at the upper left and ending at the most strongly screened case with $10^{-7}$ at the lower right. The colour indicates the frequency that particles appear with given standard gravity (horizontal axis) and fifth force (vertical axis) values. Comparing amongst the different panels and comparing simulation results with analytical linear perturbation prediction (red solid line), we observe the following features:

• •

In high-density regions, where the magnitude of the standard gravity force is large, the fifth force is generally strongly screened, and the points are well below the red line, which represents the case that the fifth force has $1/3$ of the strength of standard gravity.

• •

In the regime of intermediate magnitudes of standard gravity, representative for smaller haloes and filaments, the fifth force ratio agrees with linear theory prediction well for the weakly screened models. However, as $|f_{R0}|$ decreases, stronger screening shows up even in this regime; for example, in the last row, we can see clearly that the red dots are well below the red solid line.

• •

In the regime of weak standard gravity, i.e, the left end of each panel, which is representative of void regions, the fifth force ratio falls below the red solid line again. This is because the fifth force, unlike standard Newtonian gravity, is a short-ranged force that decays exponentially beyond the Compton wavelength of the scalar field. This implies that the standard gravity exerted by particles outside the void regions can reach the inner part of these voids, while the fifth force cannot, leading to a suppressed force ratio between the latter and the former. This can also be understood through the observation that in void regions, e.g., Fig. 4, the scalar field is nearly homogeneous and so the fifth force becomes weak. Note that the deviation from linear (anaytic) prediction of the fifth-force-to-standard-gravity ratio does not happen in Vainshtein screening models, e.g., Fig. 3 of Falck et al. (2014) since there the fifth force is long range.

Figure 6 is similar to Figure 5, but instead of the force ratio, it shows the cosine of the angle $\theta$ between the fifth and standard gravity forces. If linear theory works perfectly, the fifth force should be $1/3$ of the strength of standard gravity and the directions of the two forces would be the same. While this is the case for most particles from the intermediate gravity regime (the red and orange regions), we can see that for the strong and weak gravity regimes this is not true. The reason is the same as for the behaviour of the force ratio shown in Fig. 5, namely in regions of deep Newtonian potential (and therefore strong standard gravity) the fifth force is suppressed by the chameleon screening mechanism, while inside voids the fifth force from matter in surrounding regions suffers from the Yukawa exponential decay – because of that, and given the irregular matter distribution, the fifth force and Newtonian gravity on a particle in a weak-gravity region can receive contributions from different neighbouring particles as vector additions, and so they do not necessarily have the same direction; this is particular true because the fifth force has much smaller amplitude (note that even in the weak-gravity regions there are still particles, so the fifth force is not exactly zero).

Finally, we are interested to check how the force ratio depends on large-scale environment. For this we used the nexus+ method (Cautun et al., 2013) to identify the various cosmic web environments: nodes, filaments, sheets and voids. The nodes correspond to the densest regions, filaments to 1D linear structures, sheets to 2D wall-like planar densities and voids to underdense regions. These morphological environments have been found by first using the Delaunay Tessellation Field Estimator (Schaap & van de Weygaert, 2000; Cautun & van de Weygaert, 2011) to calculate the density field on a regular grid with a $1~{}h^{-1}\rm Mpc$ grid spacing. Then, nexus+ calculates the eigenvalues, $\lambda_{i}$ with $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$, of the Hessian matrix of the smoothed density field, which are used to classify the web environments. The exact procedure is based on some rather complex functions of the Hessian eigenvalues, however the result can be qualitatively understood as: nodes correspond to regions with $\lambda_{1}\approx\lambda_{2}\approx\lambda_{3}<0$, filaments to regions with $\lambda_{1}\approx\lambda_{2}<0$ and $\lambda_{2}\ll\lambda_{3}$, sheets to $\lambda_{1}<0$ and $\lambda_{1}\ll\lambda_{2}$, and voids to everything else. For a detailed comparison of the nexus+ technique to other web finders, please see Libeskind et al. (2018).

The resulting cosmic web is dominated in terms of volume by voids, which occupy ${\sim}80\%$ of the volume but contain only ${\sim}15\%$ of the total mass budget. In terms of mass, the filaments are the most important environment, containing over half of the mass budget but filling only $6\%$ of the cosmic volume (Cautun et al., 2014). Most of the massive haloes, with $M_{200}\gtrsim 5\times 10^{13}h^{-1}{~{}M_{\odot}}$, are found in nodes, while filaments contain the majority of lower mass haloes with mass $M_{200}\gtrsim 10^{11}h^{-1}{~{}M_{\odot}}$ (Ganeshaiah Veena et al., 2018). In contrast, sheets and especially voids correspond to below average densities and are mostly devoid of haloes with masses above $10^{12}h^{-1}{~{}M_{\odot}}$. This means that the majority of bright galaxies, that is with stellar masses above $10^{9}h^{-1}{~{}M_{\odot}}$, are found in either the filaments or nodes of the cosmic web (Ganeshaiah Veena et al., 2019).

Figure 7 shows the same fifth-force-to-standard-gravity ratio as in Figure 5, but for particles found in voids (upper left), sheets (upper right), filaments (lower left) and nodes (lower right). To increase the clarity of the plots, we have only shown the results for $|f_{R0}|=10^{-6}$ and neglected pixels which represent particles that are smaller than $0.2$ thousandth of the total particle number. The overall behaviour is similar to what Figure 5 shows, but there is also a clear distinction between the various web environments. For example, the long drop-off tail with small force ratio but strong standard gravity forces seen in Figure 5 is mainly due to particles from nodes (high-density environments), while the drop-off from the analytical line at weak standard gravity forces is dominated by low-density environments such as voids and sheets, as explained above.

3.3 The Coma Cluster

The constrained initial condition used in our simulations has a limited volume, with objects such as the Local Group and Virgo Cluster not included. Therefore, here we select the object corresponding to the Coma cluster in our simulation volume, to illustrate the behaviour of the modified gravity force in massive objects.

Coma is a cluster at a distance of about $100$ Mpc from us, with over 1000 member galaxies and a total mass of $\sim 10^{15}M_{\odot}$. The dark matter halo we identify from our simulation as the counterpart of Coma333In what follows we shall refer to both the real Coma cluster and the counterpart dark matter halo from our simulation as ‘Coma’; the context should make clear what we mean. is found to have a mass of $M_{\rm 200}=7.7\times 10^{14}h^{-1}{~{}M_{\odot}}$ and halo radius $R_{\rm 200}=1.5h^{-1}$Mpc. As a first visual inspection, in Figure 8 we show the projected density in a $40\times 40$ $\left(h^{-1}{\rm Mpc}\right)^{2}$ field of view centred around the Coma halo, with a projection depth of $5h^{-1}$Mpc. On top of this, the observed Coma member galaxy groups are also shown as black open circles. We can see that the galaxy groups broadly follow the same clustering pattern of high-density regions in the projected map.

In Figure 9 we show the fifth-force-to-standard-gravity ratio in the same region as Figure 8, for four different $f_{R0}$ parameter values as indicated in the legends of the four panels. As expected, in the inner regions of the cluster screening is more efficient, due to the deeper Newtonian potential there. As $|f_{R0}|$ decreases, screening becomes more efficient; for $|f_{R0}|=10^{-6}$, which is the model with the weakest screening, the fifth force is strongly suppressed (with force ratio $F_{5\rm th}/F_{\rm standard}\lesssim 0.01$) only up to $\sim 2~{}h^{-1}$Mpc from the cluster; as $|f_{R0}|$ decreases, this screened region (blue or red in colour) expands outwards, with the nearby filamentary structures and some smaller haloes scattered around now also featuring a strongly suppressed fifth force.

Finally, Figure 10 shows the density (top panel) and force ratio (right) profiles in the Coma halo. The density profile is obtained by computing the spherically averaged densities within logarithmic radial bins from the halo centre found by rockstar, and we show the result out to $5h^{-1}$Mpc from the halo centre, with the halo radius $R_{200}$ indicated by the dashed vertical line. The profile can be well fitted by the Navarro-Frenk-White (Navarro et al., 1996, NFW) formula,

[TABLE]

in which $\rho_{0}$ is a characteristic density and $R_{s}$ the scale radius, and the best-fit value of $R_{s}$ is found to be $0.65~{}h^{-1}\rm Mpc$, so that the halo concentration is

[TABLE]

The best-fit NFW profile for this halo is plotted as the black dotted line. Here, we use the particles within $R_{200}$ to fit the NFW profile, and the concentration would be greater if we extend the fitting to larger radii.

The lower panel of Fig. 10 shows the force ratio profiles in the same halo for the different $|f_{R0}|$ values, decreasing from top to bottom. This is obtained similarly as the density profiles, but the spherical average is now over the force ratio at the positions of all simulation particles for each radial bin. As is typical for haloes of this mass, the fifth force is efficiently suppressed inside $R_{200}$ even for the model which deviates most from GR ($|f_{R0}|=10^{-6}$). Another interesting feature is that the shapes of the force ratio profiles are similar for all $f_{R0}$ values, and the only difference is in the amplitudes. This is a natural consequence of using the same density profile for all our fifth force calculation in all models. Note that baryonic processes associated with the cluster and galaxy evolution may also play a role in the redistribution of matter and in the suppression of the fifth force. It would be beneficial to test in a zoom-in hydrodynamics simulations, which will be left as a potential project for the future.

4 Summary and conclusions

We have developed a new methodology for testing modified gravity theories using astrophysical probes (Jain et al., 2013a; Sakstein, 2018) based on constrained simulations of the local Universe. This method takes advantage of the recent developments in reconstructing the density field (and its initial conditions) of the local Universe (e.g., Wang et al., 2014; Lavaux & Jasche, 2016; Sorce et al., 2016; Carlesi et al., 2016; Wang et al., 2016), which provide a way to realistically include the environmental effects that are often important in quantitatively determining the behaviour of gravity in MG models. Our method combines the large-scale density field from these reconstruction schemes with the fully nonlinear numerical solution to the MG equations achieved through the ecosmog code. Assuming that the matter field at the low-$z$ Universe444We have argued that this is a good approximation for most cases, but note that this approximation is not needed: full simulations with MG are possible though more time consuming. behaves similarly in realistic MG models and $\Lambda$CDM, this will make it possible to create screening maps for a large number of MG models and parameter choices at a relatively low cost.

This is the first of a series of papers, where we have presented the methodology and, as a proof of concept, shown screening maps and some statistical properties that one can extract. As demonstrated in Figures 3 and 4, the simulated halo distributions and the resulting screening maps show good visual agreements with the distribution of sdss galaxies and groups, indicating that the method is capable of telling, for a given MG model, which parts of the local Universe and how well they are screened. The force behaviours displayed in Figures 5, 6 and 7 also agree with expectations based on the properties of chameleon screening, with smaller $|f_{R0}|$ values generally corresponding to more strongly suppressed fifth forces. In particular, Figs. 5 and 7 show that the fifth force is suppressed in not only high-density regions (where the Newtonian force is strong) but also in low-density regions. This seemingly counter-intuitive effect is due to the fact that the Compton wavelength of the scalar field in the models is small, such that there are few particles the fifth forces produced by which could propagate into deep voids (Paillas et al., 2019).

As a specific example, we have analysed in greater detail the dark matter halo from our simulation box which is the counterpart to the Coma cluster. Figure 9 shows that for all models considered here, the central region within $R\sim 2~{}h^{-1}$Mpc is well screened and so gravity there should behave like GR. On the other hand, in the stronger screening cases, where $|f_{R0}|\rightarrow 10^{-7}$, the screened region becomes larger, showing that the presence of a massive body can screen its smaller neighbours. This can be seen more clearly in the lower panel of Figure 10, which shows that within the virial radius the fifth force has never exceeded $\sim 0.01\%$ of the Newtonian force for all models considered.

Screening maps as shown in this paper can be invaluable for astrophysical tests (e.g., Cabre et al., 2012; Desmond et al., 2018a; Desmond et al., 2018b, c; Desmond et al., 2019), and they will enable these tests to become more reliable. However, the application of these maps in real tests are beyond the scope of this paper and will be left as future work. Also, one slight limitation of the current maps is that the Local Group is not included in the sdss field, but this is not a practical restriction for our method considering that constrained realisations that include the Local Group have now been produced by various groups. One interesting possibility is to use such constrained initial conditions to run very-high-resolution zoom-in simulations, possibly with baryons, which realistically reproduce the basic observational properties of the Milky Way Galaxy, and use that to quantify the screening inside the Milky Way and in the Solar system.

Another point that merits further investigation is related to the uncertainty in the estimation of galaxy host halo mass. The reconstruction method used here, from Wang et al. (2016), can effectively trace the $z=0$ massive haloes (\lower 3.22916pt\hbox{\sim}\hbox to0.0pt{\hss\raise 1.1625pt\hbox{>}}10^{13.5}{~{}M_{\odot}}) back to their initial condition, but its accuracy in recovering the matter distribution is poorer on smaller scales: the uncertainty was found to be respectively 0.23 dex on $2~{}h^{-1}\rm Mpc{}$ and 0.1 dex on $4~{}h^{-1}\rm Mpc{}$ scales. Therefore, in the screening map, while the environmental screening effect caused by large-scale structure can be reliably modelled, the uncertainty related to the galaxy host halo mass estimation can be another source of error, in particular for models with small Compton wavelengths of the scalar field. The impact of this uncertainty on the screening map and consequently on the error of model parameter constraints, however, can be assessed by creating a large number of screening maps using modified density fields in regions that correspond to dark matter haloes (to reflect the mass error of those haloes) before calculating the scalar field, and using these maps to quantify the scatter.

Acknowledgements

We thank the anonymous referee for detailed comments that have helped us improve the paper. SS and BL are supported by the European Research Council (ERC) via grant ERC-StG-716532-PUNCA. MC and BL are supported by STFC Consolidated Grants ST/P000541/1, ST/L00075X/1. MC is also supported by the EU Horizon 2020 research and innovation programme under a Marie Skłodowska-Curie grant agreement 794474 (DancingGalaxies) and by the ERC grant DMIDAS [GA 786910]. HYW is supported by the National Key R&D Program of China (grant No. 2018YFA0404503) and NSFC 11421303. JW is supported by STRD grants (2015CB857005, 2017YFB0203300) and NSFC grant 11873051. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1, ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.