# Cosmological gravity on all scales III: non-linear matter power spectrum in
phenomenological modified gravity

Sankarshana Srinivasan  Daniel B Thomas  and Richard Battye

(April 2020)

###### Abstract

Model-independent tests of gravity with cosmology are important when testing
extensions to the standard cosmological model. To maximise the impact of these
tests one requires predictions for the matter power spectrum on non-linear
scales. In this work we validate the ReACT approach to the non-linear matter
power spectrum against a suite of phenomenological modified gravity N-body
simulations with a time-varying gravitational constant, covering a wider range
of parameter space than previously examined. This vanilla application of ReACT
has limited range and precision due to the different concentration-mass
relation $c(M)$ that occurs when gravity is modified. We extend this approach
with a fitting function for a modified concentration-mass relation, allowing
for accurate (1$\%$) computation of the matter power spectrum up $k=2\,h\,{\rm
Mpc}^{-1}$ across a substantial range of parameter space. This fitting
function allows precision model-independent tests of modified gravity to be
carried out using the data from upcoming large scale structure surveys.

##  1 Introduction

Over the last two decades, the field of cosmology has entered an era of
precision measurements. The current cosmological paradigm, the $\Lambda$CDM
model contains Cold Dark Matter (CDM) and dark energy (sourced by the
cosmological constant $\Lambda$) as its main components, but the underlying
nature of both remains unknown. The validity of General Relativity (GR) is a
crucial assumption in this picture since both CDM and $\Lambda$ are inferred
from solely gravitational observations. In addition, efforts to detect dark
matter directly or indirectly have not yet been successful, while quantum
field theory predicts a $\Lambda$ value that is many orders of magnitude
larger than that inferred from cosmological measurements. These issues form
part of the motivation for modifying the law of gravity and the establishment
of a vast space of modified gravity models [1, 2] in order to explain the
physics attributed to the dark sector. Efficiently testing this model space is
a major challenge in cosmology.

The assumption that GR is the correct theory of gravity involves an
extrapolation of several orders of magnitude from small scales (compact object
mergers and solar-system tests) to the cosmological scales probed by
measurements of the Cosmic Microwave Background (CMB). Upcoming experiments
such as the Euclid satellite [3] 111https://www.euclid-ec.org, the Nancy Roman
telescope (NRT) [4] 222https://roman.gsfc.nasa.gov, the Vera Rubin Observatory
(VRO) [5] 333https://www.lsst.org and Square Kilometre Array (SKA)
444https://www.skatelescope.org[6] will generate an unprecedented volume of
data in the so-called non-linear regime, where the density contrast $\delta\gg
1$, i.e., the fluctuation away from the average density of the Universe is
large where GR has not been adequately tested. In the context of modified
gravity, one requires a parameterisation of modified gravity models [7, 8]
that encompasses a large region of the parameter space. The lack of a robust
tool that takes advantage of the additional constraining power from including
non-linear scales to predict cosmological observables across the range of
scales for general models of modified gravity results in strict data-cuts that
excise non-linear scales from the data analysis, severely limiting
constraining power (see, for example [9]).

Previous studies of modified gravity in the non-linear regime have primarily
been $N$-body simulations of specific models such as $f(R)$ or Dvali-
Gabadadze-Porrati (DGP) gravity (see e.g. [10, 11]). Such studies are useful
to understand structure formation as a function of model-parameters, and
explore observational signatures in specific theories [12]. However, it is not
practical to run dedicated $N$-body simulations for every modified gravity
model in order to predict cosmological observables, as the associated
computational cost would be too large. In addition, more model-independent
approaches allow for a more comprehensive null test of $\Lambda$CDM and don’t
rely on us having already come up with the correct alternative.

One approach for testing gravity in a model-independent way on all scales is
described in [13]. In this approach, the standard parameterisation in linear
theory involving two parameters is extended to all cosmological scales: $\mu$
is a change in the strength of gravity, and $\eta$ is the “slip”, a possible
difference between the two Newtonian gauge scalar metric potentials. Both of
these can in principle be almost arbitrary functions of time and space. In our
simulation work we follow “the maximally-phenomenological pixels” approach
described in [13] in which the modified gravity parameters consist of
independent piece-wise constant functions composed of independent bins in time
555The extension for $\mu$ and $\eta$ to be functions of space as well is not
conceptually difficult, however we focus on time dependence as a first step
and to ensure we understand this simpler case well before generalising
further. We leave the generalisation of our N-body simulations to scale
dependence to future work. as this is most suited to our philosophy of model-
independence.

In [14] we studied the phenomenology of $N$-body simulations with binned time-
dependent modified gravity parameters. We improved on the previous attempt to
study such models [15, 16] in which $\mu$ was set to be a constant value at
$z<50$, i.e., constant throughout the simulation redshift range. We focused on
the matter power spectrum $P(k)$ measured from the simulations and the weak-
lensing observables computed from $P(k)$. We found that the ReACT formalism
[17, 18, 19] is most successful among those available in the literature in
reproducing the cosmological observables computed in or from our simulations.
We designed a pipeline to compute cosmological observables in the non-linear
regime using ReACT. Here, we validate this pipeline in detail across a much
wider region of the (phenomenological) modified gravity parameter-space. In
particular, we quantify the range of scales any cosmological observable can
reliably be computed using ReACT. In addition, there were indications in [14]
that this implementation would not be precise enough over a wide enough range
of scales to maximise the gain from upcoming surveys, and also that this issue
could be alleviated by modifying the concentration-mass relation. We
investigate this issue in detail and provide an extension to the fit of [14]
that substantially improves the range of precision of the fitting function.
Our goal is to validate ReACT in order to run modified gravity forecasts for
binned $\mu(z)$ and $\eta(z)$.

The structure of this paper is as follows. In section 2 we review our
parameterisation, define the terms we will use to quantify the results from
our simulations and present our adapted implementation of ReACT with binned
modified gravity functions. We define our pipeline for computing $P(k)$ and
the weak-lensing convergence spectrum using ReACT. We also discuss how one can
improve the accuracy of ReACT by modifying halo concentrations. We then
present the suite of $N$-body simulations that we have run to validate our
ReACT pipeline in section 3. In section 4, we quantify the accuracy and
precision of ReACT in predicting $P(k)$ and weak-lensing observables against
our simulations both with and without the concentration modification. We
discuss the relevance of these results in the literature and the impact they
may have for cosmological parameter forecasts in the context of modified
gravity for surveys such as Euclid, NRT, VRO and SKA in section and conclude
in section 5.

##  2 Modified gravity framework

As described in [14], our parameterisation may be formulated as

$\displaystyle\frac{1}{c^{2}}k^{2}\tilde{\phi}_{\rm P}$ $\displaystyle=$
$\displaystyle-\frac{1}{c^{2}}4\pi a^{2}\bar{\rho}G_{\rm
N}\mu(a,k)\tilde{\Delta}\,,$ (2.1) $\displaystyle\tilde{\psi}_{\rm P}$
$\displaystyle=$ $\displaystyle\eta(a,k)\tilde{\phi}_{\rm P}\,,$ (2.2)

where $\bar{\rho}$ is the background density,
$\tilde{\Delta}=\tilde{\delta}-\frac{\dot{a}}{a}\frac{3}{c^{2}k^{2}}ik_{i}\tilde{v}_{i}$
is the gauge-invariant density contrast in Fourier space, $G_{\rm N}$ is
Newton’s constant, $\phi_{P}$ and $\psi_{P}$ are the standard Newtonian
gravitational potentials666$g_{00}=-1-\frac{2\phi_{P}}{c^{2}}$;
$g_{ij}=a^{2}\delta_{ij}\left(1-\frac{2\psi_{P}}{c^{2}}\right)$. (normally
found to be equal in GR) and $\mu(a,k)$ is a dimensionless function of space
and time representing a change to the strength of gravity. These equations can
be derived from the FLRW metric expanded up the first Post-Friedmann order
(which is one step beyond leading order) at which the equations describe
structure formation on all scales (see [13] for details). An important
consequence of this parameterisation is that the dynamics of massive particles
are purely governed by (2.1), i.e., one can build and run dark matter only
$N$-body simulations where the structure formation is completely specified by
a single modified gravity parameter, $\mu$. The second parameter $\eta$, the
so-called slip parameter purely affects photon geodesics and, therefore, can
be modelled in post-processing.

###  2.1 Applying the ReACT framework to modified gravity with binned $\mu$

We showed in [14] that the ReACT framework performed significantly better than
other candidates from the literature including the standard halofit
prescription and the standard halo model in predicting $P(k)$ measured from
our simulations. We showed how we can take the publicly available
implementation of the ReACT code and apply it to our binned models. We call
this implementation the ‘vanilla’ implementation throughout this paper.

The fundamental idea at the centre of ReACT is that the accuracy of the halo
model can be significantly improved by measuring the so-called ‘non-linear
response of the power spectrum’ [20] i.e., the ratio of matter power spectra
between two cosmologies that share the same linear power spectrum. The halo
model reaction formalism computes the non-linear matter power spectrum for
beyond $\Lambda$CDM cosmologies from the so-called ‘pseudo spectrum’, which is
the $\Lambda$CDM power spectrum with the the initial conditions modified to
match the growth in the modified gravity model. The use of the pseudo spectrum
ensures that the transition from the 2-halo term to the 1-halo term is smooth.
This transition was known to be a problematic region to model accurately for
the standard halo model prescription. In this case, the similarities in the
mass functions for the two cosmologies, i.e., since the two cosmologies have
the same linear clustering, allows the inaccuracy in the transition region to
be overcome. This was first shown in [20] and later validated for $f(R)$ and
DGP gravity in [21]. This idea was implemented into a numerical fitting
function dubbed the “halo model reaction” or ReACT [19, 17, 22].

A key point to note is that in the publicly available implementation of ReACT,
two separate parameters, $\mu$ and $\mathcal{F}$ are used to model the linear
and non-linear gravitational clustering, respectively. The linear
modifications to the Poisson equation are encoded in $\mu(a,k)$, while
$\mathcal{F}$ is computed by assuming a specific model and applying the
spherical collapse formalism to specific modified gravity models with
specified screening mechanisms. Given this prescription, the ReACT code
computes the reaction from the following input information:

  * •

The linear power spectrum at the redshift(s) of interest for the target
cosmology.

  * •

The function $\mu$ which encodes the linear modification to gravity via the
Poisson equation applied to scales over which perturbation theory may be
applied.

  * •

The function $\mathcal{F}$ that encodes the modification to the non-linear
Poisson equation that is relevant for $N$-body simulations.

  * •

Additional input parameters connected to specific models of modified
gravity/established linear theory paramterisations such as the Effective Field
Theory (EFT) of dark energy.

However, in our parameterisation the Poisson equation is derived making no
distinction between the clustering of matter due to gravity on linear vs non-
linear scales (see sec. IV A in [13]). Therefore, in order to make contact
with our simulations within the the ReACT framework, one sets
$\mathcal{F}\equiv\mu(a,k)-1$. In other words, we ensure that the parameter
$\mu(a,k)$ specifies both the linear and non-linear clustering in our
simulations, which in general is not true for many models considered in the
basic ReACT framework. We also switch off/set to their respective $\Lambda$CDM
values any additional model-specific parameters.

The corrections from the modifications to the Poisson equation are propagated
through to the halo model quantities such as the mass function via the
spherical collapse model implemented within ReACT. This introduces a slightly
different source term in the differential equation for the radius of the halo
within spherical collapse. We describe this procedure in appendix A.

We quantify the performance of ReACT in later sections by comparing the matter
power spectrum $P(k)$ measured from our simulations against that predicted by
the fitting function. In practice, we will compute the so-called halo model
‘reaction’ term, a dimensionless ratio of halo model quantities given by

$R(z,k)=\frac{P(z,k)}{P_{\rm pseudo}(z,k)}$ (2.3)

where $P(z,k)$ is the matter power spectrum in the specified cosmology and
$P_{\rm pseudo}(z,k)$ is the $\Lambda$CDM matter power spectrum with the
initial conditions modified such that the linear growth factor is identical to
the the numerator at the redshift of interest. Using eq. 2.3, we define the
$k_{\rm fail}$ parameter as the smallest wavenumber at which ReACT failed to
predict $R(k)$ in our simulations at a specified accuracy threshold of 1%.

###  2.2 Modifying halo model ingredients in ReACT

We found evidence in [14] to suggest that the accuracy of ReACT relative to
our simulations may be improved by modifying the halo model parameters within
the code. In particular, we modified the halo concentration-mass amplitude
parameter. In this work we systematically develop this finding into a
modification to ReACT that substantially improves the ‘vanilla’ ReACT
approach. In this subsection we describe this procedure to model $c(M)$ as a
function of $\mu(z)$ to maximise the accuracy of ReACT relative to our
simulations. We call this implementation the ‘extended’ implementation of
ReACT.

A key ingredient in the halo-model approach in ReACT is the density profile of
dark matter halos. The standard profile is the Navarro-Frenk-White (NFW)
density profile [23]

$\rho(r)=\frac{\rho_{\rm s}}{\frac{r}{r_{\rm s}}\left(1+\frac{r}{r_{\rm
s}}\right)}\,,$ (2.4)

where $r_{\rm s}$ and $\rho_{\rm s}$ are the scale radius and density,
respectively. These are both free parameters within the simple
parameterization. In general, the NFW profile is divergent, but one can impose
an upper limit on the size of the halo given by the virial radius $r=r_{\rm
vir}$. Integrating the density profile gives the halo mass given by

$M=\frac{4\pi\rho_{\rm s}r_{\rm
s}^{3}}{3}\left(\ln(1+c)-\frac{c}{1+c}\right)\,,$ (2.5)

where the NFW concentration parameter $c=r_{\rm vir}/r_{\rm s}$. The
parameterization is usually complemented by a mass-concentration relation,
$c(M)$. This relation is measured from simulations, and it reduces the number
of free parameters to one that is defined by the total mass of the halo given
a specific background density.

Within the vanilla implementation of the ReACT algorithm, the $c(M)$ relation
is specified by

$c(M)=\frac{c_{0}}{1+z}\left(\frac{M}{M_{*}}\right)^{-\alpha}\,,$ (2.6)

where $c_{0}=9.0$, $M_{*}$ is the characteristic mass computed by imposing a
condition on the peak height $\nu(M_{*})=1$ of the Sheth-Tormen mass function
and $\alpha=0.13$.

The $c(M)$ relation for the $\Lambda$CDM model, has been measured from a suite
of $N$-body simulations [24] and is parameterized by

$c(M)=A\left(\frac{M}{10^{12}h^{-1}M_{\odot}}\right)^{-\gamma}\left[1+\left(\frac{M}{M_{0}}\right)^{0.4}\right]\,,$
(2.7)

where $A_{\rm\Lambda CDM}=10.2,\gamma_{\rm\Lambda
CDM}=0.1,M_{0}=10^{17}h^{-1}\,{\rm M_{\odot}}$. We compare the concentration-
mass relation above for Planck 2018 (TT+TE+low P+BAO) [25] parameters with the
$c(M)$ measured from our simulations and the $c(M)$ relation in ReACT in fig.
1. Moreover, we compare (2.6) and (2.7). Although we find that our simulations
are quite noisy - they have not been optimised to measure the $c(M)$ relation
- all are in agreement with each other. In what follows we will use (2.7) to
modify the $c(M)$ relation within ReACT, keeping $\gamma$ and $M_{0}$ constant
and varying $A$.

Figure 1: In the left panel, we present the concentration-mass relation
measured from our $\Lambda$CDM simulations and a fit for the Planck 2018
(TT+TE+low P+BAO) [25] parameters from ([24]). We note that there is
significant scatter as our simulation suite has not been optimised for
measuring $c(M)$. In the right panel, we illustrate that the ReACT $c(M)$
using the expression in (2.7) agrees with that from (2.6). This allows us to
express the concentration-mass relation in our modified gravity simulations as
a ratio relative to the $\Lambda$CDM ReACT case.

In general, the concentration scales with mass as a weak power law at low
redshifts. This trend can be traced back to the hierarchical nature of
structure formation. At late times, when smaller halos have already formed and
are merging to form larger halos, the inner regions of these large halos are
relatively undisturbed by mass accretion. Therefore, as halos become larger,
their virial radii become bigger while the core radius remains relatively
unchanged. This then leads to more massive halos having lower concentrations
in comparison to less massive halos [24]. As such, hierarchical structure
formation leads one to expect that increasing $\mu$ at late redshifts will
enhance structure formation at late times, affecting massive halos and thus
reducing the ratio of $A/A_{\rm\Lambda CDM}$ and vice-versa for when $\mu$ is
increased at early times. We found evidence for this behaviour in [14] and in
this work. In addition, we also find additional complexity arising at very low
$z_{\rm mp}$ that we couldn’t have resolved earlier due to the larger width of
our $\mu$ bins relative to this work.

In this work, we show that the computation of the full non-linear $P(k)$ with
high accuracy in ReACT for generic modified gravity requires modification of
the concentration-mass relation in section 4.2. We do this by finding the
value of $A$ in (2.7) that gives the best representation of the non-linear
power spectrum from the simulation, where this is defined as the
representation that has the largest $k_{\rm fail}$. By performing a range of
simulations with different modifications to gravity, we provide a fitting
function for the amplitude $A$ of the $c(M)$ relation as a function of
$\mu(z)$ validated against the simulations.

We primarily use an accuracy threshold of 1% to compute $k_{\rm fail}$.
However, we note that baryonic feedback has been studied in a recent release
of ReACT [18] quoting an accuracy of 3%. Therefore, we also present, both in
our abstract and conclusions, the 1% and 3% $k_{\rm fail}$ values. This is
particularly relevant for our weak lensing results, where we also compute the
corresponding values of $\ell_{\rm fail}$ for the weak-lensing convergence
power spectrum (again for 1% and 3% accuracy thresholds), and we reflect on
how robust these numbers are across our parameter space (see sec. 4.3).

We have optimised our simulation suite to investigate and fit the non-linear
matter power spectrum, which is sub optimal for measuring the $c(M)$ relation
from the simulations. As such, we cannot perform a consistency check between
the $c(M)$ in our modified gravity simulations and our fitting function
prediction. Since we are only interested in a accurate way of creating the
non-linear power spectrum from the linear one, this test is not essential, and
we leave it to future work to examine this in detail.

We note that modifications to gravity have been shown to affect the
$\Lambda$CDM $c(M)$ relation. It has been shown that in cosmologies with
modified growth histories such as smooth dark energy models that the $c(M)$
relation is modified in order to accurately predict $P(k)$ using the halo
model [26]. In $f(R)$ gravity, it has been shown that the $c(M)$ is modified
from a simple power law due to the fact that $M_{*}$ contains implicit
environmental dependence on halo mass [27, 28, 29]. These findings further
justify our choice to improve the fit from “vanilla” ReACT by modifying the
$c(M)$ relation.

##  3 N-body simulation suite & coverage of parameter space

In [14] we studied the phenomenology of models with a time-dependent $\mu$ by
implementing piece-wise constant bins777Up until section 4.3, all references
to “bins” refer to the piecewise constant parts of the $\mu$ function, not to
tomographic weak-lensing analysis bins. in redshift for $\mu$ in our
simulations. In particular, we focused on the matter power spectrum $P(k)$ in
our simulations. We showed that non-linear clustering statistics can exhibit
significant variation as a function of $\mu(z)$ for the same linear $P(k)$.
This analysis revealed that the ReACT formalism was able to accurately predict
$P(k)$ in our simulations, albeit for a small region of parameter space. In
this work we quantitatively validate the ReACT formalism to understand the
range of scales over which it can be used to predict $P(k)$ with specified
accuracy and precision. To do this, we need to explore a wide range of
modified gravity parameter space. Since we bin $\mu$ in redshift, this
involves sampling a three dimensional parameter space where the degrees of
freedom are the centre of the bin, the bin width, and the value of $\mu$.

As shown in [14], the linear growth factor $D(z)$ between the start and end of
a redshift bin is a good quantity to use to characterise the bin width (rather
than working directly with a difference in $a$ or $z$). As such, we continue
to work with redshift bins of constant $D(z)$ in $\Lambda$CDM (this can be
computed analytically in linear theory). We now describe the procedure we
followed to fix $\mu(z)$ across our simulations.

We designed two sets of simulations to extensively test the performance of
ReACT. In simulation set $\alpha$ (presented in table 1) we fixed $D(z)$ and
obtained 8 redshift bins. All the bins in these simulations are chosen such
that the $\Lambda$CDM growth factor associated to them is constant and is
$D(z_{\rm max})/D(z_{\rm min})=1.26$. We choose the initial value of $D(z)$ to
be similar to that associated with the binning in [30] so that we may
contextualise our results (and any future results derived from the data
products in this work) with parameterised Euclid forecasts.

In set $\alpha$ we pick two different values of $\mu$ on each side of $\mu=1$,
with the expectation from our findings in [14] that our results are likely to
be symmetric about the line $\mu=1$. In other words, the matter power spectrum
exhibited equal and opposite effects when $\mu=1.1$ and $\mu=0.9$,
respectively in the same redshift bin. We found both in our previous work and
in this analysis that this behaviour persists and is independent of the choice
of $\mu(z)$, (see fig. 3 for more info).

We also noted in [14] that the performance of ReACT varied with position and
width of the bin in question. In particular, we found that the ‘maximal’
departure from $\Lambda$CDM was in redshift bins at the extreme end of the
redshift range associated to our simulations. This motivates a higher
resolution analysis at low $z$ in set $\beta$ 888Ideally, one would also like
to sample high redshifts in the simulations. However, we find that the ReACT
formalism fails to accurately reproduce $P(k)$ at $z>7$ across the range of
values of $D(z)$ that we explored in our simulation suite. Exploration of
ReACT at $z>7$ is left for future work. It is pertinent to point out that the
range of scales that are non-linear at $z>7$ are too small to be relevant in
this analysis.

In simulation set $\beta$, we vary of $\mu,D(z)$ and over a range values to
further validate ReACT in our approach. We restrict ourselves to $\mu>1$
(taking advantage of the aforementioned symmetry to save on computational
resources) and introduce additional bins at low redshift (see fig. 2). We
implement finer sampling of $\mu$ at the lower range of $z$ in set $\beta$,
which reveals that the shape of the $P(k)$ undergoes a dramatic change as the
midpoint of the redshift bin approaches $z=0$ (see figs. 3, 4 and associated
discussion). In order to develop and validate our fitting function for the
concentration-mass amplitude parameter $A$, we ran additional simulations in
which we vary the linear growth factor associated with our redshift bins;
again, we restrict ourselves to $\mu>1$. We illustrate the full range of bins
we have run in fig. 2.

We compare simulations across different values of $\mu$ and $D(z)$ by keeping
the bin midpoints constant across the parameter-space. In other words, the
endpoints of each redshift bin are varied in order to keep the midpoint the
same for all $D(z)$. We will return to this point in the next subsection where
we show our fitting function. This procedure allows us to systematically
investigate the performance of ReACT for varying bin-widths, epochs and
$\mu(z)$ in as orthogonal a manner as feasibly possible. Furthermore, as we
will see in section 4, keeping the midpoint constant allows a simple
interpretation of the variation of $A(\mu,D(z))$. Our full suite of
simulations represents an unprecedented exploration of phenomenological
modified gravity parameter space on non-linear scales.

$\Delta z_{\rm bin}$ | $\mu_{1}$ | $\mu_{2}$ | $\mu_{3}$ | $\mu_{4}$ | $\sigma^{1}_{8}$ | $\sigma^{2}_{8}$ | $\sigma^{3}_{8}$ | $\sigma^{4}_{8}$
---|---|---|---|---|---|---|---|---
$7\geq z\geq 5.4$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7896 | 0.8002 | 0.8215 | 0.8321
$5.4\geq z\geq 4.1$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7896 | 0.8002 | 0.8216 | 0.8322
$4.1\geq z\geq 3.0$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7889 | 0.7999 | 0.8218 | 0.8329
$3.0\geq z\geq 2.1$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7900 | 0.8004 | 0.8213 | 0.8318
$2.1\geq z\geq 1.47$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7915 | 0.8012 | 0.8205 | 0.8302
$1.47\geq z\geq 0.91$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7937 | 0.8023 | 0.8194 | 0.8281
$0.91\geq z\geq 0.43$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.7977 | 0.8043 | 0.8175 | 0.8241
$0.43\geq z\geq 0$ | 0.8 | 0.9 | 1.1 | 1.2 | 0.8052 | 0.8081 | 0.8136 | 0.8165
Table 1: The simulation parameters for set $\alpha$ of our simulation suite in
order to validate ReACT. This choice was designed to mimic the binning in [30]
which is likely to be representative of the bin-width typically chosen for a
Euclid-like survey. For this set of simulations, we varied $\mu$ symmetrically
about $\mu=1$, to test if the symmetry we found in [14] persists. Note that
these are all bins with equal $\Lambda$CDM growth, i.e.,
$D(z_{i-1})/D(z_{i})=D(z_{i})/D(z_{i+1})=\text{constant}=1.26$. This set does
not represent the full suite of simulations we run (see fig. 2). Due to the
symmetry in the shape of the power spectrum, we only use $\mu>1$ for set
$\beta$, to make optimal use of computational resources.  Figure 2: Locations
of the set $\beta$ simulations in the $\Delta z-z_{\rm mp}$ plane: the width
of each redshift bin $\Delta z$ as a function of the bin midpoint $z_{\rm mp}$
for the different values of $D(z)$ we used in our simulations. We present the
redshift bins across all of our simulations and we see that for $z>1$ all the
bin midpoints coincide. This allows us to use $z_{\rm mp}$ as a reference
point to compare results across different values of $D(z)$. This figure
contains both set $\alpha$ (red data-points) (see table 1) and set $\beta$
(blue and brown data points). In set $\beta$, we varied $D(z)$ and increased
our redshift resolution at $z<1$ to sample the resulting change in shape in
$P(k)$ for these low-$z_{\rm mp}$ bins. Note that we run simulations with
$\mu=\\{1.1,1.15,1.2\\}$ for all these bins.

##  4 Results

###  4.1 Vanilla Implementation

In the previous sections, we have explained the reasons that motivated our
choice of using the $\mu-D(z)$ parameter space to explore the modified gravity
phenomenology in our simulations. In this section, we show the results of
measuring the matter power spectrum from our suite of simulations and quantify
the performance of ReACT, both without any modifications and with
modifications to the halo concentrations. We present $R(k)$ for the different
$\mu(z)$ presented in table 1 in fig. 3. We note that the symmetry about
$\mu=1$ is preserved irrespective of $D(z)$, $\mu(z)$ and the midpoint of the
redshift bin. As a consequence of this symmetry, we restrict ourselves to
$\mu>1$ for the rest of our analysis.

Figure 3: $R(k)$ at $z=0$ from all of our simulations with redshift bins
associated with $\Lambda$CDM linear growth $D=1.26$. We observe symmetry in
the power spectra about $\mu=1$. We also note that the extent to which the
power spectra differ from $\Lambda$CDM increases as the redshift at which
modified gravity is switched on decreases. Consistently with this, we see
later that the accuracy of ReACT is better for the power spectra closer to
$\Lambda$CDM, i.e. for smaller $|\mu-1|$ and higher $z_{\rm mp}$.

We also observe in fig. 3 that the ‘departure’ from $\Lambda$CDM seems to be
more dramatic the later $\mu$ is switched on. This behaviour is consistent
with our earlier results in [14] (but now with greater redshift resolution in
the transition regime), where the relative difference to $\Lambda$CDM
decreases with increase in the redshift at which the modification occurs.
Furthermore, the shape of the matter power spectrum changes drastically when
$\mu$ is switched on in the bin $0\leq z\leq 0.43$ compared to the general
trend seen in the other bins. This regime is examined in greater detail in set
$\beta$ see fig. 2). We present the matter power spectra from this set of
simulations in fig. 4, where it is clear that this transition seems to occur
around the epoch where the transition to dark energy domination takes place.
We will return to this point shortly.

In [14] we hypothesised that the shape of the power spectra can be traced back
to halo formation time as a function of its mass. In other words, due to the
hierarchical nature of structure formation, the average mass of the halos
forming in our simulations is large at low redshift. Therefore, the range of
scales affected relative to the opposite case are when $\mu$ is modified at
low redshift is larger. In other words, there is a transfer of power (from
small scales to large scales) that takes place when $\mu$ is switched on at
late times. We will come back to this point when we discuss our fitting
function for halo concentrations, since this tends to reduce the typical
concentration of the halos.

In fig. 5, we quantify this point by computing the value of $k_{\rm fail}$ for
our simulations in set $\beta$. We find that the departure from $\Lambda$CDM
and, therefore, the accuracy of ReACT is monotonic with $\mu$. We don’t see
such a clear trend for $k_{\rm fail}$ as a function of $D(z)$. Interestingly,
ReACT does much better for the simulations in which $\mu$ is modified at
$2\leq z\leq 7$, across parameter space. As mentioned earlier (see sec. 2.2,
the additional complexity in structure formation arising from mergers of halos
drives down $k_{\rm fail}$ for the lower redshift bins, with exceptions
arising in cases where the bin-midpoint is such that the bin straddles the
transition period between the two separate regimes of structure
formation.999This is due to an artificial cancellation of errors.

Figure 4: $R(k)$ from the simulations with the $\mu$ modified in the redshift
range containing the transition from matter to dark energy domination. The
general trend indicated in fig. 3 is clearly broken once the Universe enters
dark energy domination. This fact will influence the performance of our
fitting function (see fig. 6). In the context of maximising the constraining
power of future surveys, ensuring that $k_{\rm fail}>1\,h\,{\rm Mpc}^{-1}$
thus requires modifying the halo concentrations consistently with $\mu(z)$
(see sec. 2).

We note that the lowest $k_{\rm fail}$ is associated to the smallest redshift
bin, i.e., $0\leq z\leq 0.43$. As mentioned before, the transition into dark
energy domination seems to coincide with the change in shape of the power
spectra at $z<0.5$ as we previously mentioned (see fig. 4), we see that the
trend of how $P(k)$ is modified as a function of $\mu(z)$ is broken and indeed
reversed in the range $0\leq z\leq 0.9$. Due to the accelerated expansion at
these late times, structure formation appears to be ‘frozen in’, and
increasing $\mu$ simply adds power much more uniformly across a range of
quasi-linear to non-linear scales relative to the pseudo-$\Lambda$CDM case as
seen in fig. 4.

From our analysis, we see that the main challenge in predicting $P(k)$ across
our simulations is ensuring a consistent value of $k_{\rm fail}$ across all
our simulations. Particularly in the context of maximising the constraining
power of future surveys, ensuring that $k_{\rm fail}>1\,h\,{\rm Mpc}^{-1}$
would require a factor $>2$ improvement in $k_{\rm fail}$ for the lowest
redshift bins. We have already mentioned in sec. 2 that the ReACT $P(k)$
prediction can be improved by modifying the halo concentrations consistently
with $\mu(z)$. We will now describe how this can be achieved.

Figure 5: The $k_{\rm fail}$ values as a function of $z_{\rm mp}$ i.e., the
midpoint of the redshift bin at which $\mu$ is ‘switched on’ for our vanilla
implementation of ReACT. We see that, as indicated in fig. 3, the accuracy of
ReACT when $\mu$ is closer to GR, i.e., $k_{\rm fail}$ monotonically decreases
as $\mu$ increases. We note that $k_{\rm fail}$ tends to increase with $z_{\rm
mp}$, but this trend is broken when $z_{\rm mp}$ is in the transition regime
from matter domination to dark energy domination.

###  4.2 Extended Implementation

We perform a fit for the ratio of amplitude of the $c(M)$ parameter relative
to the $\Lambda$CDM case, that is $A/A_{\Lambda{\rm CDM}}$, required to
maximise the $k_{\rm fail}$ for all $\mu(z)$. This allows us to capture the
varying phenomenology of our simulations more accurately using ReACT. Note
that while we are fitting for a single parameter, i.e., the ratio of the
amplitudes of $c(M)$, the actual fit contains multiple parameters (see eq.
(4.4)). Following the ReACT philosophy of computing quantities as a ratio
relative to $\Lambda$CDM, we compute our fitting function for the ratio
$A/A_{\rm\Lambda CDM}$, thus ensuring our analysis is insensitive to
realisation-dependent effects [31]. We remind the reader that in order to
compute $R(k)$, one needs the modified gravity $P(k)$ as well as the
$\Lambda$CDM $P(k)$ with equalised linear growth.

As earlier, we use the $k_{\rm fail}$ parameter to quantify the accuracy with
which ReACT captures the simulation results. We vary the amplitude $A$ of the
$c(M)$ relation to maximise $k_{\rm fail}$. We then calculate the ratio
$A/A_{\rm\Lambda CDM}$ as a function of the bin midpoint for each set of
simulations in fig. 2. In other words, we model the parameter space across the
range of values of $1\leq\mu\leq 1.2$ and $1.26\leq D\leq 1.45$. With the
variation of the amplitude of the concentration-mass relation, we are able to
achieve $k_{\rm fail}>1.5\,h^{-1}{\rm Mpc}$ for an accuracy threshold of $1\%$
and $k_{\rm fail}>3\,h\,{\rm Mpc}^{-1}$ for an accuracy threshold of $3\%$.

We wish to create a fit for the concentration modification as a function of
bin midpoint, and do this as a function of $\mu$ and $D$, thus covering the
full phenomenological parameter space we are considering. The general
functional form for the fit is given by

$A/A_{\Lambda\rm CDM}=\begin{cases}f_{1}(\bar{z})&\text{if }\bar{z}\geq 1.2\\\
f_{2}(\bar{z}),&\text{if }\bar{z}\leq 1.2\end{cases}$ (4.1)

where $\bar{z}=z_{\rm mp}-z_{\rm pk}$ and $z_{\rm pk}$ is the redshift at
which the power spectrum is computed. We divide the redshift range into two
regimes. In the first regime, i.e., $1.2\leq\bar{z}\leq 7$, we use an
exponential function

$\displaystyle f_{1}(\bar{z})=C_{1}e^{-\bar{z}}+1\,,$ (4.2)

where $C_{1}<0$ when $\mu>1$ (and vice-versa for $\mu<1$). This form ties back
to the qualitative prediction that we made from hierarchical structure
formation, i.e., modifying $\mu$ at later times leads to larger halos and
therefore smaller concentrations relative $\Lambda$CDM. Note that for
calculating the power spectrum at $z_{\rm pk}$ when $\mu\neq 1$ in a
particular bin, $A/A_{\rm\Lambda CDM}$ has a single constant value that is a
function of $z_{\rm pk}$ and the bin parameters $\\{z_{mp},\mu,D\\}$ (see
appendix C for details). In particular, $A/A_{\rm\Lambda CDM}$ is not taken as
a function of redshift during the integration in redshift that occurs in the
halo model calculation within ReACT. Our general fit for $z_{\rm pk}$ could be
used to make $A/A_{\rm\Lambda CDM}$ a function of the redshift in the halo
model calculation; we leave an investigation of whether this further improves
the accuracy to future work.

We find that this trend is broken at very low redshift, i.e., at $z_{\rm
mp}<1$. We fit this upturn in the behaviour of $A/A_{\rm\Lambda CDM}$ using a
cubic polynomial

$f_{2}(z)=(C_{2}+C_{3}\bar{z}+C_{4}\bar{z}^{2}+C_{5}\bar{z}^{3})\,.$ (4.3)

In the transition region between these two cases ($0.7\leq z_{\rm mp}\leq
1.2$) we use a cubic spline to interpolate between $f_{1}$ and $f_{2}$. Thus,
we have a list of parameters $C_{i}=\\{C_{1},C_{2},C_{3},C_{4},C_{5}\\}$ that
vary as a function of $\\{\mu,D\\}$.

We find that the $C_{i}$ scale linearly with $D$ at fixed values of $\mu$, and
the dependence on $\mu$ is a second order polynomial at fixed values of $D$
(see appendix C for more details). We consider $\mu$ first and derive the
fitting parameters $\tilde{C}_{i}$ and $\bar{C}_{i}$, as a function of
arbitrary $\mu$ at the extreme values of $D=\\{1.26,1.35\\}$, respectively. We
then linearly interpolate in $D$ to obtain the general form $\hat{C}_{i}$ in
the two-dimensional $(\mu,D)$ space. This allows us to explicitly construct
the final fitting function for $A$ at $z_{\rm pk}=0$ to be

$\displaystyle x$ $\displaystyle=$ $\displaystyle 1-(1.45-D)/(0.19)\,,$ (4.4)
$\displaystyle\hat{C}_{i}$ $\displaystyle=$
$\displaystyle\tilde{C}_{i}(1-x)+\bar{C}_{i}x\,,$ (4.5) $\displaystyle
f_{1}(\bar{z})$ $\displaystyle=$ $\displaystyle\hat{C}_{1}e^{-\bar{z}}+1\,,$
(4.6) $\displaystyle f_{2}(\bar{z})$ $\displaystyle=$
$\displaystyle\hat{C}_{2}+\hat{C}_{3}\bar{z}+\hat{C}_{4}\bar{z}^{2}+\hat{C}_{5}\bar{z}^{3}\,,$
(4.7)

where we present the explicit functional forms for $\tilde{C}_{i}$ and
$\bar{C}_{i}$ in appendix C. Note that equation gives an accurate calculation
of the non-linear matter power spectrum at z=0 for arbitrary $\mu$, $D(z)$ in
the range $0\leq|\mu-1|\leq 0.2$ and $1.26\leq D(z)\leq 1.45$. See appendix C
for how this is generalised to power spectra at non-zero z and for $\mu<1$.
This general form is implemented in our extended ReACT code.

We show the performance of the fitting function described in this section at
$z=0$ as a function of $\mu$ and $D(z)$ in fig.6. The two redshift regimes are
clearly visible in both the panels. Each data point is obtained from the
simulations, and the $k_{\rm fail}$ values associated with this ‘extended’
implementation can be seen in fig. 10. The robustness of this fitting function
to simulation and cosmological parameters is discussed in appendix B. We also
present a qualitative justification of the general observation of the
modification at low redshift being more substantial than when $\mu$ is
modified at high redshift in appendix A.

Interestingly, the sharp upturn in $A$ coincides with the redshift at which
$\Omega_{\Lambda}$ starts becoming comparable to the critical density of the
Universe. Once dark energy domination is firmly established, structure
formation is frozen in and each modification of gravity affects structure
formation uniformly across the range of scales of convergence in our
simulations. This numerical recipe is to our knowledge the only validated one
in the literature for model-independent modified gravity.

We have already presented the 1% $k_{\rm fail}$ for all of our simulations in
fig. 10. In particular, we present the $k_{\rm fail}$ for varying accuracy
thresholds in appendix C. We find that across all of our simulations, we
achieve a 1% $k_{\rm fail}\geq 2\,h\,{\rm Mpc}^{-1}$. Note that this is
considerably improved for the bins with $z_{\rm mp}>1$, where $k_{\rm
fail}>3\,h\,{\rm Mpc}^{-1}$ can be achieved. Most notably, we see that a
factor $\gtrsim 2$ improvement in $k_{\rm fail}$ with our concentration
modification compared to the vanilla ReACT. Finally, $k_{\rm fail}$ values of
$\gtrsim 3\,h\,{\rm Mpc}^{-1}$ can be achieved to within 3% consistently
across all our simulations.

Figure 6: The fitting function for the ratio of the amplitude of the
concentration-mass parameter relative to $\Lambda$CDM as a function of $z_{\rm
mp}$. The data points are computed from our simulations at the midpoint of the
redshift bin at which $\mu$ is switched on in the simulation. In the left
panel, we present our fit for varying $D(z)$ and in the right panel, for
varying $\mu$. We see that the shape of the fit remains consistent across
parameter space, which ties back to our discussion regarding the two regimes
(see sec. 2). The exponential functional form at $z>1$ follows from the
hierarchical structure formation (see appendix A) while the strong upturn in
the concentration at low redshifts is likely due to the freezing in of
structure formation as dark energy domination occurs.

###  4.3 Quantifying the performance of the fit for lensing observations

Figure 7: We show the ratio of the weak-lensing convergence spectrum computed
using the power spectra directly measured from our simulations to that
computed from ReACT. The vertical lines are the cut-off in $\ell$ where the
disagreement exceeds 1% for each case. In the left panel, we present the most
pessimistic scenario, i.e., the $\mu(z)$ corresponding to the lowest value of
$k_{\rm fail}$. In the right panel, we present a more ‘realistic’ scenario
with a $k_{\rm fail}$ is around the typical value in our simulations. For
larger values of $\mu$ the improvement is from $\ell\sim 200$ to $\ell\sim
500-1000$. For smaller $\mu$ the improvement is from $\ell\sim 400$ to
$\ell>1000$.  Figure 8: The weak-lensing convergence power spectrum cut-off
computed for the bin with the smallest $k_{\rm fail}$ compared to the
$\Lambda$CDM spectrum with Euclid-like error bars. As presented in table 2, we
show the cutoff that DES employed in their beyond-$\Lambda$CDM analysis which
was restricted to linear scales, 1% and 3% values of $\ell_{\rm fail}$ for
$\mu=1.2$ (red) and $\mu=1.1$ (blue). We also indicate the regions where
baryons are most likely to play a significant role [32]. Here, we see
explicitly here that the amount of information scales as $\ell^{2}$, and
therefore cutting off the $C_{\ell}$ at low value of $\ell$ leads to
information loss and therefore smaller constraining power. Note also that the
error bars become smaller as one goes from small $\ell$ to large $\ell$, i.e.,
from large scales to small scales.

We now quantify the accuracy and range of validity of our fit for future
lensing surveys and deduce $\ell_{\rm fail}$, analogous to $k_{\rm fail}$ in
angular space101010Note that the angular projection means that there is not a
one-to-one correspondence between $k_{\rm fail}$ and $\ell_{\rm fail}$. We
compute the weak-lensing convergence power spectrum from the matter power
spectrum and the modified gravity parameters

$C_{\ell}=\frac{9H_{0}^{4}\Omega_{\rm m}^{2}}{4c^{4}}\int_{0}^{\chi_{\rm
max}}\mathrm{d}\chi\frac{1}{a^{2}(\chi)}g^{2}(z)\frac{\mu^{2}(1+\eta)^{2}}{4}P_{\delta}(\ell/\chi)\,,$
(4.8)

where $P_{\delta}$ is the matter power spectrum, $\eta$ is the second modified
gravity parameter that affects the photon geodesics and $\chi$ is the comoving
angular diameter distance to the source along the line of sight. We do this
for the two conservative and realistic scenarios from our simulations.

The function $g(z)$ is a filter function that depends on the redshift
distribution of the background galaxies and is typically written as

$g_{i}(z)=\int_{z}^{\infty}\mathrm{d}z^{\prime}\left(1-\frac{\chi(z)}{\chi(z^{\prime})}\right)n_{i}(z^{\prime})\,.$
(4.9)

We use the standard expression for the source galaxy redshift distribution
given by [33]

$n(z)\propto\left(\frac{z}{z_{0}}\right)^{\alpha}\exp\left[-\left(\frac{z}{z_{0}}\right)^{\beta}\right]\,.$
(4.10)

We adopt a Euclid-like binning of the source number density into 10 equi-
populated bins according to eq. (4.10) with $z_{0}=0.9/\sqrt{2}$, $\alpha=2$
and $\beta=3/2$, where we have assumed an average source density $\bar{n}_{\rm
g}=30\,{\rm arcmin}^{-2}$ [34, 35, 30, 36, 37]. The tomographic bins are such
that our fitting function may be used across the entire weak-lensing data-set.
The projected errorbars are obtained by computing the following [37]

$\delta C_{\ell}^{ij}=\sqrt{\frac{2}{f_{\rm
sky}(2\ell+1)}}\left(C_{l}^{ij}+\frac{\sigma_{\epsilon}^{2}}{\bar{n}_{i}}\delta^{ij}\right)\,,$
(4.11)

where $f_{\rm sky}=0.7$ is the fractional sky coverage,
$\sigma_{\epsilon}=0.21$ is the variance of the observed ellipticities and
$\bar{n}_{i}$ is the surface galaxy density of each bin.

We compute the $C_{\ell}$ from the matter power spectra using both our
simulation snapshots as well as from ReACT. We compare the prediction from
vanilla ReACT and with our concentration fitting function implemented in the
code with that obtained from the simulations to obtain an $\ell_{\rm fail}$
value up to which we formally validate ReACT for weak-lensing.

While we cannot associate a one-to-one mapping between $k_{\rm fail}$ and
$\ell_{\rm fail}$, we would expect naively that simulations with small values
of $k_{\rm fail}$ will be correlated to the ones with low $\ell_{\rm fail}$
values. We report the 1% $\ell_{\rm fail}$ value for this case as our most
conservative estimate of the scale up to which one can reliably compute the
convergence power spectrum $C_{\ell}$. However, we also quote the 3%
$\ell_{\rm fail}$ as a more ‘realistic’ estimate of the physical scales that
can be accurately modelled with ReACT. Note that our 1% accuracy results may
be viewed to be conservative because the expected level of accuracy that can
be achieved across the parameter space validated in our simulations accounting
for baryonic feedback [17] in ReACT is 3%. We show an example of these two
cases in fig. 7. Note that applying our extended implementation of ReACT
always leads to a significant improvement in the range of scales we can
accurately model. In addition, For smaller values of $\mu$, even our vanilla
implementation of ReACT is able to achieve sub-3% accuracy.

Table 2 shows that with our vanilla implementation, we are able to
consistently model the weak-lensing convergence power spectrum only up to
$\ell\sim 500$. In comparison, with our extended implementation, we are able
to reliably model up to $\ell\sim 1300$. To put this in context, the DES
analysis that reported 3x2pt constraints [9] on extensions to $\Lambda$CDM
employed stringent scale-cuts that only data on linear scales survived. This
corresponds to an $\ell_{\rm fail}\sim 100$ in our analysis. We note that this
leads to a significant loss in constraining power. In a recent study [32], it
was found that the modelling of baryons limited the constraining power that
can be gained on scales corresponding to $\ell>1500$. We condense all these
results into fig. 8, where it is clear the that using our extended formalism
allows one to achieve values of $\ell_{\rm fail}$ close to the limit quoted
above which baryons will likely dominate the constraining power. We note that
we are able to achieve $\ell_{\rm fail}$ within a factor of $\sim 2$ of this
limit even in the most conservative cases, and indeed achieve $\ell_{\rm
fail}>1500$ for the vast majority of our simulations. Without our
concentration fitting function, we consistently achieve $\ell_{\rm fail}<500$.

$|\mu-1|$ | Vanilla $\ell_{\rm fail}$ | Extended $\ell_{\rm fail}$ | Vanilla $k_{\rm fail}$ | Extended $k_{\rm fail}$
---|---|---|---|---
1% | 3% | 1% | 3% | 1% | 3% | 1% | 3%
0.2 | 400 | 900 | 900 | 1200 | 0.4 | 0.8 | 1.2 | 3.0
0.1 | 600 | 2100 | 2400 | 3500 | 0.7 | 1.0 | 1.4 | 4.3
Table 2: We report the value of $\ell_{\rm fail}$ that we find across the
range of values of $\mu$ in our simulations. We find that the variation in
$\ell_{\rm fail}$ is much more significant with $\mu$ compared to the linear
growth factor $D(z)$. Therefore, we report a conservative 1% $\ell_{\rm fail}$
corresponding to the simulation with the smallest 1% $k_{\rm fail}$ and a
‘realistic’ 3% $\ell_{\rm fail}$ that corresponds to a more typical value of
$k_{\rm fail}$ that can be achieved to within 3% (see appendix C for more
details on $k_{\rm fail}$ as a function of $\mu$ and accuracy threshold). We
pick these numbers to correspond to the bin corresponding to the smallest
values of $k_{\rm fail}$ as the limiting value for given value of $\mu$, which
is $0.91\leq z\leq 1.47$.

##  5 Conclusion

The wealth of data provided by upcoming surveys such as the Euclid satellite,
NRT, VRO and SKA creates an unprecedented opportunity to test modified gravity
in a model-independent way with cosmology, which could provide crucial
evidence for the nature of the dark sector. The majority of this survey data
is on non-linear scales, so robust theoretical predictions for this regime are
crucial for maximising the gain from these surveys.

In this work we extend the preliminary investigation in [14] to perform a
thorough validation of our ReACT approach for predicting the non-linear matter
power spectrum in phenomenological modified gravity. We cover the modified
gravity parameter space by treating $\mu$ as a piece-wise constant function,
and we create a substantial suite of N-body simulations varying the value in
each piece, the width of each piece, and the location (in redshift) of each
piece.

We analyse the performance of our “vanilla” implementation of ReACT across
this simulation suite, finding that we are able to model the weak lensing
convergence power spectrum measured in our simulations across parameter space
reliably up to $\ell\sim 500$. There is a significant loss in constraining
power associated to this cutoff, as seen in fig. 8. Therefore, we provide a
fitting function to improve this.

We find that the concentration-mass relation is one of the key ingredients in
ReACT and have derived a fit for the ratio of its amplitude relative to
$\Lambda$CDM across the MG parameter space. The key results are in equations
(4.1), (4.4), generalised to arbitrary redshift in appendix C, which we
implement as a modification to the $\Lambda$CDM $c(M)$ relation in our
extended version of ReACT. See table 2 for the range of validity of this
expression for the matter power spectrum and weak-lensing calculations, given
different choices of accuracy level and the region of phenomenological MG
parameter space. We are able to achieve 1% accuracy in predicting the matter
power spectrum up to $k=1\,h\,{\rm Mpc}^{-1}$ for the extreme case of
$\mu=1.2$ in each of our bins and 3% accuracy up to $k=3\,h\,{\rm Mpc}^{-1}$.
We achieve larger values of $k_{\rm fail}$ for $\mu=1.1$ and this can be seen
in fig. 10. We show for a Euclid-like survey, how this translates to a 1% (3%)
cut-off of $\ell=900\,(1200)$ for bins with $\mu=1.2$ and $\ell=1300\,(2000)$
for $\mu=1.1$.

Most importantly, the tools developed here allow the data from existing and
future surveys to be used to constrain phenomenological modified gravity.
Performing this analysis robustly will be facilitated by the existing work
demonstrating how baryonic effects and massive neutrions can be included in
the ReACT approach [22].

Our simulations present a useful dataset from which one can recover the
peculiar velocity field and therefore compute the small scale redshift-space-
distortions of galaxy clustering [38]. Cross-correlating the galaxy-clustering
and lensing data can prove to be a powerful tool in the future to probe
gravity on non-linear scales. Such an analysis would provide valuable insight
into the nature of the pixelised $\mu$ approach, which will inform the full
time and scale dependent analysis. In other words, we provide a way to produce
3x2pt forecasts and constraints in a model-independent way.

The time-dependent case studied here and in [14] exhibits rich phenomenology
on non-linear scales, and some theories (such as those encompassed by the
Parameterised-Post-Newtonian Cosmology approach [39, 40, 41]), have no scale
dependence in $\mu$ on non-linear scales. As such a detection or tight
constraints on this purely time-dependent case would be a significant step
forward in terms of testing gravity on cosmological scales. Nonetheless, now
that the time-dependent analysis has been understood and fully validated, this
approach can be extended to a wider analysis that includes the extra
complexity from additionally allowing scale dependent behaviour of $\mu$.

## Acknowledgments

This work was performed partly using the DiRAC Data Intensive service at
Leicester, operated by the University of Leicester IT Services, which forms
part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was
funded by BEIS capital funding via STFC capital grants ST/K000373/1 and
ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of
the National e-Infrastructure. We would like to acknowledge the assistance
given by Research IT, and the use of The HPC Pool funded by the Research
Lifecycle Programme at The University of Manchester. DBT acknowledges support
from the Science and Technology Facilities Council (STFC, grant numbers
ST/P000592/1 and ST/X006344/1). SS was supported formerly by a George Rigg
Scholarship and more recently by the UK Science and Technology Facilities
Council (STFC).

## References

  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physical Reports 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. Odintsov, and V. Oikonomou, Physics Reports 692, 1 (2017).
  * Laureijs  _et al._ [2011] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, B. Garilli, P. Gondoin, L. Guzzo, J. Hoar, H. Hoekstra, R. Holmes, T. Kitching, T. Maciaszek, Y. Mellier, F. Pasian, W. Percival, J. Rhodes, G. Saavedra Criado, M. Sauvage, R. Scaramella, L. Valenziano, S. Warren, R. Bender, F. Castander, A. Cimatti, O. Le Fèvre, H. Kurki-Suonio, M. Levi, P. Lilje, G. Meylan, R. Nichol, K. Pedersen, V. Popa, R. Rebolo Lopez, H. W. Rix, H. Rottgering, W. Zeilinger, F. Grupp, P. Hudelot, R. Massey, M. Meneghetti, L. Miller, S. Paltani, S. Paulin-Henriksson, S. Pires, C. Saxton, T. Schrabback, G. Seidel, J. Walsh, N. Aghanim, L. Amendola, J. Bartlett, C. Baccigalupi, J. P. Beaulieu, K. Benabed, J. G. Cuby, D. Elbaz, P. Fosalba, G. Gavazzi, A. Helmi, I. Hook, M. Irwin, J. P. Kneib, M. Kunz, F. Mannucci, L. Moscardini, C. Tao, R. Teyssier, J. Weller, G. Zamorani, M. R. Zapatero Osorio, O. Boulade, J. J. Foumond, A. Di Giorgio, P. Guttridge, A. James, M. Kemp, J. Martignac, A. Spencer, D. Walton, T. Blümchen, C. Bonoli, F. Bortoletto, C. Cerna, L. Corcione, C. Fabron, K. Jahnke, S. Ligori, F. Madrid, L. Martin, G. Morgante, T. Pamplona, E. Prieto, M. Riva, R. Toledo, M. Trifoglio, F. Zerbi, F. Abdalla, M. Douspis, C. Grenet, S. Borgani, R. Bouwens, F. Courbin, J. M. Delouis, P. Dubath, A. Fontana, M. Frailis, A. Grazian, J. Koppenhöfer, O. Mansutti, M. Melchior, M. Mignoli, J. Mohr, C. Neissner, K. Noddle, M. Poncet, M. Scodeggio, S. Serrano, N. Shane, J. L. Starck, C. Surace, A. Taylor, G. Verdoes-Kleijn, C. Vuerli, O. R. Williams, A. Zacchei, B. Altieri, I. Escudero Sanz, R. Kohley, T. Oosterbroek, P. Astier, D. Bacon, S. Bardelli, C. Baugh, F. Bellagamba, C. Benoist, D. Bianchi, A. Biviano, E. Branchini, C. Carbone, V. Cardone, D. Clements, S. Colombi, C. Conselice, G. Cresci, N. Deacon, J. Dunlop, C. Fedeli, F. Fontanot, P. Franzetti, C. Giocoli, J. Garcia-Bellido, J. Gow, A. Heavens, P. Hewett, C. Heymans, A. Holland, Z. Huang, O. Ilbert, B. Joachimi, E. Jennins, E. Kerins, A. Kiessling, D. Kirk, R. Kotak, O. Krause, O. Lahav, F. van Leeuwen, J. Lesgourgues, M. Lombardi, M. Magliocchetti, K. Maguire, E. Majerotto, R. Maoli, F. Marulli, S. Maurogordato, H. McCracken, R. McLure, A. Melchiorri, A. Merson, M. Moresco, M. Nonino, P. Norberg, J. Peacock, R. Pello, M. Penny, V. Pettorino, C. Di Porto, L. Pozzetti, C. Quercellini, M. Radovich, A. Rassat, N. Roche, S. Ronayette, E. Rossetti, B. Sartoris, P. Schneider, E. Semboloni, S. Serjeant, F. Simpson, C. Skordis, G. Smadja, S. Smartt, P. Spano, S. Spiro, M. Sullivan, A. Tilquin, R. Trotta, L. Verde, Y. Wang, G. Williger, G. Zhao, J. Zoubian, and E. Zucca, arXiv e-prints , arXiv:1110.3193 (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Eifler  _et al._ [2021] T. Eifler, H. Miyatake, E. Krause, C. Heinrich, V. Miranda, C. Hirata, J. Xu, S. Hemmati, M. Simet, P. Capak, A. Choi, O. Doré, C. Doux, X. Fang, R. Hounsell, E. Huff, H.-J. Huang, M. Jarvis, J. Kruk, D. Masters, E. Rozo, D. Scolnic, D. N. Spergel, M. Troxel, A. von der Linden, Y. Wang, D. H. Weinberg, L. Wenzl, and H.-Y. Wu, Monthly Notices of the Royal Astronomical Society 507, 1746 (2021), arXiv:2004.05271 [astro-ph.CO] .
  * LSST Dark Energy Science Collaboration [2012] LSST Dark Energy Science Collaboration, arXiv e-prints , arXiv:1211.0310 (2012), arXiv:1211.0310 [astro-ph.CO] .
  * Square Kilometre Array Cosmology Science Working Group  _et al._ [2018] Square Kilometre Array Cosmology Science Working Group, D. J. Bacon, R. A. Battye, P. Bull, S. Camera, P. G. Ferreira, I. Harrison, D. Parkinson, A. Pourtsidou, M. G. Santos, L. Wolz, F. Abdalla, Y. Akrami, D. Alonso, S. Andrianomena, M. Ballardini, J. L. Bernal, D. Bertacca, C. A. P. Bengaly, A. Bonaldi, C. Bonvin, M. L. Brown, E. Chapman, S. Chen, X. Chen, S. Cunnington, T. M. Davis, C. Dickinson, J. Fonseca, K. Grainge, S. Harper, M. J. Jarvis, R. Maartens, N. Maddox, H. Padmanabhan, J. R. Pritchard, A. Raccanelli, M. Rivi, S. Roychowdhury, M. Sahlen, D. J. Schwarz, T. M. Siewert, M. Viel, F. Villaescusa-Navarro, Y. Xu, D. Yamauchi, and J. Zuntz, arXiv e-prints (2018), arXiv:1811.02743 .
  * Battye and Pearson [2012] R. A. Battye and J. A. Pearson, Journal of Cosmology and Astro-Particle Physics 2012, 019 (2012), arXiv:1203.0398 [hep-th] .
  * Gleyzes  _et al._ [2014] J. Gleyzes, D. Langlois, and F. Vernizzi, International Journal of Modern Physics D 23, 1443010 (2014), arXiv:1411.3712 [hep-th] .
  * DES Collaboration [2018] DES Collaboration, Physical Review D 98, 043526 (2018), arXiv:1708.01530 .
  * Arnold  _et al._ [2019] C. Arnold, M. Leo, and B. Li, Nature Astronomy 3, 945 (2019), arXiv:1907.02977 [astro-ph.CO] .
  * Arnold  _et al._ [2022] C. Arnold, B. Li, B. Giblin, J. Harnois-Déraps, and Y.-C. Cai, Monthly Notices of the Royal Astronomical Society 515, 4161 (2022), arXiv:2109.04984 [astro-ph.CO] .
  * Hassani and Lombriser [2020] F. Hassani and L. Lombriser, arXiv e-prints , arXiv:2003.05927 (2020), arXiv:2003.05927 [astro-ph.CO] .
  * Thomas [2020] D. B. Thomas, Physical Review D 101 (2020), 10.1103/physrevd.101.123517.
  * Srinivasan  _et al._ [2021] S. Srinivasan, D. B. Thomas, F. Pace, and R. Battye, Journal of Cosmology and Astroparticle Physics 2021, 016 (2021), arXiv:2103.05051 [astro-ph.CO] .
  * Cui  _et al._ [2010] W. Cui, P. Zhang, and X. Yang, Physical Review D 81, 103528 (2010), arXiv:1001.5184 [astro-ph.CO] .
  * Zhang  _et al._ [2013] Y. Zhang, P. Zhang, X. Yang, and W. Cui, Physical Review D 87, 023521 (2013), arXiv:1301.3255 [astro-ph.CO] .
  * Bose  _et al._ [2021] B. Bose, B. S. Wright, M. Cataneo, A. Pourtsidou, C. Giocoli, L. Lombriser, I. G. McCarthy, M. Baldi, S. Pfeifer, and Q. Xia., Monthly Notices of the Royal Astronomical Society 508, 2479 (2021), arXiv:2105.12114 [astro-ph.CO] .
  * Bose  _et al._ [2020a] B. Bose, M. Cataneo, T. Tröster, Q. Xia, C. Heymans, and L. Lombriser, Monthly Notices of the Royal Astronomical Society 498, 4650 (2020a), arXiv:2005.12184 [astro-ph.CO] .
  * Bose  _et al._ [2020b] B. Bose, M. Cataneo, T. Tröster, Q. Xia, C. Heymans, and L. Lombriser, arXiv e-prints , arXiv:2005.12184 (2020b), arXiv:2005.12184 [astro-ph.CO] .
  * Mead [2017] A. J. Mead, Monthly Notices of the Royal Astronomical Society 464, 1282 (2017), arXiv:1606.05345 [astro-ph.CO] .
  * Cataneo  _et al._ [2019] M. Cataneo, L. Lombriser, C. Heymans, A. J. Mead, A. Barreira, S. Bose, and B. Li, Monthly Notices of the Royal Astronomical Society 488, 2121 (2019), arXiv:1812.05594 [astro-ph.CO] .
  * Carrilho  _et al._ [2022] P. Carrilho, K. Carrion, B. Bose, A. Pourtsidou, J. C. Hidalgo, L. Lombriser, and M. Baldi, Monthly Notices of the Royal Astronomical Society 512, 3691 (2022), arXiv:2111.13598 [astro-ph.CO] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophysical Journal 462, 563 (1996), astro-ph/9508025 .
  * Klypin  _et al._ [2016] A. Klypin, G. Yepes, S. Gottlöber, F. Prada, and S. Heß, Monthly Notices of the Royal Astronomical Society 457, 4340 (2016), http://oup.prod.sis.lan/mnras/article-pdf/457/4/4340/18515365/stw248.pdf .
  * Planck Collaboration  _et al._ [2018] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J.-F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, D. Herranz, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J.-M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M.-A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A.-S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, and A. Zonca, ArXiv e-prints (2018), arXiv:1807.06209 .
  * Dolag  _et al._ [2004] K. Dolag, M. Bartelmann, F. Perrotta, C. Baccigalupi, L. Moscardini, M. Meneghetti, and G. Tormen, Astronomy and Astrophysics 416, 853 (2004), arXiv:astro-ph/0309771 [astro-ph] .
  * Lombriser  _et al._ [2014] L. Lombriser, K. Koyama, and B. Li, Journal of Cosmology and Astroparticle Physics 2014, 021 (2014), arXiv:1312.1292 [astro-ph.CO] .
  * Shi  _et al._ [2015] D. Shi, B. Li, J. Han, L. Gao, and W. A. Hellwing, Monthly Notices of the Royal Astronomical Society 452, 3179 (2015), arXiv:1503.01109 [astro-ph.CO] .
  * Ruan  _et al._ [2023] C.-Z. Ruan, C. Cuesta-Lazaro, A. Eggemeier, B. Li, C. M. Baugh, C. Arnold, S. Bose, C. Hernández-Aguayo, P. Zarrouk, and C. T. Davies, arXiv e-prints , arXiv:2301.02970 (2023), arXiv:2301.02970 [astro-ph.CO] .
  * Casas  _et al._ [2017] S. Casas, M. Kunz, M. Martinelli, and V. Pettorino, Physics of the Dark Universe 18, 73 (2017), arXiv:1703.01271 [astro-ph.CO] .
  * McDonald  _et al._ [2006] P. McDonald, H. Trac, and C. Contaldi, Monthly Notices of the Royal Astronomical Society 366, 547 (2006), arXiv:astro-ph/0505565 [astro-ph] .
  * Mancini and Bose [2023] A. S. Mancini and B. Bose, “On the degeneracies between baryons, massive neutrinos and f(r) gravity in stage iv cosmic shear analyses,” (2023), arXiv:2305.06350 [astro-ph.CO] .
  * Smail  _et al._ [1994] I. Smail, R. S. Ellis, and M. J. Fitchett, Monthly Notices of the Royal Astronomical Society 270, 245 (1994), arXiv:astro-ph/9402048 [astro-ph] .
  * Hu [1999] W. Hu, Astrophysical Journal, Letters 522, L21 (1999), arXiv:astro-ph/9904153 [astro-ph] .
  * Castro  _et al._ [2005] P. G. Castro, A. F. Heavens, and T. D. Kitching, Physical Review D 72, 023516 (2005), arXiv:astro-ph/0503479 [astro-ph] .
  * Spurio Mancini  _et al._ [2018] A. Spurio Mancini, R. Reischke, V. Pettorino, B. M. Schäfer, and M. Zumalacárregui, Monthly Notices of the Royal Astronomical Society 480, 3725 (2018), arXiv:1801.04251 [astro-ph.CO] .
  * Euclid Collaboration  _et al._ [2020] Euclid Collaboration, A. Blanchard, S. Camera, C. Carbone, V. F. Cardone, S. Casas, S. Clesse, S. Ilić, M. Kilbinger, T. Kitching, M. Kunz, F. Lacasa, E. Linder, E. Majerotto, K. Markovič, M. Martinelli, V. Pettorino, A. Pourtsidou, Z. Sakr, A. G. Sánchez, D. Sapone, I. Tutusaus, S. Yahia-Cherif, V. Yankelevich, S. Andreon, H. Aussel, A. Balaguera-Antolínez, M. Baldi, S. Bardelli, R. Bender, A. Biviano, D. Bonino, A. Boucaud, E. Bozzo, E. Branchini, S. Brau-Nogue, M. Brescia, J. Brinchmann, C. Burigana, R. Cabanac, V. Capobianco, A. Cappi, J. Carretero, C. S. Carvalho, R. Casas, F. J. Castander, M. Castellano, S. Cavuoti, A. Cimatti, R. Cledassou, C. Colodro-Conde, G. Congedo, C. J. Conselice, L. Conversi, Y. Copin, L. Corcione, J. Coupon, H. M. Courtois, M. Cropper, A. Da Silva, S. de la Torre, D. Di Ferdinando, F. Dubath, F. Ducret, C. A. J. Duncan, X. Dupac, S. Dusini, G. Fabbian, M. Fabricius, S. Farrens, P. Fosalba, S. Fotopoulou, N. Fourmanoit, M. Frailis, E. Franceschi, P. Franzetti, M. Fumana, S. Galeotta, W. Gillard, B. Gillis, C. Giocoli, P. Gómez-Alvarez, J. Graciá-Carpio, F. Grupp, L. Guzzo, H. Hoekstra, F. Hormuth, H. Israel, K. Jahnke, E. Keihanen, S. Kermiche, C. C. Kirkpatrick, R. Kohley, B. Kubik, H. Kurki-Suonio, S. Ligori, P. B. Lilje, I. Lloro, D. Maino, E. Maiorano, O. Marggraf, N. Martinet, F. Marulli, R. Massey, E. Medinaceli, S. Mei, Y. Mellier, B. Metcalf, J. J. Metge, G. Meylan, M. Moresco, L. Moscardini, E. Munari, R. C. Nichol, S. Niemi, A. A. Nucita, C. Padilla, S. Paltani, F. Pasian, W. J. Percival, S. Pires, G. Polenta, M. Poncet, L. Pozzetti, G. D. Racca, F. Raison, A. Renzi, J. Rhodes, E. Romelli, M. Roncarelli, E. Rossetti, R. Saglia, P. Schneider, V. Scottez, A. Secroun, G. Sirri, L. Stanco, J. L. Starck, F. Sureau, P. Tallada-Crespí, D. Tavagnacco, A. N. Taylor, M. Tenti, I. Tereno, R. Toledo-Moreo, F. Torradeflot, L. Valenziano, T. Vassallo, G. A. Verdoes Kleijn, M. Viel, Y. Wang, A. Zacchei, J. Zoubian, and E. Zucca, Astronomy and Astrophysics 642, A191 (2020), arXiv:1910.09273 [astro-ph.CO] .
  * Ruan  _et al._ [2022] C.-Z. Ruan, C. Cuesta-Lazaro, A. Eggemeier, C. Hernández-Aguayo, C. M. Baugh, B. Li, and F. Prada, Mon. Not. Roy. Astron. Soc. 514, 440 (2022), arXiv:2110.10033 [astro-ph.CO] .
  * Thomas  _et al._ [2023] D. B. Thomas, T. Clifton, and T. Anton, Journal of Cosmology and Astroparticle Physics 2023, 016 (2023).
  * Clifton and Sanghai [2019] T. Clifton and V. A. Sanghai, Physical Review Letters 122 (2019), 10.1103/physrevlett.122.011301.
  * Sanghai and Clifton [2015] V. A. Sanghai and T. Clifton, Physical Review D 91 (2015), 10.1103/physrevd.91.103532.
  * Troxel  _et al._ [2018] M. Troxel, N. MacCrann, J. Zuntz, T. Eifler, E. Krause, S. Dodelson, D. Gruen, J. Blazek, O. Friedrich, S. Samuroff, J. Prat, L. Secco, C. Davis, A. Ferté, J. DeRose, A. Alarcon, A. Amara, E. Baxter, M. Becker, G. Bernstein, S. Bridle, R. Cawthon, C. Chang, A. Choi, J. D. Vicente, A. Drlica-Wagner, J. Elvin-Poole, J. Frieman, M. Gatti, W. Hartley, K. Honscheid, B. Hoyle, E. Huff, D. Huterer, B. Jain, M. Jarvis, T. Kacprzak, D. Kirk, N. Kokron, C. Krawiec, O. Lahav, A. Liddle, J. Peacock, M. Rau, A. Refregier, R. Rollins, E. Rozo, E. Rykoff, C. Sánchez, I. Sevilla-Noarbe, E. Sheldon, A. Stebbins, T. Varga, P. Vielzeuf, M. Wang, R. Wechsler, B. Yanny, T. Abbott, F. Abdalla, S. Allam, J. Annis, K. Bechtol, A. Benoit-Lévy, E. Bertin, D. Brooks, E. Buckley-Geer, D. Burke, A. C. Rosell, M. C. Kind, J. Carretero, F. Castander, M. Crocce, C. Cunha, C. D’Andrea, L. da Costa, D. DePoy, S. Desai, H. Diehl, J. Dietrich, P. Doel, E. Fernandez, B. Flaugher, P. Fosalba, J. García-Bellido, E. Gaztanaga, D. Gerdes, T. Giannantonio, D. Goldstein, R. Gruendl, J. Gschwend, G. Gutierrez, D. James, T. Jeltema, M. Johnson, M. Johnson, S. Kent, K. Kuehn, S. Kuhlmann, N. Kuropatkin, T. Li, M. Lima, H. Lin, M. Maia, M. March, J. Marshall, P. Martini, P. Melchior, F. Menanteau, R. Miquel, J. Mohr, E. Neilsen, R. Nichol, B. Nord, D. Petravick, A. Plazas, A. Romer, A. Roodman, M. Sako, E. Sanchez, V. Scarpine, R. Schindler, M. Schubnell, M. Smith, R. Smith, M. Soares-Santos, F. Sobreira, E. Suchyta, M. Swanson, G. Tarle, D. Thomas, D. Tucker, V. Vikram, A. Walker, J. Weller, and Y. Zhang, Physical Review D 98 (2018), 10.1103/physrevd.98.043528.

##  Appendix A Spherical Collapse in ReACT

In this appendix, we describe the physics behind the ReACT approach and
provide some qualitative justification for the shape of the fitting function
in our extended implementation. The ReACT method augments the halo model with
the spherical collapse formalism. The modifications to gravity enter via the
Poisson equation (2.1). Halo formation is modelled by following the evolution
of a spherical top-hat overdensity of radius $R_{\rm TH}$ in terms of the
gravitational potential and $\mu(z)$. The equation solved within the code is
given by [19]

$y^{\prime\prime}+\frac{H^{\prime}}{H}y^{\prime}-\left(1+\frac{H^{\prime}}{H}\right)y+\left[\frac{H_{0}}{H}\right]^{2}\frac{\Omega_{\rm
m,0}}{2a^{3}}\mu(z)\delta\left(\frac{a}{a_{i}}+y\right)=0$ (A.1)

where y = $R_{\rm TH}/R_{i}-a/a_{i}$, $a_{i}$ is the initial scale factor,
$R_{i}$ is the initial radius of the overdensity and the prime denotes a
derivative w.r.t $\ln(a)$. The initial conditions are set deep into matter
domination when an Einstein-de-Sitter Universe is an accurate description. The
initial overdensity $\delta_{i}$ is iteratively adjusted until the collapse
condition corresponding to $R_{\rm TH}=0$ at characterised by the scale factor
at the time of collapse $(a_{\rm coll})$ is reached. One can then apply the
virial theorem (see appendix A in [21] for details) to determine the viral
overdensity parameter given by

$\Delta_{\rm vir}=[1+\delta(a_{\rm vir})]\left(\frac{a_{\rm coll}}{a_{\rm
vir}}\right)^{3}\,.$ (A.2)

The results of the spherical collapse calculation are used to compute a
modified Sheth-Tormen mass function given by

$n_{\rm vir}=\frac{dn}{d\ln M_{\rm vir}}=\frac{\bar{\rho}_{\rm m,0}}{M_{\rm
vir}}\nu f(\nu)\frac{d\ln\nu}{d\ln M_{\rm vir}}\,,$ (A.3)

where $\nu=\delta_{\rm c}/\sigma$ is the peak height with $\sigma$ is the
variance in the amplitude of the density fluctuations within a top-hat of a
fixed radius (typically set to $8\,h^{-1}{\rm Mpc}$). Note that the inclusion
of the spherical collapse calculation leads to a different peak height
entering the halo model. The concentration-mass relation measured in
simulations have a weak negative power law, meaning the largest halos have the
smallest $c$. Since larger halos typically form at late times, the spherical
collapse results shown in fig. 9 indicate that one would expect that variation
in the fitting function for the amplitude $\hat{c}=c/c_{\rm\Lambda CDM}$ of
the concentration-mass relation in ReACT will be more dramatic at small
redshifts.

Effectively, the change in the redshift at which halo-collapse and
virialisation occur is propagated into the mass function. The non-linear
matter power spectrum may be computed using the halo model via the sum of the
1-halo and 2-halo terms given by

$\displaystyle P_{1h}$ $\displaystyle=$ $\displaystyle\int{\rm d}\ln M_{\rm
vir}n_{\rm vir}\left(\frac{M_{\rm vir}}{\bar{\rho}_{\rm
m,0}}\right)^{2}\left|u(k,M_{\rm vir})\right|^{2}\,,$ (A.4) $\displaystyle
P_{2h}$ $\displaystyle=$ $\displaystyle\left(\int{\rm d}\ln M_{\rm vir}n_{\rm
vir}\frac{M_{\rm vir}}{\bar{\rho}_{\rm m,0}}u(k,M_{\rm vir})b_{\rm L}(M_{\rm
vir})\right)^{2}P_{\rm L}(k)\,,$ (A.5)

where $u(k,m)$ is the Fourier transform of the density profile, $b_{\rm L}$ is
the linear bias parameter and $P_{\rm L}(k)$ is the linear power spectrum.

ReACT also implements second order perturbation-theory in order to get better
accuracy for the specific behaviour seen in $f(R)$ gravity and DGP gravity.
However, in our formalism, all modifications to gravity are captured by the
parameter $\mu(z)$. Therefore, we switch off the flags that implement these
routines within the code (see appendix B in [14]).

Figure 9: We show the solution to the differential equation (A.1) for
different redshift bins for a constant value of $\mu$. Clearly, the spherical
collapse parameter is more dramatically affected when $\mu$ is modified at
late times.

##  Appendix B Sensitivity of analysis to simulation parameters

In our analysis in this work as well as in [14], we ran our $N$-body
simulations with Planck-2018 [25] (TT+TE+low P +BAO) cosmological parameters
and fixed random seed for the initial conditions. Since we compute the ratio
of matter power spectra with the identical initial conditions, our results are
realisation-independent [31]. We verified this by computing the $R(z,k)$ for
two of our simulations after averaging over three realisations and found
identical results. We also investigated the effects of varying the standard
cosmological parameters. We obtain the same $R(z,k)$ at both $z=0$ and $z=0.5$
when we vary cosmological parameters, to the DES Y1 cosmological parameters
[42] and the joint Planck-DES Y1 analysis [25]. These results justify that the
results in the main text are not dependent on the realisations that were
examined, or on the $\Lambda$ CDM parameters. This ensures that the fitting
function has a wide range of validity and is suitable for future data
analyses.

We make use of the $k_{\rm fail}$ statistic in order to derive the fitting
parameters in our extended implementation of ReACT (see sec. 4.2). We derive
our fitting parameters by imposing 1% accuracy threshold in computing $k_{\rm
fail}$. This might be an overly conservative choice considering that recent
release of ReACT that incorporated baryonic modelling [17] reported an
accuracy of 3%. In fig. 10, we show the value of $k_{\rm fail}$ we achieve for
1%, 2% and 3% thresholds, respectively for across the range of our simulation
parameters $(\mu,D(z))$. We see that our choice of using $k_{\rm fail}$ means
that we are most sensitive to variation in $\mu$ rather than $D(z)$. This
motivates our choice of presenting our weak-lensing results for the extreme
values of $\mu$ in our simulation suite.

Figure 10: We present the $k_{\rm fail}$ parameter with different accuracy
thresholds for our simulations across the our parameter-space of ($\mu,D(z)$).
This is particularly relevant for the lensing prediction since we would cut-
off $P(k)$ at $k_{\rm fail}$ in the computation of the convergence power
spectrum. In particular, our choice of $1\%$ in the main text may be
conservative, given that baryonic effects may limit the accuracy of ReACT to
$3\%$. We also note that the variation in $k_{\rm fail}$ is substantial
between $\mu=1.1$ and $\mu=1.2$, whereas the difference is much less
significant between the two cases $D=1.26$ and $D=1.45$. These observations
motivate out choice of presenting the corresponding weak-lensing $\ell_{\rm
fail}$ as a function of $\mu$ for ‘conservative’ and ‘realistic’ cases,
respectively.

##  Appendix C Concentration Fitting function

In the main text, we derive the fitting function for the amplitude of the
concentration-mass relation as a function of $\bar{z}$ to be an exponential at
$\bar{z}>1$ and a cubic polynomial at $\bar{z}<1$. We noted that the fitting
parameters vary linearly as a function of $D$ and quadratically as a function
of $\mu$ (see fig. 11. We thus derived the general form of the fitting
parameters at $z_{\rm pk}=0$. We now give the explicit functional forms we use
to construct the fitting parameters at $z_{\rm pk}=0$.

We derive the fitting parameters $\tilde{C}_{i}$ as a function of $\mu$ at
$D=1.26$ to be $C_{i}=\gamma_{i}(\mu-1)+\delta_{i}(\mu-1)^{2}$. We show in
table 3 the values of these parameters for the extreme values of $D(z)$. Note
that these have been derived for $\mu>1$, but can be used to to derive the fit
for arbitrary $\mu$. We will return to this point shortly. We derive the
$\hat{C}_{i}$ at arbitrary values of $D(z)$ by linearly interpolating between
$\tilde{C}_{i}$ and $\bar{C}_{i}$ see eq. (4.4) in main text.

$D(z)$ | $C_{i}$ | $\gamma_{i}$ | $\delta_{i}$
---|---|---|---
1.26 | $\tilde{C}_{1}$ | 1.9 | 1.3
$\tilde{C}_{2}$ | -0.6 | -122.9
$\tilde{C}_{3}$ | 172.7 | 653.0
$\tilde{C}_{4}$ | -450.6 | -1579.0
$\tilde{C}_{5}$ | 328.5 | 1145.6
1.45 | $\bar{C}_{1}$ | 4.9 | 14.0
$\bar{C}_{2}$ | -154.7 | 634.0
$\bar{C}_{3}$ | 904.5 | -3610.1
$\bar{C}_{4}$ | -1438.7 | 5700.0
$\bar{C}_{5}$ | 669.6 | -2648.0
Table 3: The fitting parameters $\tilde{C}_{i}$ and $\bar{C}_{i}$ as a
function of $\mu$. We linearly interpolate between these values for the
general case of arbitrary $D(z)$.

Figure 11: The variation of the fitting parameters as a function of $\mu$ at
$z=0$.

While this fit is derived for $\mu>1$, we find that there is a symmetry in the
shape of the matter power spectrum across $\mu=1$ (see figs. 3 and 4 and the
relevant discussion in the main text). Therefore, $A/A_{\rm\Lambda CDM}$ for
$\mu<1$ can be computed by taking the reciprocal of the fit at the reflected
value at $\mu>1$. We show an example of this behaviour in fig. 12.

Figure 12: We show the symmetry in the shape of the matter power spectrum
across $\mu=1$ is also seen in the fitting function. In the above figure, the
data points are derived from the simulations for $\mu=1.1$ (circles) and
$\mu=0.9$, (crosses) respectively but we only derive the fit for $\mu=1.1$. We
then derive the fit to the $\mu=0.9$ datapoints by taking the reciprocal of
the fit the $\mu=1.1$ points.

We find that the general functional form of the fit is preserved at all
redshifts. However, we find that the ratio of the amplitudes of the
concentration parameter relative to $\Lambda$CDM $A/A_{\Lambda CDM}$ is
inversely proportional to redshift. As a result, at $z_{\rm pk}\geq 1$ we find
that vanilla implementation of ReACT (see sec. 4) achieves a $k_{\rm fail}\geq
3\,h\,{\rm Mpc}^{-1}$ at an accuracy threshold of 1%. Therefore, we derive the
fitting function at arbitrary redshifts in the range $0\leq z_{\rm pk}\leq 1$
by linearly interpolating between the redshift at which we know the fitting
parameters, i.e., $z_{\rm pk}=\\{0,0.5,1\\}$ 111111We will be publicly
releasing this code with our extended ReACT implementation shortly. We show an
example of this for $D=1.26$ in fig. 13.

Figure 13: We show the fitting function at $z=0$ and $z_{\rm pk}=0.5$, where
the circular and cross data points represent the measurement of the fit at
$z_{\rm mp}$ corresponding to the simulations at $z_{\rm pk}=0$ and $z_{\rm
pk}=0.5$, respectively. We linearly interpolate between the these two curves
to obtain the fitting parameters at $z_{\rm pk}=0.25$, and indeed at other
redshifts to produce the weak-lensing convergence spectra in the main text
(see fig. 7). In other words, we do not use any simulation data to produce the
curve at $z_{\rm pk}=0.25$, and still obtain $\ell_{\rm fail}>800$
consistently across our simulation set as well as agreement with the
simulations (blue data points).