# Stage-IV Cosmic Shear with Modified Gravity and Model-independent Screening

M. Tsedrik,11footnotetext: Corresponding author. B. Bose  P. Carrilho  A.
Pourtsidou  S. Pamuk  S. Casas  J. Lesgourgues

###### Abstract

We forecast constraints on minimal model-independent parametrisations of
several Modified Gravity theories using mock Stage-IV cosmic shear data. We
include nonlinear effects and screening, which ensures recovery of General
Relativity on small scales. We introduce a power spectrum emulator to
accelerate our analysis and evaluate the robustness of the growth index
parametrisation with respect to two cosmologies: $\Lambda$CDM and the normal
branch of the DGP model. We forecast the uncertainties on the growth index
$\gamma$ to be of the order $\sim 10\%$. We find that our halo-model based
screening approach demonstrates excellent performance, meeting the precision
requirements of Stage-IV surveys. However, neglecting the screening transition
results in biased predictions for cosmological parameters. We find that the
screening transition shows significant degeneracy with baryonic feedback,
requiring a much better understanding of baryonic physics for its detection.
Massive neutrinos effects are less prominent and challenging to detect solely
with cosmic shear data.

##  1 Introduction

Stage-IV cosmological surveys, such as Euclid222http://euclid-ec.org [1] and
the Vera C. Rubin Observatory’s Legacy Survey of Space and Time
(LSST)333https://www.lsst.org/ [2], will provide us with high precision data,
enabling us to constrain deviations from General Relativity (GR). GR serves as
the foundational theory of gravity in the standard cosmological model,
$\Lambda$CDM. Enhancing the robustness of our modelling has become imperative,
as statistical uncertainties will no longer be the limiting obstacle; instead,
Stage-IV data analyses will be limited by the accuracy of the theoretical
modelling of the observables.

Cosmic shear, the weak gravitational lensing effect that systematically
distorts the shapes of galaxy images, can be directly connected to the spatial
distribution of all gravitating matter along the line-of-sight (for
comprehensive reviews, see Refs. [3, 4, 5]). For Stage-IV surveys, the small
scales of the cosmic shear measurement exhibit the highest signal-to-noise
ratio [6]. Their inclusion will be of the utmost importance to distinguish
between competing gravity models with distinct nonlinear behaviors. Nontrivial
effects that should be taken into account include nonlinear structure
formation [7], baryonic effects on the matter distribution [8], intrinsic
alignments [9] and the impact of massive neutrinos on structure formation [10,
11]. For beyond-$\Lambda$CDM cosmologies, which include extensions of GR, the
screening mechanism [12, 13, 14, 15] is an additional important nonlinear
effect. Screening is a key phenomenological aspect of all viable Modified
Gravity (MG) theories, allowing recovery of GR on small scales. At small
scales, notably within our Solar System, GR exhibits highly accurate
observational consistency (see Ref. [16] for a review).

There is a plethora of non-standard cosmologies (for reviews, see e.g. Refs.
[17, 18]), and a few of them have been tested with real photometric data (see,
for example, Refs. [19, 20, 21, 22]) and simulations [23, 24, 25, 26], or have
been subject to forecasts for Stage-IV weak lensing surveys [27, 28, 29]. An
alternative approach to studying a broader class of MG models involves
Horndeski theories [30], the corresponding $\alpha$-parametrisation [31, 32,
33] or direct parametrisation of the linear relation between matter and the
gravitational potential (see, for example, Ref. [34]). Typically, in such
parameterised analyses, screening has either been neglected [35, 36] or
accounted for through conservative scale cuts [37, 38, 22, 39].

To our knowledge, attempts to take the screening effects into account have
been made in Refs. [40, 41, 42, 43]. Therein, screening was modeled as an
exponential cut-off in the $\alpha$-parameters with a characteristic screening
scale at which the screening mechanism becomes effective. This scale was
either varied as a free parameter or fixed to the value motivated by
simulations [44]. However, the performance of this approximation across
various MG scenarios remains to be tested, in order to understand if the
corresponding accuracy suffices for Stage-IV surveys.

In Ref. [45], we introduced an approximate phenomenological model based on the
error function to model the screening effects within the halo model reaction
framework [46]. Various extended scenarios were shown to be in good agreement
(within a few percent) with the exact solution for the screened regime. In the
same paper, we described the minimal extension model, where the background is
assumed to be flat $\Lambda$CDM, and the deviation of linear growth from GR is
parameterised by the growth index $\gamma$ [47, 48, 49]. The growth index is
defined via the growth rate, $f={\rm d}\ln D/{\rm d}\ln a$, as an exponent of
the time-dependent matter density relative to critical:

$f(z)=\Omega_{\rm m}^{\gamma}(z)\,.$ (1.1)

This phenomenological parameterisation was developed with galaxy clustering
probes in mind and has been used to look for deviations from standard
cosmology [50, 51, 52]. In Ref. [53] it was demonstrated that $\gamma$ can be
constrained with weak lensing probes as well. Constraining the growth index is
one of the main objectives of Stage-IV surveys [54, 55, 6]. We should note
that this parametrisation assumes that the linear growth factor and rate are
scale-independent which is not the case for some modified gravity models, such
as the well studied $f(R)$ model of Ref. [56]. To account for such scale
dependencies one would need to extend this parametrisation as done in Ref.
[57]. We will not consider this here and restrict ourselves to the scale-
independent case.

In this paper, we investigate whether we could provide constraints on the
minimal model-independent parametrisation of MG, including nonlinear scales
with forthcoming Stage-IV cosmic shear data. We present an emulator, trained
with cosmopower444github.com/alessiospuriomancini/cosmopower on power spectra
generated using the halo model reaction code ReACT555github.com/nebblu/ACTio-
ReACTio, which will be used in our analyses. We have made a public repository
for this and related emulators 666github.com/nebblu/ReACT-emus/emulators.

In Section 2, we introduce the modelling of cosmic shear and relevant
nonlinear effects. In Section 3, we present our analysis setup. In Section 4,
we study how well the nonlinear extension of the growth index performs on mock
data in GR and an example MG model. Then, we investigate the importance of
including the screening scale, the degeneracy between the screening scale and
baryonic feedback, and the effects from massive neutrinos. Finally, we
summarise our results and outline the next necessary steps, paving the way to
an optimal model-independent parametrisation on all scales.

##  2 Cosmic Shear Modelling

In this section, we describe the modelling of our observable, the cosmic shear
power spectrum, detailing its main components, the MG scenarios under
consideration, and their signatures on the weak lensing power spectrum.

###  2.1 Cosmic Shear

For a generic theory of gravity, the cosmic shear power spectrum can be
derived from the lensing potential, also known as the Weyl potential. The Weyl
potential is intricately linked to the underlying matter density fluctuations
via a corresponding Poisson equation. Therefore, we can obtain the cosmic
shear power spectrum by integrating the nonlinear matter power spectrum along
the line-of-sight under the Limber approximation [41]:

$C^{\gamma\gamma}_{ij}(\ell)=\int_{z_{\rm min}=0}^{z_{\rm
max}}\mathrm{d}z\frac{W^{\gamma}_{i}(z)W^{\gamma}_{j}(z)}{H(z)\chi^{2}(z)}\Sigma^{2}[k_{\ell}(z),z]P_{\rm
NL}[k_{\ell}(z),z]\,,$ (2.1)

where we set the speed of light to $c=1$, $k_{\ell}=(\ell+1/2)/\chi(z)$,
$\chi(z)$ signifies the radial comoving distance from an observer at $z=0$ to
an object at redshift $z$, $H(a)$ represents the Hubble function,
$\Sigma(k,z)$ denotes the modification to the Poisson equation for the Weyl
potential in MG theories, and $z_{\rm max}$ stands for the maximum redshift of
the source distribution in a survey. The weighting functions $W^{\gamma}_{i}$
are defined as:

$W^{\gamma}_{i}(z)=\frac{3}{2}H_{0}^{2}\Omega_{\rm m}\chi(z)(1+z)\int_{z_{\rm
min}=0}^{z_{\rm
max}}\mathrm{d}z^{\prime}n_{i}(z^{\prime})\left[1-\frac{\chi(z)}{\chi(z^{\prime})}\right]\,,$
(2.2)

where $\Omega_{\rm m}=\Omega_{\rm m}(z=0)$ denotes the total matter density
fraction today, $H_{0}$ represents the Hubble constant in $\rm{Mpc}^{-1}$, and
$n_{i}$ is the redshift distribution for bin $i$. The observed lensing signal
$C^{\rm WL}_{ij}(\ell)$ is contaminated by the intrinsic alignment
contribution (IA, see Ref. [58] for a recent review). IA represent
correlations in the orientation between galaxies which are not caused by
lensing but rather by the same gravitational field in which these galaxies
evolve. We model this using the redshift-dependent nonlinear alignment (zNLA)
model [59, 60]. Additionally, we account for the shape-noise contribution due
to the intrinsic ellipticity field (unlensed). These effects introduce the
following additional contribution to the signal:

$C^{\rm
WL}_{ij}(\ell)=C^{\gamma\gamma}_{ij}(\ell)+C^{I\gamma}_{ij}(\ell)+C^{II}_{ij}(\ell)+\frac{\sigma_{\epsilon}^{2}}{\bar{n}}\delta^{\rm
K}_{ij}\,,$ (2.3)

where $\delta^{\rm K}_{ij}$ is the Kronecker’s delta, $\sigma_{\epsilon}$
stands for the variance of the intrinsic ellipticity distribution,
$\bar{n}=n/N_{\rm bins}$ with $n$ being number of galaxies per radians
squared, and the superscript $I$ denotes the contribution from the IA with the
following kernel:

$W^{I}_{i}(k,z)=-A_{\rm IA}C_{\rm IA}(1+z)^{\eta_{\rm IA}}\frac{\Omega_{\rm
m}}{D(k,z)/D(k,0)}n_{i}(z)H(z)\,,$ (2.4)

where we set $C_{\rm IA}=0.0134$ (its conventional value, as it is degenerate
with $A_{\rm IA}$), while $A_{\rm IA}$ and $\eta_{\rm IA}$ are additional
nuisance parameters in our modelling. $D(k,z)$ denotes the growth
factor777This is found by solving the linearised growth equation in the
cosmological scenario of interest, which we compute using MGrowth:
github.com/MariaTsedrik/MGrowth. – the general definition for MG theories
includes a scale dependence, but in this work we only study models with scale-
independent linear growth. In general, IA occurs on small astrophysical
scales, where we assume that all modifications of gravity are screened [33].
The zNLA model allows for amplitude and redshift dependence of the IA signal
from the gravity model through $D(k,z)$ in Eq. 2.4 and $P_{\rm NL}$ in Eq.
2.1. Unlike other IA models where tidal physics is taken into account and
perturbatively re-derived for MG scenarios, this does not affect NLA [22].
While there is still debate on whether or not zNLA captures all IA effects
sufficiently well for Stage-IV surveys [61], we leave a detailed exploration
of more complex models for future work.

The main component from Eq. 2.1 that we model and modify in our analysis is
the nonlinear power spectrum:

$P_{\rm NL}(k,z)=P^{\mathrm{MG}+\nu}_{\rm NL}(k,z)\times B^{\rm
baryons}(k,z)=B^{\mathrm{MG}+\nu}\times B^{\rm baryons}\times P^{\Lambda\rm
CDM}_{\rm NL}\,,$ (2.5)

where the nonlinear $\Lambda$CDM power spectrum, $P^{\Lambda\rm CDM}_{\rm
NL}$, is combined with the emulated boosts due to effects from baryonic
feedback, $B^{\rm baryons}$, and MG with massive neutrinos,
$B^{\mathrm{MG}+\nu}$. This prescription is described in detail in the
following subsections.

###  2.2 Nonlinear Power Spectrum Modelling

In the halo model reaction framework888 The series of papers “On the road to
percent accuracy” covers the developments of the framework including the basic
foundation [46], an emulator for the pseudo power spectrum [62], inclusion of
massive neutrinos [63], the first forecasts for an LSST-like survey and the
ReACT code release [64], a comparison with simulations and simulation-based
emulators including baryonic contribution in [65], and interactions between
dark matter and dark energy [66]. [46], based on the halo-model approach (see
Ref. [67] for a recent review), we can compute the nonlinear power spectrum in
a modified theory of gravity including massive neutrinos as:

$P^{\mathrm{MG}+\nu}_{\rm NL}(k,z)=\mathcal{R}(k,z)P^{\rm pseudo}_{\rm
NL}(k,z)\,,$ (2.6)

where $P^{\rm pseudo}_{\rm NL}(k,z)$ is called the pseudo power spectrum. This
is defined as a nonlinear power spectrum, evolved in a $\Lambda$CDM universe
with adjusted initial conditions in order to match the linear clustering of
the MG model of interest at the target redshifts, i.e., $P^{\rm pseudo}_{\rm
L}(k,z_{j})=P^{\rm MG}_{\rm L}(k,z_{j})$. In the context of N-body simulations
it can be understood as follows: initial conditions are generated from a power
spectrum computed as $P^{\rm pseudo}_{\rm L}(k,z_{\rm ini})=D^{2}_{\Lambda\rm
CDM}(z_{\rm ini})/D^{2}_{\Lambda\rm CDM}(z_{j})~{}P_{\rm L}^{\rm MG}(k,z_{j})$
while the evolution up to $z_{j}$ is done with standard gravity and
$\Lambda$CDM expansion. In the context of HMcode it can be understood as using
$D^{2}_{\Lambda\rm CDM}(z=0)/D^{2}_{\Lambda\rm CDM}(z_{j})~{}P^{\rm
MG}(k,z_{j})$ as input linear power spectrum at redshift $z=0$. Since in our
case the re-scaling of the initial conditions is scale-independent (as we are
dealing with MG theories with scale-independent linear growth), one can also
compute a pseudo power spectrum at the target redshift simply by re-scaling
the primordial amplitude as $A^{\rm pseudo}_{s}=D^{2}_{\rm
MG}(z_{j})/D^{2}_{\Lambda\rm CDM}(z_{j})A_{s}$. Instead of utilising a
$\Lambda$CDM nonlinear power spectrum with the same cosmological parameters
and modelling the reaction as an expected ratio $P^{\mathrm{MG}+\nu}_{\rm
NL}/P^{\Lambda\mathrm{CDM}}_{\rm NL}$ with halo model, we opt for the pseudo
cosmology, which implies the halo mass functions of both MG and pseudo
cosmologies are similar. As a consequence, the transition between 2- and
1-halo terms becomes smoother. For example, in Fig. 5 and 7 of Ref. [46], when
comparing $P^{\mathrm{MG}}_{\rm NL}/P^{\Lambda\mathrm{CDM}}_{\rm NL}$ and
$P^{\mathrm{MG}}_{\rm NL}/P^{\rm pseudo}_{\rm NL}$, the latter clearly appears
as a smoother function.

We compute the pseudo power spectrum using HMCode [68, 69, 70, 71], in
particular its 2020-version. Alternatively, the pseudo power spectrum can be
computed with HaloFit [72] or for MG theories with scale-independent linear
growth, $\Lambda$CDM-based emulators such as EuclidEmulator2 [73, 74] or
BACCOemu [75] can be used by adjusting the spectrum amplitude parameter,
$A_{s}$ or $\sigma_{8}$, to match the modified cosmology.

The reaction function $\mathcal{R}(k,z)$ contains all nonlinear corrections to
the pseudo spectrum coming from modifications of gravity and massive
neutrinos. We refer the reader to Refs. [46, 64] for more details on how to
model the reaction function. In a nutshell, similar to the halo-model
approach, there are three distinct regimes in scale: linear (2-halo term),
quasi-nonlinear (1-loop corrections and smoothing), and nonlinear (1-halo
term). For each we need to follow a certain prescription within our framework.
These three regimes are demonstrated clearly in the Poisson equations that
connect the gravitational potential, $\Phi$, to the matter density
fluctuations, $\delta$:

$\displaystyle-\left(\frac{k}{aH(a)}\right)^{2}\Phi_{\rm
QNL}(\mbox{\boldmath$k$},a)=$ $\displaystyle\frac{3\Omega_{\rm
m}(a)}{2}\mu_{\rm L}(k,a)\,\delta_{\rm
QNL}(\mbox{\boldmath$k$},a)+S(\mbox{\boldmath$k$},a)\,,$ (2.7)
$\displaystyle-\left(\frac{k}{aH(a)}\right)^{2}\Phi_{\rm
NL}(\mbox{\boldmath$k$},a)=$ $\displaystyle\frac{3\Omega_{\rm
m}(a)}{2}\mu_{\rm NL}(k,a)\,\delta_{\rm NL}(\mbox{\boldmath$k$},a)\,,$ (2.8)

where $\Omega_{\rm m}(a)=\Omega_{\rm{m}}a^{-3}H_{0}^{2}/H^{2}(a)$. The
subscripts QNL and NL denote “quasi-nonlinear” and “nonlinear”, respectively.

In the linear regime, the enhancement or suppression of structure formation
relative to $\Lambda$CDM is controlled by $\mu_{\rm L}(k,a)$. In the quasi-
nonlinear regime, the modifications at second and third orders are captured by
a source term $S(\mbox{\boldmath$k$},a)$ (see Eq. 2.8 in Ref. [76]). Within
the reaction function, two parameters guarantee a smooth transition between
the 2- and 1-halo terms in the quasi-nonlinear regime: $\mathcal{E}$ and
$k_{\star}$. The first parameter corresponds to the 1-halo power spectrum
ratio in the modified and pseudo cosmologies at very large scales, tuning the
similarities in their halo mass functions. The second parameter controls the
rate of the transition, and for that parameter the one-loop corrections given
by $S(\mbox{\boldmath$k$},a)$ are essential.

In the nonlinear regime, the modification to gravity is governed by $\mu_{\rm
NL}(k,a)$, which should recover GR at very nonlinear scales $(k\gg
10\,h/\mathrm{Mpc})$ in the screened MG theories $\mu_{\rm NL}\rightarrow 1$
or $\mu_{\rm NL}=\mu_{\rm L}$ in the unscreened MG theories. This function is
then considered when solving the gravitational collapse equation for the top-
hat radius. From the solution of the gravitational collapse equation and
Virial theorem, we derive the density at the collapse, virial scale factor and
virial density, which are then used for computing the halo-mass function [77,
78] and halo density profile [79] – two essential properties of the 1-halo
term power spectrum.

All in all, to fully describe $P^{\mathrm{MG}+\nu}$ we require the following
information: the total neutrino mass $M_{\nu}$ (or equivalently the neutrino
density parameter, related to the mass by $M_{\nu}\approx
93.14~{}\Omega_{\nu}h^{2}$ eV), and 5 functions – the expansion history
$H(a)$, the modification of gravity on linear scales $\mu_{\rm L}(k,a)$, two
functions $\gamma_{2}$ and $\gamma_{3}$ included in the source term
$S(\mbox{\boldmath$k$},a)$ that modifies quasi-nonlinear scales, and the
modification of gravity in the nonlinear regime $\mu_{\rm NL}(k,a)$. In the
subsequent subsection, we discuss these functions for the MG theories of
interest.

###  2.3 Modified Gravity Scenarios

In our previous work [45], we studied Horndeski theories employing the
Effective Field Theory of Dark Energy (EFTofDE). Horndeski theories [80]
encompass the most general class of scalar-tensor theories of gravity in 4
dimensions that are Lorentz-covariant, ghost-free and have second-order
equations of motion. However, in this work, we specifically focus on a popular
example of scalar-tensor theories with Vainshtein screening [81] – the normal
branch of DGP theories (nDGP, [82]). We aim to recover its main features with
the phenomenological growth index parameterisation [47, 48, 49] extended to
the nonlinear regime (see the “minimal parametrisation" in Ref. [45]).

Considering linear scalar perturbations on a Friedmann-Robertson-Walker
metric, one can write the line element in the conformal Newtonian gauge and a
spatially flat background as:

$\mathrm{d}s^{2}=-(1+2\Phi)\mathrm{d}t^{2}+a^{2}(1-2\Psi)\mathrm{d}\mbox{\boldmath$x$}^{2}\,,$
(2.9)

with the Newtonian potential $\Phi$ and the intrinsic spatial curvature
potential $\Psi$. We can then modify the linear Poisson equations neglecting
the anisotropic shear as

$\displaystyle-\left(\frac{k}{aH(a)}\right)^{2}\Phi_{\rm
L}(\mbox{\boldmath$k$},a)=$ $\displaystyle\frac{3\Omega_{\rm m}(a)}{2}\mu_{\rm
L}(k,a)\,\delta_{\rm L}(\mbox{\boldmath$k$},a)\,,$ (2.10)
$\displaystyle\Psi_{\rm L}=$ $\displaystyle\eta_{\rm L}(k,a)\Phi_{\rm L}\,,$
(2.11)

where $\eta_{\rm L}$ is the slip parameter, and which simplifies to the
standard GR case for $\mu_{\rm L}=\eta_{\rm L}=1$. We now obtain a Poisson
equation for the Weyl potential, $\Psi_{\rm W}=(\Phi_{\rm L}+\Psi_{\rm L})/2$
as

$-\left(\frac{k}{aH(a)}\right)^{2}\Psi_{\rm
W}(\mbox{\boldmath$k$},a)=\frac{3\Omega_{\rm m}(a)}{2}\Sigma_{\rm
L}(k,a)\,\delta_{\rm L}(\mbox{\boldmath$k$},a)\,,$ (2.12)

with $\Sigma_{\rm L}(k,a)=\mu_{\rm L}(1+\eta_{\rm L})/2$, which equals unity
in the GR case. Generally, most MG theories exhibit a strong preference to a
no-slip condition so that $\Phi_{\rm L}\approx\Psi_{\rm L}$ (see Fig. 2-4 in
Ref. [83]), hence $\Sigma_{\rm L}=\mu_{\rm L}$. We can further extend the
equality of the potentials into the nonlinear regime, which implies
$\Sigma(k,z)=\mu_{\rm NL}(k,z)$ in Eq. 2.1: forcing not only $P_{\rm NL}$ but
also $\Sigma$ to converge to a GR-limit on small scales. However, for our
theories of interest we have $\Sigma_{\rm L}=1$, i.e., lensing is not
affected: for nDGP this is derived in Refs. [84, 85] and for the
$\gamma$-parametrisation we set it by hand. Therefore, in this work, in Eq.
2.1 we set $\Sigma(k,z)=1$ and modify $P_{\rm NL}$ only.

Figure 1: Upper panel: left – linear modification to the Poisson equation for
nDGP gravity with strong modification of $\Omega_{rc}=0.25$ (red solid line),
GR (black line) and different growth index values specified in the colorbar;
right – nonlinear modification to the Poisson equation for the same nDGP
scenario, and growth index with $\gamma=0.4$ and various values of the
screening scale $q_{1}$ at $z=1.25$ as a function of the normalised halo
radius $y_{\rm h}$, defined in Eq. 2.18. Lower panel: ratio between the growth
factors computed for various values of $\gamma$; left – with respect to the
growth in GR, middle and right – with respect to nDGP for time-independent and
time-dependent growth index respectively. Dashed black lines denote the $1\%$
range.

Below is a summary of the functions required for the reaction calculations and
our assumptions:

  1. 1.

The background expansion is set to be equal the standard cosmology: $H(a)=\\\
H_{0}\sqrt{\Omega_{\rm m}a^{-3}+(1-\Omega_{\rm m})}$.

  2. 2.

Modifications in the linear regime:

$\mu_{\rm L}(a)=\begin{cases}1+\frac{1}{3\beta(a)}&\text{ for nDGP}\\\
\frac{2}{3}\Omega_{\rm m}^{\gamma-1}(a)\left[\Omega_{\rm
m}^{\gamma}(a)+2-3\gamma+3(\gamma-1/2)\Omega_{\rm m}(a)\right]&\text{ for
$\gamma=const$}\\\ \frac{2}{3}\Omega_{\rm
m}^{\gamma-1}(a)\gamma_{1}(a-1/a)\log{\Omega_{\rm
m}(a)}+\mu^{\gamma=const}_{\rm L}(a)&\text{ for $\gamma(a)$}\end{cases}$
(2.13)

with $\beta(a)$ defined as

$\beta(a)=1+\frac{H(a)}{H_{0}}\frac{1}{\sqrt{\Omega_{rc}}}\left(1+\frac{\dot{H}(a)}{3H^{2}(a)}\right)\,,$
(2.14)

where dot denotes a derivative with respect to the time coordinate and
$\Omega_{rc}$ corresponds to the strength of modification in the nDGP
cosmology ($\Omega_{rc}=0$ in GR). In addition to the standard (constant)
growth index parameter, we consider a time-dependent version of the form
$\gamma(a)=\gamma_{0}+\gamma_{1}(a+1/a-2)$ from the recent work of Ref. [86].
We add this model to MGCAMB999github.com/MariaTsedrik/MGCAMB. The latter
parametrisation was recently proposed because it captures the time-evolution
of growth in Horndeski theories better than the constant growth index. We
demonstrate this in the lower right panel of Fig. 1. Additionally, in the same
figure, we illustrate the impact different choices of $\mu_{\rm L}$ have on
the growth factor $D$ for the time-independent growth index. Indeed, for nDGP
the time-dependent parametrisation of the growth index provides a better fit
than the constant one. However, for the standard cosmology, the Linder gamma
parametrisation remains excellent, whereas, for a scalar-tensor theory like
nDGP, it does not. Further discussion on the correspondence between the
constant growth index and DGP theories is provided in Appendix A. In all
models, we ensure a match with General Relativity at high redshifts. In other
words, we introduce modifications of the standard cosmology only in the late
universe, while the early universe remains $\Lambda$CDM.

  3. 3.

In the quasi-nonlinear regime for nDGP $\gamma_{2}$\- and
$\gamma_{3}$-functions within $S(\mbox{\boldmath$k$},a)$ are specified in Ref.
[76], while for the scale-independent growth index parametrisation we set them
to zero. In Ref. [45] it was demonstrated that the 1-loop corrections impact
the reaction function very weakly for nDGP, less than $1\%$. However, this
might become a consideration when scale-dependency in the linear growth is
present.

  4. 4.

Modifications in the fully nonlinear regime:

$\mu_{\rm
NL}(a)=\begin{cases}\frac{2}{3\beta(a)}\frac{\sqrt{1+s^{3}}-1}{s^{3}}+1&\text{
for nDGP}\\\ (\mu^{\gamma}_{\rm L}(a)-1)\mathrm{erf}(a\,y_{\rm
h}10^{q_{1}})+1&\text{ for $\gamma$ parametrisations}\,.\end{cases}$ (2.15)

The expression for nDGP is derived in Ref. [87] following the solution of a
spherically symmetrical overdensity with

$s=\left[\frac{2\Omega_{\rm
m}(\delta+1)}{9a^{3}\beta(a)^{2}\Omega_{rc}}\right]^{\frac{1}{3}},$ (2.16)

and $\delta$ being the nonlinear over-density given by

$\delta=y_{\rm h}^{-3}(1+\delta_{\rm ini})-1\,,$ (2.17)

with $\delta_{\rm ini}$ being the initial over-density and

$y_{\rm h}\equiv\frac{R_{\rm TH}/a}{R_{i}/a_{\rm ini}}\,,$ (2.18)

$R_{\rm TH}$ and $R_{i}$ being the physical halo top-hat radius at the target
scale factor $a$ and at the initial scale factor $a_{\rm ini}$, respectively.
The second expression is phenomenological and was demonstrated in Ref. [45] to
reproduce the Vainshtein screening behaviour within $1\%$ up to
$k=5~{}h/\mathrm{Mpc}$ when compared against the exact solution. The parameter
$q_{1}$ corresponds to the screening scale, and its impact is demonstrated in
the right upper panel of Fig. 1. From this figure, we see that the higher
$q_{1}$ the deeper into the halo the transition from modification to no
modification happens. The opposite is also true, the lower $q_{1}$ the farther
outside the halo $\mu_{\rm NL}=1$. The GR-limit is fully recovered at all
scales when $q_{1}\rightarrow-\infty$. Values of $q_{1}\gtrsim 2$ correspond
to a screening scale deep within a halo, indicating that we will observe no
effect of screening in the nonlinear power spectrum at the scales of our
interest.

###  2.4 Baryonic Effects

The last missing component from Eq. 2.6 is the boost from the baryonic
feedback. It is well-known that astrophysical processes significantly impact
the matter power spectrum in the nonlinear regime. For instance, active
galactic nuclei (AGN), supernovae, and stellar winds repulse matter from
clustered centralised clumps into intergalactic medium, causing a suppression
of the observed structures on small scales. On even smaller astrophysical
scales, the observed structures get boosted by gas cooling and star formation.
The effects of baryonic processes have been extensively studied in various
hydrodynamical simulations, including COSMO-OWLS [88], BAHAMAS [89],
Illustris-TNG [90], and more recent FLAMINGO [91], among others. However,
these effects are sensitive to the sub-grid physics and can vary significantly
between different simulations (as depicted in Fig. 9 of Ref. [92]).

Instead of running computationally expensive hydro-simulations and fine-tuning
them to match observations, the impact of baryons can be modelled using a
baryonification model [93, 92]. In this approach, halo profiles, the outputs
of a dark-matter-only simulations, are modified by slightly displacing
particles around halo centres in a spherically symmetric way. This
displacement is characterised by 7 baryonic parameters:
${\log_{10}M_{c},\,\theta_{\rm
ej},\,\mu,\gamma,\,\delta,\,\eta_{\delta},\,\eta}$ (the first five describe
the gas distribution, while the last two are related to the stellar
abundances), along with one cosmological parameter, the baryonic fraction
$f_{\rm b}=\Omega_{\rm b}/\Omega_{\rm cdm+b}$. In our analysis we use the
baryonic emulator BCEMU [94], which accurately reproduces the power spectra of
state-of-the art hydrodynamical simulations even with a reduced number of
parameters. This emulator has been employed to analyse cosmic shear data in
Ref. [95], and provides a robust and flexible tool for incorporating baryonic
effects into our analyses. In Section 4.3 we discuss the cosmology dependence
in the context of MG theories for the baryonification approach.

The importance of baryonic feedback becomes evident on scales comparable to
the size of a halo, which coincides with the scales where screening is crucial
for MG theories. In this work, we aim to investigate whether an independent
parametrisation of the screening scale can be decoupled from the effects of
baryons. To achieve this, when exploring the generalised parametrisation in
Sections 4.3 and 4.5, we only vary $\log_{10}M_{c}$ (similarly to the DES Y3
shear analysis in Ref. [96]), while keeping the others fixed at their fiducial
values. Exploring degeneracies with the full set of baryonic parameters will
be explored in future work. The $\log_{10}M_{c}$ parameter controls the slope
of the gas density distribution: smaller haloes with masses less than $M_{c}$
have shallower profiles. Both the screening and baryonic feedback effects
display weak cosmological dependence and result in the suppression of growth
on approximately the same scales. Already from this fact alone we can
highlight the urgency of accurately measuring baryonic parameters.

Fortunately, there exists a well-established connection between the baryonic
effects on the power spectrum and the gas and stellar fractions in haloes
[93]. This suggests that additional CMB and/or X-ray measurements of galaxy
groups and clusters can help constrain these parameters or impose physically
motivated priors. For instance, in Ref. [97] authors measured
$\log_{10}{M_{c}}=14.53\pm 0.20$ with a compilation of Bayesian population
studies of galaxy groups and clusters and with cluster gas density profiles
derived from deep, high-resolution X-ray observations. While Ref. [98] found
from the combination of the cosmic shear data from KiDS-1000 with the gas
profiles from kinematic Sunyaev–Zeldovich (kSZ) observations and X-ray data
constrained $\log_{10}{M_{c}}\approx 13.2\pm 0.4$. Recently, Ref. [99] found a
significantly stronger baryonic feedback suppression than in previous studies
using the DES Y3 cosmic shear and CMASS + ACT DR5 kSZ observations with
$\log_{10}{M_{c}}=13.22^{+0.42}_{-0.29}$. In these three works all 7
parameters of the model were varied. Interestingly, it has been demonstrated
in Ref. [8] that with the addition of mock gas fractions from eROSITA a
Euclid-like setup finds
$\sigma_{\log_{10}{M_{c}}}/\log_{10}{M_{c}}\sim\mathcal{O}(0.001)$. This
implies that future surveys can decrease baryonification model errors by more
than an order of magnitude.

Alternative approaches to mitigate the baryonic effects include a) a slightly
different model of baryonification called Baryon Correction Model [100, 101]
that was applied to DES Y3 data [96, 9], b) halo-model [102, 71], and c)
Principal Component Analysis (PCA) [103, 104, 105].

##  3 Analysis Setup

The survey specifications in our pipeline are chosen to mimic Stage-IV
specifications: sky fraction $f_{\rm sky}=0.4$, number density of galaxies per
arcminute squared $n=30$ arcmin-2, and per-component dispersion in the
intrinsic galaxy ellipticities $\sigma_{\epsilon}=0.3$. We take 10
equipopulated tomographic bins between redshifts $z\in[0.001,2.5]$. The
photometric redshift distribution of galaxies is following
$n_{i}(z)\propto(z/z_{0})^{2}\exp{-\left(z/z_{0}\right)^{3/2}}$ with
$z_{0}=0.9/\sqrt{2}$ [106] . The uncertainty on the photometric redshift is
modeled as the sum of two Gaussian distributions: one for the well determined
photometric redshifts and another for the outliers (we follow the prescription
in Section 3.3.1 of Ref. [6]).

As input we use synthetic data, namely $C_{\ell}$ computed according to Eq.
2.1, with the fiducial cosmology in two scenarios:

  1. 1.

$\Lambda$CDM scenario: nonlinear power spectrum with HMcode2020 in which
linear input is computed with BACCOemu;

  2. 2.

nDGP scenario with strong deviation from the standard cosmology
$\Omega_{rc}=0.25$: nonlinear power spectrum computed with reaction from ReACT
and pseudo-power spectrum from HMcode2020.

We use these scenarios to validate the pipeline and then test the performance
of the growth index parametrisation. In total, we emulate 3 models using
cosmopower [43]: nDGP, $\gamma+q_{1}$ and $\gamma(z)+q_{1}$ in the ranges
specified in Table 2 and the accuracy of emulators provided in the
corresponding code repository101010nebblu/ReACT-emus/emulators/Accuracy Plots.
We compute our models in the wavenumber range of
$k\in[0.01,5]\,h/\mathrm{Mpc}$. For smaller $k$-values we extrapolate the
boost to be constant and equal to its value at $k=0.01\,h/\mathrm{Mpc}$, while
for larger $k$-values we extrapolate the boost as a power law.

The likelihood is computed as (see Appendix B.2 of Ref. [11]):

$\mathcal{L}=-\frac{1}{2}\sum_{\ell=10}^{\ell_{\rm max}}(2\ell+1)f_{\rm
sky}\left(\frac{d_{\ell}^{\rm mix}}{d_{\ell}^{\rm
model}}+\ln{\frac{d_{\ell}^{\rm model}}{d_{\ell}^{\rm data}}}-N_{\rm
bins}\right)$ (3.1)

where $N_{\rm bins}$ is the number of photometric bins and $d_{\ell}$ denote
the determinants of $C_{\ell}^{ij}$: the mock data and model are both computed
with Eq. 2.3, and $d_{\ell}^{\rm mix}$ is defined as

$d_{\ell}^{\rm mix}=\sum_{k=1}^{N_{\rm bins}}\det\left[\begin{cases}C^{\rm
model}_{ij}(\ell)&\text{if }j\neq k\\\ C^{\rm data}_{ij}(\ell)&\text{if
}j=k\end{cases}\right]\,.$ (3.2)

The modelling and likelihood function we use were validated in Ref. [107]. The
priors are given in Table 1, and the fiducial values for the fixed baryonic
parameters are: $\theta_{\rm ej}=4.235$, $\mu=0.93$, $\gamma=2.25$,
$\delta=6.4$, $\eta_{\delta}=0.14$, $\eta=0.15$.

Effect | Parameter | Prior | Fiducial
---|---|---|---
Cosmological | $\omega_{c}$ | $\mathcal{U}(0.11,0.13)$ | 0.12
$\omega_{b}$ | $\mathcal{N}(0.02268,0.00038)$ | 0.02268
$h$ | $\mathcal{U}(0.63755,0.7305)$ | 0.68
$\log{\left(10^{10}~{}A_{s}\right)}$ | $\mathcal{U}(2.7081,3.2958)$ | 3.044
$n_{s}$ | $\mathcal{N}(0.97,0.004)$ | 0.97
Intrinsic Alignments | $A^{\rm IA}$ | $\mathcal{U}(0,12.1)$ | 1.72
$\eta^{\rm IA}$ | $\mathcal{U}(-7,6.17)$ | -0.41
Modified Gravity | $\log_{10}\Omega_{rc}$ | $\mathcal{U}(-3,2)$ | 0.25
$\gamma,\gamma_{0}$ | $\mathcal{U}(0,1)$ | -
$\gamma_{1}$ | $\mathcal{U}(-0.7,0.7)$ | -
$q_{1}$ | $\mathcal{U}(-2,2)$ | -
Baryonic | $\log_{10}M_{c}$ | $\mathcal{U}(11.4,14.6)$ | 13.32
Massive Neutrinos | $M_{\nu}$ | $\mathcal{U}(0,0.1576)$ | 0.06
Table 1: Sampling parameters, prior ranges and fiducial values in the model and mock data respectively. The spectral index $n_{s}$ has a tight Planck priors and the baryonic density component has a BBN-prior, both motivated by weak constraints on these parameters from cosmic shear data [108].  | $\Omega_{\rm m}$ | $\Omega_{\rm b}$ | $\Omega_{\nu}$ | $h$ | $n_{s}$ | $A_{s}$ |  $\gamma$, $\gamma_{0}$ | $\gamma_{1}$ | $q_{1}$ | $\log_{10}\Omega_{rc}$ | $z$
---|---|---|---|---|---|---|---|---|---|---|---
lower | 0.2899 | 0.04044 | 0. | 0.629 | 0.9432 |  1.5 $\times 10^{-9}$ | 0. | -0.7 | -2 | -3 | 0.
upper | 0.3392 | 0.05686 | 0.00317 | 0.731 | 0.9862 |  2.7 $\times 10^{-9}$ | 1. | 0.7 | 3 | 2. | 2.4
Table 2: Prior ranges for the emulators.

The posterior distribution is sampled with Nautilus [109] for $3\times 10^{3}$
live points. This sampler is based on the importance nested sampling (INS)
technique combined with deep learning via neural network regression. Similarly
to other nested sampling-based samplers, Nautilus tends to under-predict the
size of the posterior contours when compared to much slower but more accurate
samplers. This is a known feature (see, for example, Appendix D of Ref. [110])
and not relevant for the purpose of this work. We varied the number of live
points to guarantee the convergence and to control the differences in
posterior-volumes when comparing between models with different number of
parameters.

##  4 Results

We first validate our models on synthetic data in $\Lambda$CDM, i.e. data with
GR as the model of gravity. Then we test how the index growth models – with
and without screening – perform when the true cosmology is nDGP, i.e., with
enhancement of structure on large scales and recovering of GR on small-scales
via Vainshtein screening. After that, we produce synthetic data with the
screened $\gamma+q_{1}$ model, and study degeneracies between the screening
scale and baryonic feedback, and massive neutrinos. Finally, we compare the
performance of the exact nDGP model versus the model-independent growth index
approach on a mock data vector assuming nDGP gravity, baryonic feedback and
massive neutrinos.

###  4.1 Validation with GR

Figure 2: Marginalised posterior distribution for the cosmological parameters
for the two model choices, as detailed in the legend. We fit the cosmic shear
power spectra with $\ell_{\rm max}=3000$. Solid lines and filled contours
correspond to $B^{\rm baryons}=1$ from Eq. 2.5 (for both – data and models),
while dashed lines and empty contours correspond to the baryonic feedback
contribution with varying $M_{c}$ and fixed values of the other baryonic
parameters as detailed in the main text. Grey lines mark the true values of
the synthetic data.

Figure 3: Left panel: variation in the ratio $P^{\rm MG}_{\rm
NL}/P^{\Lambda\rm CDM}_{\rm NL}$ at redshift zero for the screened,
$\gamma+q_{1}$, and unscreened, $\gamma$, models drawn from the posterior
distributions in Fig. 2 within the $1\sigma$ deviation around the posterior
maximum for $\gamma$. Right panel: marginalised posterior distribution in
$\Omega_{\rm m}-\gamma$ for the same models on the mock GR data with
$\ell_{\rm max}=3000$ (filled contours) and $\ell_{\rm max}=1000$ (unfilled
contours).

In order to validate our likelihood pipeline, we perform MCMC analyses on mock
data in the standard cosmology (without massive neutrinos). For the scale-cut
we take $\ell_{\rm max}=3000$, which is considered a middle value between
“pessimistic” ($\ell_{\rm max}=1500$) and “optimistic” ($\ell_{\rm max}=5000$)
scenarios for Stage-IV surveys [6]. We test two models – gamma without
screening via a pseudo-power spectrum ($\gamma$) and gamma with screening
($\gamma+q_{1}$). In Fig. 2 we show posterior distributions for the main
parameters of interest in weak lensing as well as the expansion rate today,
$h$. The latter is added to demonstrate its anti-correlation with the extended
parameter $\gamma$, while the $h$ degeneracy with $\omega_{b}$ and $n_{s}$ is
broken by informative priors on the latter two parameters. The full posterior
distribution is demonstrated in Fig. 10 of Appendix B together with a detailed
discussion on the various degeneracies and the parameters controlling the
amplitude111111For the amplitude of the power spectrum we present both – the
sampled parameter of the primordial amplitude $A_{s}$, and its derived late-
Universe counterpart $\sigma_{8}=D_{\gamma}(z=0)/D_{\Lambda\rm
CDM}(z=0)\,\sigma^{\Lambda\rm CDM}_{8}$, with $\sigma_{8}$ being the r.m.s.
density variation when smoothed with a tophat-filter of radius
$8~{}\mathrm{Mpc}/h$, and $S_{8}=\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$. We also
show the derived matter density parameter $\Omega_{\rm
m}=(\omega_{c}+\omega_{b})/h^{2}$..

In Fig. 2, both models recover all cosmological parameters correctly, and the
value of gamma is $\sim 0.55$ as expected, with uncertainties
$\sigma_{\gamma}=0.03$ for $\gamma+q_{1}$ and $\sigma_{\gamma}=0.04$ for
$\gamma$ without screening. The fact that a model with more parameters results
in tighter constraints is explained by the functional form of $\mu_{\rm NL}$
(see Eq. 2.15): any deviation from $\gamma\sim 0.55$ or $\mu_{\rm L}=1$ opens
the possibility of constraining $q_{1}$, which leads to a worse fit or smaller
likelihood due to nonlinear signatures that are very different from the GR
predictions. We illustrate this in the left panel of Fig. 3, where we plot the
power spectra ratio $P^{\rm MG}_{\rm NL}/P^{\Lambda\rm CDM}_{\rm NL}$ at $z=0$
from the posterior distribution’s $1\sigma$ $\gamma$ constraints. We see that
broader uncertainty on $\gamma$ in the pseudo model manifests on linear
scales, however for $k>1~{}h/\mathrm{Mpc}$ the variation in the power spectrum
becomes larger for $\gamma+q_{1}$ than for $\gamma$ only. This also implies
that for more conservative scale-cuts we should not obtain differences in
$\sigma_{\gamma}$, this effect is purely due to the inclusion of nonlinear
scales and high sensitivity to these scales in a Stage-IV-like setup. As an
additional proof of this conclusion we demonstrate the $\Omega_{\rm m}-\gamma$
contours for stricter scale-cuts of $\ell_{\rm max}=1000$ in the right panel
of Fig. 3. The $1\sigma$ ranges are visibly in better agreement between the
models now. The inclusion of a baryonic feedback model increases the errors
for the pseudo model to $\sigma_{\gamma}=0.04$ and $\sigma_{\gamma}\approx
0.05$ for the screened model. In this case the strong degeneracy between the
screening scale and baryonic parameters explains a more significant broadening
of the contours. We study this degeneracy in detail in Section 4.3.

###  4.2 Screening Scale with nDGP Gravity

Figure 4: Left panel: Power spectrum ratios at redshift $z=0$ for screened
(green lines) and unscreened (orange line) models with nonlinear prescriptions
described in the text. Grey line denotes $k_{\rm max}=5~{}h/\mathrm{Mpc}$ that
we compute and emulate with ReACT, for higher values we interpolate the MG
boost as a power law. Right panel: Corresponding ratios of the shear angular
power spectra. Central redshifts: for bin 1 $z_{c}=0.21$, for bin 2
$z_{c}=0.49$. We fix all cosmological parameters to the fiducial values.

Next we validate our model on a theory with the Vainshtein-screening mechanism
on nonlinear scales and a strong deviation from GR on linear scales ($\sim
14\%$ more structure than in GR on linear scales at lower redshifts). In the
left panel of Fig. 4 we show the ratio of power spectra $P^{\rm MG}_{\rm
NL}/P^{\Lambda\rm CDM}_{\rm NL}$ at $z=0$ for the nDGP model with
$\Omega_{rc}=0.25$ (black solid line). For the same fiducial cosmological
parameters, we show also this ratio for the screened (green line) and pseudo
(orange line) growth index models. The pseudo approach captures the bump-
feature in the range of $k\sim 0.2-0.9~{}h/$Mpc but overpredicts its
amplitude. The screened model reproduces this feature extremely well, at the
percent level for $k\lesssim 3~{}h/\rm{Mpc}$. In the right panel, we show
ratios of the corresponding angular power spectra, $C_{\ell}$, for the 3
models with respect to GR for the same set of cosmological parameters. The
bump-feature in the nonlinear regime is visible in cross-correlated redshift-
bins, while in the auto-correlated bins we see characteristic curves that
tightly converge for all models at smaller scales. However, we notice a shift
at linear scales for the $\gamma$ models. The reason for this shift is the
difference in the time evolution between nDGP and the growth index
parametrisation as discussed in Section 2.3: while the growth modification is
matched between different models at $z=0$, this is no longer the case for
higher redshifts.

Figure 5: Marginalised posterior distribution for the cosmological parameters
for the nDGP cosmology and the two model and scale-cut choices, as detailed in
the legend. The baryonic boost is set to unity. Grey solid lines mark the true
values of the synthetic data, grey dashed lines mark the GR-limit. Green
dashed lines denote the best-fit values from the MCMC chain for $\gamma+q_{1}$
model.

We consider the same two models: models with screening in the reaction
function ($\gamma+q_{1}$) and without screening via a pseudo-power spectrum
($\gamma$) and repeat our analysis on the mock nDGP data with $\ell_{\rm
max}=3000$ and using the same setup. We show the marginalised posterior
distribution in Fig. 5. The full posterior distribution for the sampled
cosmological parameters is shown in Fig. 13 of Appendix D with a more detailed
description. Overall, our model with screening (green contours and lines)
correctly captures the fiducial cosmology, detects a $\gamma$-value lower than
its GR-limit as expected, namely $\gamma=0.38\pm 0.03$, and detects a
screening transition with $q_{1}=0.29\pm 0.10$. We notice a $1\sigma$-bias in
$\Omega_{\rm m}$ and a $2\sigma$-bias in $S_{8}$. Both are consequences of
fitting a model with the different time-evolution evolution of growth to data
with small error-bars. We check whether this can also be connected to
projection effects [52]. In other words whether the likelihood maxima (green
dashed lines) are at the true values of cosmological parameters and not at the
posterior maximum values. From the MCMC chain, we derive the best-fit value
$\Omega^{\rm best-fit}_{\rm m}\approx 0.306$, which is closer to the fiducial
value of $\Omega^{\rm fid}_{\rm m}\approx 0.309$ than the posterior mean of
$\Omega^{\rm mean}_{\rm m}\approx 0.304$. Also the best-fit values for
$\gamma^{\rm best-fit}\approx 0.42$ and $S_{8}^{\rm best-fit}\approx 0.898$
very close to the true value of $S_{8}^{\rm nDGP}\approx 0.895$. We use the
likelihood minimiser minuit121212github.com/jpivarski/pyminuit and obtain
similar values for the best-fit parameters. This implies a presence of
projection or prior-volume effects. Lastly, we note that this bias vanishes
when more parameters are varied and the contours are broadened (see Section
4.5). Overall, the agreement within $2\sigma$ satisfies the “precision versus
accuracy” test for our model.

Figure 6: Marginalised posterior distribution for the cosmological and
extended parameters for the analysis with the fiducial nDGP cosmology using
the cosmic shear power spectrum with $\ell_{\rm max}=3000$. We fit the time-
dependent growth index with the screening scale, with and without fixing a
value of $\gamma_{1}$, as detailed in the legend. Grey solid lines mark the
true values of the synthetic data, grey dashed lines mark the GR-limit. Orange
dashed lines denote the best fit values from the MCMC chain for the second
model.

In Fig. 5, we notice that even for scale-cuts below the pessimistic scenario
for Stage-IV surveys, namely for $\ell_{\rm max}=1000$, there is a significant
bias in $S_{8}$ and $\Omega_{\rm m}$ for the model without screening (orange
contours and lines), with $S_{8}$ being closer to its GR value. Moreover, the
value of $\gamma$ does not correspond to the expected behaviour on large
scales but rather points towards suppression of structure formation relative
to the standard cosmology. This signals that ignoring a correct screening
implementation will result in the wrong extracted cosmology.

In the nonlinear regime our model-independent screening scale seems to be in a
good agreement with nDGP theory which includes Vainshtein screening
mechanisms. The disadvantage of this model is the inaccurate representation of
time evolution of growth on linear scales. We can improve this by introducing
the time-dependent gamma-parametrisation from Ref. [86] (also see Section
2.3). The results are shown in Fig. 6 and Fig. 14 for the same setup as
before. From the green contours in Fig. 6 we see that while fiducial values
lay within $1\sigma$ of all parameters, the screening scale is not properly
constrained with cosmic shear information alone. Note how in general, extended
parameters make the posterior distribution non-Gaussian by bringing additional
degeneracies with cosmological parameters and with each other. Two extended
parameters $\gamma_{0}$ and $\gamma_{1}$ are strongly anti-correlated, so we
can aim to break degeneracies by fixing $\gamma_{1}$. We find from fitting to
the boost at lower redshifts that $\gamma_{1}=-0.19$ is a good fit and it
breaks the degeneracy and allows us to constrain the screening parameter
$q_{1}=0.57^{+0.07}_{-0.12}$ with $\gamma_{0}=0.47\pm 0.02$ (orange contours
and lines). In Fig. 6, we clearly see that while $\Omega_{\rm m}$ is still
$1\sigma$-biased towards lower values, the bias in $S_{8}$ towards higher
values decreases to $1\sigma$ contrary to the time-independent growth index.
We again compute and plot best-fit values (dashed orange lines) from the MCMC
chain, for fixed $\gamma_{1}$ they are $\Omega^{\rm best-fit}_{\rm m}\approx
0.304$ ($0.306$ from minuit), $\gamma^{\rm best-fit}\approx 0.49$, and
$S_{8}^{\rm best-fit}\approx 0.894$. The model with more accurate time-
evolution of growth seems to decrease the projection effects. The necessity of
tight constraints on $\gamma_{1}$ implies that additional measurements from
clustering will improve our constraints. We aim to investigate probe
combinations in future work.

###  4.3 Screening Scale versus Baryons

Figure 7: Left panel: the power spectrum boost at redshift zero for the MG
contribution, $\gamma=0.4$ and varying $q_{1}$ (green lines), and baryonic
feedback (orange lines). Right panel: marginalised posterior distributions
from the analysis with the fiducial $\gamma+q_{1}$ cosmology using the cosmic
shear power spectrum with $\ell_{\rm max}=3000$ for the MG and baryonic
parameters with different priors on the latter, as detailed in the legend.
Grey lines mark the true values of the synthetic data.

Here we explore the degeneracy between the suppression due to the screening
mechanism and baryonic feedback effects. In the left panel of Fig. 7 we
demonstrate the comparison between MG and baryonic boosts for different
parameters. The larger the value of $q_{1}$, the less is the suppression due
to screening. The larger the value of $\log_{10}{M_{c}}$, the stronger is the
suppression due to baryonic feedback at nonlinear scales. It is clear that
both effects become noticeable at the same scales $k\gtrsim
1~{}h/\mathrm{Mpc}$, since they demonstrate similar scale-dependence and
overall amplitude. Already from this we can conclude that our model-
independent screening scale, $q_{1}$, is highly degenerate with the baryonic
suppression parameter $\log_{10}{M_{c}}$.

In order to avoid biases due to different time evolution of growth structure
on linear scales, we perform an MCMC-analysis on the mock data computed with
the $\gamma+q_{1}$ model directly with $\gamma^{\rm fid}=0.4$ and $q_{1}^{\rm
fid}=0.76$. We show the marginalised posterior distribution for $\gamma-
q_{1}-\log_{10}{M_{c}}$ in the right panel of Fig. 7, while the corresponding
full posterior distribution is shown in Fig. 15. The unbiased determination of
$\gamma$, as well as other cosmological parameters, is not affected by the
inclusion of the baryonic feedback. However, the screening scale demonstrates
a very strong positive correlation with the baryonic parameter and prefers the
unscreened regime of $q_{1}>1.5$. The baryonic parameter absorbs the
suppression due to the screening and as a result gets biased towards higher
values with respect to the fiducial value. This degeneracy persists, when
testing variations in the fiducial values of $\log_{10}{M_{c}}$ within the
same $\gamma+q_{1}$ model, while also examining different fiducial values of
$\gamma+q_{1}$ while keeping $\log_{10}{M_{c}}$ constant. We perform these
checks to mitigate potential coincidences in the choice of fiducial
parameters. However, we find that when we impose a tight flat prior on
$\log_{10}{M_{c}}$, the screening scale is recovered unbiased. This means we
require a good understanding of baryonic physics for the detection of the
model-independent screening transition: for instance, tight priors on baryonic
parameters of the order $\sigma_{\log_{10}{M_{c}}}=0.1$ or
$\sigma_{\log_{10}{M_{c}}}/\log_{10}{M_{c}}\approx 1\%$.

Instead of imposing priors, for example, based on independent galaxy cluster
X-ray observations, a better approach will be the cross correlation or the
joint analyses of weak lensing and X-ray/kSZ probes like in the works
mentioned in Section 2.4. We also highlight the importance of non-standard
hydro-dynamical simulations. The BCEMU emulator we use in this work was
trained for standard cosmologies. Despite the flexibility of the
baryonification model with 7 or 3 parameters, we still do not know whether
this is accurate enough to model nonlinear baryonic physics in
beyond-$\Lambda$CDM cosmologies. Reassuringly, Refs. [69, 111] show both
within the halo model and in $N$-body simulations that baryons and MG physics
can be modelled independently to a large extent. A promising argument for this
would be that the astrophysics responsible for baryonic effects is happening
on the scales in the screened regime, hence there is no reason to believe that
it is significantly different from the standard baryon feedback processes.
Especially since these are only weakly cosmology dependent (up to several
percent for small variations in cosmology) and mainly depend on the fraction
of baryons to the total matter according to various studies (e.g., Refs. [112,
113]). However, more work is needed to confirm whether these assumptions
satisfy the requirements of Stage-IV surveys [114].

###  4.4 Screening Scale versus Neutrinos

Figure 8: Left panel: the power spectrum boost at redshift zero for the MG
contribution with $\gamma=0.4$, $M_{\nu}=0$ and $q_{1}$ varied (green lines),
and $\gamma=0.4$, $q_{1}=0.76$, $M_{\nu}$ varied (orange lines). The boosts
for various neutrino masses are normalised to the value at the largest scales
computed. Right panel: marginalised posterior distributions from the analysis
with the fiducial $\gamma+q_{1}$ cosmology using the cosmic shear power
spectrum with $\ell_{\rm max}=3000$ for models with different cosmological
parameters varied, as detailed in the legend. Grey lines mark the true values
of the synthetic data.

The last nonlinear effect to consider is the suppression of structure growth
due to the contribution of massive neutrinos. In the left panel of Fig. 8 we
demonstrate the comparison between MG boosts for different values of screening
and neutrino mass parameters. The larger the value of $q_{1}$, the smaller the
suppression due to screening. The larger the value of $M_{\nu}$, the stronger
is the suppression due to massive neutrinos. However, massive neutrinos lead
to suppression of growth already at mildly nonlinear scales $k\gtrsim
0.1~{}h/\mathrm{Mpc}$. Furthermore, the scale-dependency in the suppression
due to massive massive neutrinos is drastically different form the suppression
due to screening: there is no clear feature in the power spectrum in a
constrained scale range. We still expect a positive correlation between
$q_{1}$ and $M_{\nu}$, but this will be in addition to noticeable degeneracies
with other cosmological parameters controlling the amplitude and slope of the
power spectrum at all scales, and not only in the nonlinear regime. For
instance, the primordial amplitude $A_{s}$, the expansion rate $h$, and the
matter density $\Omega_{\rm m}$. In the right panel of Fig. 8, for the same
fiducial $\gamma+q_{1}$ cosmology as in the previous section with the addition
of massive neutrinos with $M_{\nu}^{\rm fid}=0.06$ eV, we showcase how the
more parameters we vary the more sensitivity to the neutrino mass we are
losing. For example, the prominent and expected negative correlation between
$\gamma$ and $M_{\nu}$ (pink contour) disappears even if just one additional
parameter is varied (purple contour). We conclude that Stage-IV cosmic shear
measurements alone are not sufficient to put any constraints on the neutrino
mass, which is in agreement with the findings of other works, for example Ref.
[29].

###  4.5 Combined Nonlinear Effects

Figure 9: Marginalised posterior distribution for the cosmological and MG
parameters for the analysis of the fiducial nDGP cosmology with massive
neutrinos and baryonic feedback with a tight uniform prior on the baryonic
parameter, using the cosmic shear power spectrum up to $\ell_{\rm max}=3000$.
For nDGP and $\gamma+q_{1}$ the prior on $\log_{10}{M_{c}}$ is
$\mathcal{U}(11.4,14.6)$, while for the $\gamma(z)+q_{1}$ model it is
$\mathcal{U}(13.3,13.34)$. Grey lines mark the true values of the synthetic
data.

Finally, we combine our findings in the previous sections to demonstrate how
our model-independent approach performs if the underlying cosmology is nDGP
with massive neutrinos (with $M_{\nu}=0.06$ eV) and known baryonic feedback
effects (within some prior). All background properties and the growth factor
in the IA terms are computed with the total matter density including neutrinos
via $\Omega_{\rm m}$, while for the baryonification, since it was trained
without neutrinos and is relevant on nonlinear scales, we exclude the
neutrinos contribution in $f_{b}=\Omega_{\rm b}/(\Omega_{\rm
m}-\Omega_{\nu})$.

In our main result of this work, Fig. 9 (the full posterior distribution can
be found in Appendix E), we compare our model independent approach (green
colour for time-independent $\gamma$, orange colour for time-dependent
$\gamma$) against the exact MG model (purple colour). Overall, for the
fiducial value of $\log_{10}{\Omega_{rc}}=-0.6$, we get the mean value of
$\log_{10}{\Omega_{rc}}=-0.86^{+0.69}_{-0.36}$, which corresponds to
$\gamma=0.44^{+0.06}_{-0.05}$ with $q_{1}=0.33^{+0.21}_{-0.17}$ or
$\gamma_{0}=\gamma=0.47^{+0.03}_{-0.04}$ with $\gamma_{1}=-0.2$ and
$q_{1}=0.58^{+0.08}_{-0.14}$. The differences in the $S_{8}$ uncertainties
between models are due to the strong $\gamma$-dependence of the growth at
lower redshifts (this is discussed in Appendix B). However, the inferred value
of $S_{8}$ agrees between all models. This shows that the growth index
parametrisation seems to be a reasonable approximation in the late universe.
It is in the connection to the earlier universe where we see discrepancies in
the primordial amplitude, due to different growth evolution. For the time-
dependent growth index we infer slightly biased $\Omega_{\rm m}$ but more
accurate $\log{10^{10}~{}A_{s}}$. The tighter constraint on $\Omega_{\rm m}$
with $\gamma(z)$ is due to tighter prior on the baryonic parameter. In
general, we notice that the uncertainties on cosmological parameters are in a
good agreement between the exact and model-independent approaches. Therefore,
we advocate that our model-independent approach performs as well as the exact
modelling, but is more general.

In Fig. 9 we also notice that $\log_{10}M_{c}$ is biased towards higher values
in all models. It is highly degenerate with the extended parameters
$\log_{10}{\Omega_{rc}},~{}\gamma$, as well as with $\Omega_{\rm
m},~{}h,~{}\log{10^{10}~{}A_{s}}$. We impose an informative prior on the
baryon feedback $[\log_{10}{M^{\rm fid}_{c}}-0.2,\log_{10}{M_{c}^{\rm
fid}}+0.2]$ for nDGP and $\gamma+q_{1}$ based on the investigation from
Section 2.4. For the $\gamma(z)+q_{1}$ model, we require tighter priors on the
baryonic feedback and screening parameters: $[\log_{10}{M^{\rm
fid}_{c}}-0.02,\log_{10}{M_{c}^{\rm fid}}+0.02]$ and $[-1,1]$ respectively.
Otherwise their prominent degeneracy is weakening the constraints of
cosmological parameters, and the screening transition is poorly constrained
while preferring larger values closer to the unscreened limit. The total
neutrino mass is not constrained in the range of values that we emulated.

For completeness we list forecasted constraints on the nDGP parameter
$\log_{10}{\Omega_{rc}}$ from the literature. In Ref. [115], for a setup
similar to ours and with survey area of $15,000~{}\mathrm{deg}^{2}$, the
authors get from combined cosmic shear, photometric clustering and their
cross-correlation analysis with the same fiducial value of
$\log_{10}{\Omega_{rc}}$: $\sigma_{\log_{10}{\Omega_{rc}}}=0.3$ (pessimistic)
and $\sigma_{\log_{10}{\Omega_{rc}}}=0.12$ (optimistic). However, this study
does not include baryonic feedback effects, which we demonstrated to be
degenerate with the extended parameters. In Ref. [26], for a different setup
(5 redshift bins, number of galaxies per arcminute squared per tomographic bin
$\bar{n}=6$, similar shape noise, survey area $5,000~{}\mathrm{deg}^{2}$) and
for a comparably strong modification $\Omega_{rc}=0.36$, the authors find in
the optimistic scenario for cosmic shear
$\sigma_{\log_{10}{\Omega_{rc}}}=0.07$, varying three cosmological parameters
($\Omega_{\rm m},~{}S_{8},~{}h$), and omitting baryonic effects.

##  5 Conclusions

The growth index parametrisation is a single parameter extension of the
standard cosmology which allows for deviations in the linear growth functions.
Originally developed in the context of spectroscopic measurements, it serves
as an indicator of the detection of modified gravity (MG) theories. Previous
studies have demonstrated its effectiveness for a few close-positioned
redshift bins and up to mildly nonlinear scales [52]. In this work, we
presented an analysis with this parametrisation extended to the nonlinear
regime with a model-independent screening parametrisation, and applied this to
cosmic shear forecasts. For weak lensing, most of the information comes from
the small, nonlinear scales and is integrated over a broad range of redshifts,
in our case $z\in[0.001,2.5]$. Our theoretical model for the shear angular
power spectrum is based on the halo-based parameterisation introduced in Ref.
[45] and was emulated with the cosmopower emulator. We also take advantage of
the following emulators to accelerate computation of the nonlinear matter
power spectrum model: BCEMU for baryonic contribution and BACCOemu [116] for
the linear power spectrum as an input to HMcode.

To validate our model, we first tested it on noiseless mock data generated in
the standard ($\Lambda$CDM) cosmology in a Stage-IV setup with a scale cut
$\ell_{\rm max}=3000$. We successfully recovered unbiased cosmological
parameters and found the expected uncertainty on the growth index, $\gamma$,
to be $\sigma_{\gamma}\approx 0.05$, noting a slight asymmetry in the errors
when one baryonic feedback parameter is included. For a model without
screening and employing standard structure formation rescaled to match the
modifications on large scales, the obtained constraint was
$\sigma_{\gamma}=0.04$. We further compared the performance and differences
between these two $\gamma$-models in the absence of baryonic feedback. We
found that the form of our screening model at the nonlinear scales results in
tighter constraints on the growth index.

We then conducted a similar analysis on noiseless mock data computed in an
nDGP cosmology. This MG theory exhibits enhanced structure growth on linear
scales that transitions to GR behaviour on nonlinear scales via a screening
mechanism. The analysis was performed on mock data with a rather strong
deviation from GR (by $14\%$ at low redshifts in the linear regime). We aimed
to demonstrate the robustness of our modelling in this relatively extreme
scenario to guarantee its validity in more realistic scenarios with a weaker
deviation form the standard cosmology. We found that our model with screening
$\gamma+q_{1}$ performs well: it successfully recovers cosmological parameters
within $2\sigma$, finds $\gamma\sim 0.38$ lower than its GR-limit $\gamma\sim
0.55$, and detects a screening transition with its associated parameter
$q_{1}\sim 0.29$.

We also obtained a bias towards higher values in the primordial amplitude
$A_{\rm s}$. This bias arises from differences in the time evolution of the
linear growth factors between the growth index parametrisation and the nDGP
model. This inaccurate representation of the time evolution of structure
growth for scalar-tensor theories via a constant growth index parametrisation
has been pointed out in Ref. [86]. To address this issue, we explored a time-
dependent functional form $\gamma(z)=\gamma_{0}+\gamma_{1}z^{2}/(1+z)$.
Incorporating this time-dependent growth index, we found that the bias in the
amplitude vanishes. However, the constraints on the expansion rate and
screening transition are weakened significantly when we considered this model.
By fixing $\gamma_{1}=-0.19$, we broke the corresponding degeneracies, and we
recovered unbiased cosmological parameters with $\gamma_{0}\sim 0.47$ and
$q_{1}\sim 0.57$.

We found that ignoring the screening transition leads to biases in the
expansion rate and matter density, as well as to a false detection of $\gamma$
that exceeds its GR-limit. This happens even with scale cuts lower than the
“pessimistic” scenario in Ref. [6], with $\ell_{\rm max}=1500$. This
demonstrates the importance of correct inclusion of a screening scale when
extended cosmologies are considered.

We proceeded to study the degeneracies between the screening transition
parameter $q_{1}$, baryonic feedback, and massive neutrinos. All three effects
are nonlinear and result in the suppression of structure formation at small
scales. As before, we explored the case with only one free baryonic feedback
parameter $M_{c}$, which controls the slope of the gas distribution. We found
a strong positive correlation between the screening and baryonic parameters.
The suppression due to the screening was absorbed by the baryonic feedback if
no priors on the latter were imposed. Both effects are prominent in the same
range of scales and are highly degenerate. For this reason, we conclude that
using cosmic shear alone a detection of the model-independent screening
transition is possible only if tight priors on the baryonic parameters are
imposed. We used the uniform prior $[\log_{10}{M^{\rm
fid}_{c}}-0.2,\log_{10}{M_{c}^{\rm fid}}+0.2]$ motivated by cluster
measurements in Ref. [97]. Alternatively, a cross-correlation of cosmic shear
with X-ray/kSZ observations can break this degeneracy too. With the total
neutrino mass as a free parameter, we found that it cannot be constrained with
the cosmic shear data alone. This is in agreement with findings in Stage-III
surveys [117] as well as in forecasts for Stage-IV experiments [29]. We also
found no strong degeneracy between the screening transition and neutrino mass
– while they have similar impact on structure formation, their scale-
dependence and strength of impact differ.

Combining all the aforementioned nonlinear effects, we found that our model-
independent approach performs well when compared against the exact modelling,
and derived the following constraints in the full analysis on nDGP data with
massive neutrinos and baryonic feedback with a narrow flat prior:
$\sigma_{\log_{10}{\Omega_{rc}}}\approx 0.53~{}(88\%)$,
$\sigma_{\gamma}\approx 0.06~{}(13\%)$, $\sigma_{q_{1}}\approx 0.19~{}(58\%)$.

To conclude, we outline the next necessary steps in preparation for a fully
model-independent analysis for beyond-$\Lambda$CDM cosmologies with Stage-IV
cosmic shear surveys. In order to reduce the error bars on extended and
baryonic parameters, we advocate for the combination of cosmic shear with
photometric galaxy clustering and the corresponding cross-correlation, i.e.
the 3$\times$2-point analysis. Similarly, a combination with spectroscopic
galaxy clustering can alleviate the degeneracies and yield tighter constraints
[6]. After demonstrating the robustness of the $\gamma(z)$ approach, we aim to
explore a more agnostic approach, that would not require any assumption of the
time-evolution in the linear growth. For example, we can bin $\mu(z)$ from Eq.
2.13 in redshift directly [38, 35]. The goal would be to find an optimal
binning scheme. Additionally, as demonstrated in Ref. [45], we can extend our
modelling to include not only MG theories but also exotic dark energy models.

## Acknowledgments

MT would like to thank Alessandra Silvestri for useful comments and
suggestions, as well as for the warm hospitality in Leiden University, where a
large part of this project was accomplished. MT’s research is supported by a
doctoral studentship in the School of Physics and Astronomy, University of
Edinburgh. BB was supported by a UK Research and Innovation Stephen Hawking
Fellowship (EP/W005654/2). PC’s research is supported by grant RF/ERE/221061.
AP is a UK Research and Innovation Future Leaders Fellow [grant MR/X005399/1].
For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## References

  * [1] EUCLID collaboration, _Euclid Definition Study Report_ , 1110.3193.
  * [2] LSST collaboration, _LSST: from Science Drivers to Reference Design and Anticipated Data Products_ , _Astrophys. J._ 873 (2019) 111 [0805.2366].
  * [3] M. Kilbinger, _Cosmology with cosmic shear observations: a review_ , _Reports on Progress in Physics_ 78 (2015) 086901 [1411.0115].
  * [4] M. Bartelmann and M. Maturi, _Weak gravitational lensing_ , _Scholarpedia_ 12 (2017) 32440 [1612.06535].
  * [5] R. Mandelbaum, _Weak lensing for precision cosmology_ , _Ann. Rev. Astron. Astrophys._ 56 (2018) 393 [1710.03235].
  * [6] Euclid collaboration, _Euclid preparation: VII. Forecast validation for Euclid cosmological probes_ , _Astron. Astrophys._ 642 (2020) A191 [1910.09273].
  * [7] M. Takada and B. Jain, _Cosmological parameters from lensing power spectrum and bispectrum tomography_ , _Mon. Not. Roy. Astron. Soc._ 348 (2004) 897 [astro-ph/0310125].
  * [8] A. Schneider, A. Refregier, S. Grandis, D. Eckert, N. Stoira, T. Kacprzak et al., _Baryonic effects for weak lensing. Part II. Combination with X-ray data and extended cosmologies_ , _JCAP_ 04 (2020) 020 [1911.08494].
  * [9] G. Aricò, R.E. Angulo, M. Zennaro, S. Contreras, A. Chen and C. Hernández-Monteagudo, _DES Y3 cosmic shear down to small scales: constraints on cosmology and baryons_ , _Astron. Astrophys._ 678 (2023) A109 [2303.05537].
  * [10] J. Hamann, S. Hannestad and Y.Y.Y. Wong, _Measuring neutrino masses with a future galaxy survey_ , _JCAP_ 11 (2012) 052 [1209.1043].
  * [11] B. Audren, J. Lesgourgues, S. Bird, M.G. Haehnelt and M. Viel, _Neutrino masses and cosmological parameters from a Euclid-like survey: Markov Chain Monte Carlo forecasts including theoretical errors_ , _JCAP_ 01 (2013) 026 [1210.2194].
  * [12] P. Brax, _Screening mechanisms in modified gravity_ , _Class. Quant. Grav._ 30 (2013) 214005.
  * [13] C. Burrage and J. Sakstein, _Tests of Chameleon Gravity_ , _Living Rev. Rel._ 21 (2018) 1 [1709.09071].
  * [14] K. Koyama, _Gravity beyond general relativity_ , _Int. J. Mod. Phys._ D27 (2018) 1848001.
  * [15] V. Vardanyan and D.J. Bartlett, _Modeling and Testing Screening Mechanisms in the Laboratory and in Space_ , _Universe_ 9 (2023) 340 [2305.18899].
  * [16] C.M. Will, _The Confrontation between General Relativity and Experiment_ , 2014.  10.12942/lrr-2014-4.
  * [17] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, _Modified Gravity and Cosmology_ , _Phys. Rept._ 513 (2012) 1 [1106.2476].
  * [18] A. Joyce, L. Lombriser and F. Schmidt, _Dark Energy Versus Modified Gravity_ , _Ann. Rev. Nucl. Part. Sci._ 66 (2016) 95 [1601.06133].
  * [19] X. Liu et al., _Constraining $f(R)$ Gravity Theory Using Weak Lensing Peak Statistics from the Canada-France-Hawaii-Telescope Lensing Survey_, _Phys. Rev. Lett._ 117 (2016) 051101 [1607.00184].
  * [20] KiDS collaboration, _KiDS-1000 Cosmology: Constraints beyond flat $\Lambda$CDM_, _Astron. Astrophys._ 649 (2021) A88 [2010.16416].
  * [21] DES collaboration, _Probing gravity with the DES-CMASS sample and BOSS spectroscopy_ , _Mon. Not. Roy. Astron. Soc._ 509 (2021) 4982 [2104.14515].
  * [22] DES collaboration, _Dark Energy Survey Year 3 results: Constraints on extensions to $\Lambda$CDM with weak lensing and galaxy clustering_, _Phys. Rev. D_ 107 (2023) 083504 [2207.05766].
  * [23] A. Barreira, S. Bose, B. Li and C. Llinares, _Weak lensing by galaxy troughs with modified gravity_ , _JCAP_ 02 (2017) 031 [1605.08436].
  * [24] B. Li and M. Shirasaki, _Galaxy–galaxy weak gravitational lensing in $f(R)$ gravity_, _Mon. Not. Roy. Astron. Soc._ 474 (2018) 3599 [1710.07291].
  * [25] F. Hassani, J. Adamek, M. Kunz and F. Vernizzi, _$k$ -evolution: a relativistic N-body code for clustering dark energy_, _JCAP_ 12 (2019) 011 [1910.01104].
  * [26] J. Harnois-Déraps, C. Hernandez-Aguayo, C. Cuesta-Lazaro, C. Arnold, B. Li, C.T. Davies et al., _mglens: Modified gravity weak lensing simulations for emulation-based cosmological inference_ , _Mon. Not. Roy. Astron. Soc._ 525 (2023) 6336 [2211.05779].
  * [27] F. Schmidt, _Weak Lensing Probes of Modified Gravity_ , _Phys. Rev. D_ 78 (2008) 043002 [0805.4812].
  * [28] Euclid collaboration, _Euclid: Constraints on f(R) cosmologies from the spectroscopic and photometric primary probes_ , 2306.11053.
  * [29] A. Spurio Mancini and B. Bose, _On the degeneracies between baryons, massive neutrinos and f(R) gravity in Stage IV cosmic shear analyses_ , 2305.06350.
  * [30] G.W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int.J.Theor.Phys._ 10 (1974) 363.
  * [31] E. Bellini and I. Sawicki, _Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity_ , _JCAP_ 07 (2014) 050 [1404.3713].
  * [32] L. Lombriser, C. Dalang, J. Kennedy and A. Taylor, _Inherently stable effective field theory for dark energy and modified gravity_ , _JCAP_ 01 (2019) 041 [1810.05225].
  * [33] N. Frusciante and L. Perenon, _Effective field theory of dark energy: A review_ , _Phys. Rept._ 857 (2020) 1 [1907.03150].
  * [34] A. Zucca, L. Pogosian, A. Silvestri and G.-B. Zhao, _MGCAMB with massive neutrinos and dynamical dark energy_ , _JCAP_ 05 (2019) 001 [1901.05956].
  * [35] S. Srinivasan, D.B. Thomas, F. Pace and R. Battye, _Cosmological gravity on all scales. Part II. Model independent modified gravity N-body simulations_ , _JCAP_ 06 (2021) 016 [2103.05051].
  * [36] S. Srinivasan, D.B. Thomas and R. Battye, _Cosmological gravity on all scales III: non-linear matter power spectrum in phenomenological modified gravity_ , 2306.17240.
  * [37] S. Joudaki et al., _KiDS-450 + 2dFLenS: Cosmological parameter constraints from weak gravitational lensing tomography and overlapping redshift-space galaxy clustering_ , _Mon. Not. Roy. Astron. Soc._ 474 (2018) 4894 [1707.06627].
  * [38] L. Pogosian, M. Raveri, K. Koyama, M. Martinelli, A. Silvestri, G.-B. Zhao et al., _Imprints of cosmological tensions in reconstructed gravity_ , _Nature Astron._ 6 (2022) 1484 [2107.12992].
  * [39] Z. Wang, S.H. Mirpoorian, L. Pogosian, A. Silvestri and G.-B. Zhao, _New MGCAMB tests of gravity with CosmoMC and Cobaya_ , _JCAP_ 08 (2023) 038 [2305.05667].
  * [40] D. Alonso, E. Bellini, P.G. Ferreira and M. Zumalacárregui, _Observational future of cosmological scalar-tensor theories_ , _Phys. Rev. D_ 95 (2017) 063502 [1610.09290].
  * [41] R. Reischke, A. Spurio Mancini, B.M. Schäfer and P.M. Merkel, _Investigating scalar–tensor gravity with statistics of the cosmic large-scale structure_ , _Mon. Not. Roy. Astron. Soc._ 482 (2019) 3274 [1804.02441].
  * [42] A. Spurio Mancini, R. Reischke, V. Pettorino, B.M. Schäfer and M. Zumalacárregui, _Testing (modified) gravity with 3D and tomographic cosmic shear_ , _Mon. Not. Roy. Astron. Soc._ 480 (2018) 3725 [1801.04251].
  * [43] A. Spurio Mancini, F. Köhlinger, B. Joachimi, V. Pettorino, B.M. Schäfer, R. Reischke et al., _KiDS + GAMA: constraints on horndeski gravity from combined large-scale structure probes_ , _Monthly Notices of the Royal Astronomical Society_ 490 (2019) 2155 [http://oup.prod.sis.lan/mnras/article-pdf/490/2/2155/30257339/stz2581.pdf].
  * [44] A. Barreira, B. Li, W.A. Hellwing, C.M. Baugh and S. Pascoli, _Nonlinear structure formation in the Cubic Galileon gravity model_ , _JCAP_ 10 (2013) 027 [1306.3219].
  * [45] B. Bose, M. Tsedrik, J. Kennedy, L. Lombriser, A. Pourtsidou and A. Taylor, _Fast and accurate predictions of the nonlinear matter power spectrum for general models of Dark Energy and Modified Gravity_ , 2210.01094.
  * [46] M. Cataneo, L. Lombriser, C. Heymans, A. Mead, A. Barreira, S. Bose et al., _On the road to percent accuracy: non-linear reaction of the matter power spectrum to dark energy and modified gravity_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) 2121 [1812.05594].
  * [47] P.J. Peebles, _The Large-Scale Structure of the Universe_ , Princeton University Press (11, 1980).
  * [48] E.V. Linder, _Cosmic growth history and expansion history_ , _Phys. Rev. D_ 72 (2005) 043529 [astro-ph/0507263].
  * [49] E.V. Linder and R.N. Cahn, _Parameterized beyond-Einstein growth_ , _Astroparticle Physics_ 28 (2007) 481 [astro-ph/0701317].
  * [50] S. Alam, S. Ho and A. Silvestri, _Testing deviations from $\Lambda$CDM with growth rate measurements from six large-scale structure surveys at $z=$0.06–1_, _Mon. Not. Roy. Astron. Soc._ 456 (2016) 3743 [1509.05034].
  * [51] N.-M. Nguyen, D. Huterer and Y. Wen, _Evidence for Suppression of Structure Growth in the Concordance Cosmological Model_ , _Phys. Rev. Lett._ 131 (2023) 111001 [2302.01331].
  * [52] C. Moretti, M. Tsedrik, P. Carrilho and A. Pourtsidou, _Modified gravity and massive neutrinos: constraints from the full shape analysis of BOSS galaxies and forecasts for Stage IV surveys_ , _JCAP_ 12 (2023) 025 [2306.09275].
  * [53] I. Tereno, E. Semboloni and T. Schrabback, _COSMOS weak-lensing constraints on modified gravity_ , _Astron. Astrophys._ 530 (2011) A68 [1012.5854].
  * [54] DESI collaboration, _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , 1611.00036.
  * [55] LSST Dark Energy Science collaboration, _The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document_ , 1809.01669.
  * [56] W. Hu and I. Sawicki, _Models of f(R) Cosmic Acceleration that Evade Solar-System Tests_ , _Phys. Rev. D_ 76 (2007) 064004 [0705.1158].
  * [57] N. Mirzatuny and E. Pierpaoli, _An Accurate Fitting Function For Scale-dependent Growth Rate in Hu-Sawicki $f(R)$ Gravity_, _JCAP_ 09 (2019) 066 [1906.07683].
  * [58] C. Lamman, E. Tsaprazi, J. Shi, N.N. Šarčević, S. Pyne, E. Legnani et al., _The IA Guide: A Breakdown of Intrinsic Alignment Formalisms_ , 2309.08605.
  * [59] S. Bridle and L. King, _Dark energy constraints from cosmic shear power spectra: impact of intrinsic alignments on photometric redshift requirements_ , _New J. Phys._ 9 (2007) 444 [0705.0166].
  * [60] R. Mandelbaum et al., _The WiggleZ Dark Energy Survey: Direct constraints on blue galaxy intrinsic alignments at intermediate redshifts_ , _Mon. Not. Roy. Astron. Soc._ 410 (2011) 844 [0911.5347].
  * [61] C.D. Leonard, M.M. Rau and R. Mandelbaum, _Photometric redshifts and intrinsic alignments: degeneracies and biases in 3 $\times$2pt analysis_, 2401.06060.
  * [62] B. Giblin, M. Cataneo, B. Moews and C. Heymans, _On the road to per cent accuracy – II. Calibration of the non-linear matter power spectrum for arbitrary cosmologies_ , _Mon. Not. Roy. Astron. Soc._ 490 (2019) 4826 [1906.02742].
  * [63] M. Cataneo, J. Emberson, D. Inman, J. Harnois-Deraps and C. Heymans, _On the road to per cent accuracy – III. Non-linear reaction of the matter power spectrum to massive neutrinos_ , _Mon. Not. Roy. Astron. Soc._ 491 (2020) 3101 [1909.02561].
  * [64] B. Bose, M. Cataneo, T. Tröster, Q. Xia, C. Heymans and L. Lombriser, _On the road to per cent accuracy IV: ReACT – computing the non-linear power spectrum beyond $\Lambda$CDM_, _Mon. Not. Roy. Astron. Soc._ 498 (2020) 4650 [2005.12184].
  * [65] B. Bose, B.S. Wright, M. Cataneo, A. Pourtsidou, C. Giocoli, L. Lombriser et al., _On the road to per cent accuracy – V. The non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback_, _Mon. Not. Roy. Astron. Soc._ 508 (2021) 2479 [2105.12114].
  * [66] P. Carrilho, K. Carrion, B. Bose, A. Pourtsidou, J.C. Hidalgo, L. Lombriser et al., _On the road to per cent accuracy VI: the non-linear power spectrum for interacting dark energy with baryonic feedback and massive neutrinos_ , _Mon. Not. Roy. Astron. Soc._ 512 (2022) 3691 [2111.13598].
  * [67] M. Asgari, A.J. Mead and C. Heymans, _The halo model for cosmology: a pedagogical review_ , _arXiv e-prints_ (2023) arXiv:2303.08752 [2303.08752].
  * [68] A. Mead, J. Peacock, C. Heymans, S. Joudaki and A. Heavens, _An accurate halo model for fitting non-linear cosmological power spectra and baryonic feedback models_ , _Mon. Not. Roy. Astron. Soc._ 454 (2015) 1958 [1505.07833].
  * [69] A. Mead, C. Heymans, L. Lombriser, J. Peacock, O. Steele and H. Winther, _Accurate halo-model matter power spectra with dark energy, massive neutrinos and modified gravitational forces_ , _Mon. Not. Roy. Astron. Soc._ 459 (2016) 1468 [1602.02154].
  * [70] A. Mead, _Spherical collapse, formation hysteresis and the deeply non-linear cosmological power spectrum_ , _Mon. Not. Roy. Astron. Soc._ 464 (2017) 1282 [1606.05345].
  * [71] A. Mead, S. Brieden, T. Tröster and C. Heymans, _HMcode-2020: Improved modelling of non-linear cosmological power spectra with baryonic feedback_ , 2009.01858.
  * [72] VIRGO Consortium collaboration, _Stable clustering, the halo model and nonlinear cosmological power spectra_ , _Mon. Not. Roy. Astron. Soc._ 341 (2003) 1311 [astro-ph/0207664].
  * [73] Euclid collaboration, _Euclid preparation: II. The EuclidEmulator – A tool to compute the cosmology dependence of the nonlinear matter power spectrum_ , _Mon. Not. Roy. Astron. Soc._ 484 (2019) 5509 [1809.04695].
  * [74] Euclid collaboration, _Euclid preparation: IX. EuclidEmulator2 – power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations_ , _Mon. Not. Roy. Astron. Soc._ 505 (2021) 2840 [2010.11288].
  * [75] R.E. Angulo, M. Zennaro, S. Contreras, G. Aricò, M. Pellejero-Ibañez and J. Stücker, _The BACCO simulation project: exploiting the full power of large-scale structure for cosmology_ , _Mon. Not. Roy. Astron. Soc._ 507 (2021) 5869 [2004.06245].
  * [76] B. Bose and K. Koyama, _A Perturbative Approach to the Redshift Space Power Spectrum: Beyond the Standard Model_ , _JCAP_ 1608 (2016) 032 [1606.02520].
  * [77] R.K. Sheth and G. Tormen, _Large scale bias and the peak background split_ , _Mon. Not. Roy. Astron. Soc._ 308 (1999) 119 [astro-ph/9901122].
  * [78] R.K. Sheth and G. Tormen, _An Excursion Set Model of Hierarchical Clustering : Ellipsoidal Collapse and the Moving Barrier_ , _Mon. Not. Roy. Astron. Soc._ 329 (2002) 61 [astro-ph/0105113].
  * [79] J.F. Navarro, C.S. Frenk and S.D.M. White, _A Universal density profile from hierarchical clustering_ , _Astrophys. J._ 490 (1997) 493 [astro-ph/9611107].
  * [80] G.W. Horndeski and A. Silvestri, _50 Years of Horndeski Gravity: Past, Present and Future_ , _Int. J. Theor. Phys._ 63 (2024) 38 [2402.07538].
  * [81] A. Vainshtein, _To the problem of nonvanishing gravitation mass_ , _Phys.Lett._ B39 (1972) 393.
  * [82] G. Dvali, G. Gabadadze and M. Porrati, _4-D gravity on a brane in 5-D Minkowski space_ , _Phys.Lett._ B485 (2000) 208 [hep-th/0005016].
  * [83] S. Peirone, K. Koyama, L. Pogosian, M. Raveri and A. Silvestri, _Large-scale structure phenomenology of viable Horndeski theories_ , _Phys. Rev. D_ 97 (2018) 043519 [1712.00444].
  * [84] K. Koyama and R. Maartens, _Structure formation in the dgp cosmological model_ , _JCAP_ 01 (2006) 016 [astro-ph/0511634].
  * [85] K. Koyama and F.P. Silva, _Non-linear interactions in a cosmological background in the DGP braneworld_ , _Phys. Rev. D_ 75 (2007) 084040 [hep-th/0702169].
  * [86] Y. Wen, N.-M. Nguyen and D. Huterer, _Sweeping Horndeski canvas: new growth-rate parameterization for modified-gravity theories_ , _JCAP_ 09 (2023) 028 [2304.07281].
  * [87] F. Schmidt, W. Hu and M. Lima, _Spherical Collapse and the Halo Model in Braneworld Gravity_ , _Phys. Rev. D_ 81 (2010) 063005 [0911.5178].
  * [88] A.M.C.L. Brun, I.G. McCarthy, J. Schaye and T.J. Ponman, _Towards a realistic population of simulated galaxy groups and clusters_ , _Mon. Not. Roy. Astron. Soc._ 441 (2014) 1270 [1312.5462].
  * [89] I.G. McCarthy, J. Schaye, S. Bird and A.M.C. Le Brun, _The BAHAMAS project: Calibrated hydrodynamical simulations for large-scale structure cosmology_ , _Mon. Not. Roy. Astron. Soc._ 465 (2017) 2936 [1603.02702].
  * [90] V. Springel et al., _First results from the IllustrisTNG simulations: matter and galaxy clustering_ , _Mon. Not. Roy. Astron. Soc._ 475 (2018) 676 [1707.03397].
  * [91] J. Schaye et al., _The FLAMINGO project: cosmological hydrodynamical simulations for large-scale structure and galaxy cluster surveys_ , _Mon. Not. Roy. Astron. Soc._ 526 (2023) 4978 [2306.04024].
  * [92] A. Schneider, R. Teyssier, J. Stadel, N.E. Chisari, A.M.C. Le Brun, A. Amara et al., _Quantifying baryon effects on the matter power spectrum and the weak lensing shear correlation_ , _JCAP_ 03 (2019) 020 [1810.08629].
  * [93] A. Schneider and R. Teyssier, _A new method to quantify the effects of baryons on the matter power spectrum_ , _JCAP_ 12 (2015) 049 [1510.06034].
  * [94] S.K. Giri and A. Schneider, _Emulation of baryonic effects on the matter power spectrum and constraints from galaxy cluster data_ , _JCAP_ 12 (2021) 046 [2108.08863].
  * [95] A. Schneider, S.K. Giri, S. Amodeo and A. Refregier, _Constraining baryonic feedback and cosmology with weak-lensing, X-ray, and kinematic Sunyaev–Zeldovich observations_ , _Mon. Not. Roy. Astron. Soc._ 514 (2022) 3802 [2110.02228].
  * [96] DES collaboration, _Constraining the baryonic feedback with cosmic shear using the DES Year-3 small-scale measurements_ , _Mon. Not. Roy. Astron. Soc._ 518 (2023) 5340 [2206.08591].
  * [97] S. Grandis, G. Arico’, A. Schneider and L. Linke, _Determining the baryon impact on the matter power spectrum with galaxy clusters_ , _Mon. Not. Roy. Astron. Soc._ 528 (2024) 4379 [2309.02920].
  * [98] A. Schneider, S.K. Giri, S. Amodeo and A. Refregier, _Constraining baryonic feedback and cosmology with weak-lensing, X-ray, and kinematic Sunyaev–Zeldovich observations_ , _Mon. Not. Roy. Astron. Soc._ 514 (2022) 3802 [2110.02228].
  * [99] DES collaboration, _Weak lensing combined with the kinetic Sunyaev Zel’dovich effect: A study of baryonic feedback_ , 2404.06098.
  * [100] G. Aricò, R.E. Angulo, C. Hernández-Monteagudo, S. Contreras, M. Zennaro, M. Pellejero-Ibañez et al., _Modelling the large scale structure of the Universe as a function of cosmology and baryonic physics_ , 1911.08471.
  * [101] G. Aricò, R.E. Angulo, S. Contreras, L. Ondaro-Mallea, M. Pellejero-Ibañez and M. Zennaro, _The BACCO simulation project: a baryonification emulator with neural networks_ , _Mon. Not. Roy. Astron. Soc._ 506 (2021) 4070 [2011.15018].
  * [102] E. Semboloni, H. Hoekstra, J. Schaye, M.P. van Daalen and I.J. McCarthy, _Quantifying the effect of baryon physics on weak lensing tomography_ , _Mon. Not. Roy. Astron. Soc._ 417 (2011) 2020 [1105.1075].
  * [103] T. Eifler, E. Krause, S. Dodelson, A. Zentner, A. Hearin and N. Gnedin, _Accounting for baryonic effects in cosmic shear tomography: Determining a minimal set of nuisance parameters using PCA_ , _Mon. Not. Roy. Astron. Soc._ 454 (2015) 2451 [1405.7423].
  * [104] I. Mohammed and N.Y. Gnedin, _Baryonic Effects in Cosmic Shear Tomography: PCA Parameterization and the Importance of Extreme Baryonic Models_ , _Astrophys. J._ 863 (2018) 173 [1707.02332].
  * [105] H.-J. Huang, T. Eifler, R. Mandelbaum and S. Dodelson, _Modelling baryonic physics in future weak lensing surveys_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) 1652 [1809.01146].
  * [106] B. Joachimi and S.L. Bridle, _Simultaneous measurement of cosmology and intrinsic alignments using joint cosmic shear and galaxy number density correlations_ , _Astron. Astrophys._ 523 (2010) A1 [0911.2454].
  * [107] Euclid collaboration, _Euclid: Validation of the MontePython forecasting tools_ , 2303.09451.
  * [108] A. Hall, _Cosmology from weak lensing alone and implications for the Hubble tension_ , _Mon. Not. Roy. Astron. Soc._ 505 (2021) 4935 [2104.12880].
  * [109] J.U. Lange, _nautilus: boosting Bayesian importance nested sampling with deep learning_ , _Mon. Not. Roy. Astron. Soc._ 525 (2023) 3181 [2306.16923].
  * [110] Kilo-Degree Survey, Dark Energy Survey collaboration, _DES Y3 + KiDS-1000: Consistent cosmology combining cosmic shear surveys_ , _Open J. Astrophys._ 6 (2023) 2305.17173 [2305.17173].
  * [111] C. Arnold and B. Li, _Simulating galaxy formation in f(R) modified gravity: Matter, halo, and galaxy-statistics_ , _Mon. Not. Roy. Astron. Soc._ 490 (2019) 2507 [1907.02980].
  * [112] M.P. van Daalen, J. Schaye, C.M. Booth and C.D. Vecchia, _The effects of galaxy formation on the matter power spectrum: A challenge for precision cosmology_ , _Mon. Not. Roy. Astron. Soc._ 415 (2011) 3649 [1104.1174].
  * [113] A. Schneider, N. Stoira, A. Refregier, A.J. Weiss, M. Knabenhans, J. Stadel et al., _Baryonic effects for weak lensing. Part I. Power spectrum and covariance matrix_ , _JCAP_ 04 (2020) 019 [1910.11357].
  * [114] W. Elbers et al., _The FLAMINGO project: the coupling between baryonic feedback and cosmology in light of the $S_{8}$ tension_, 2403.12967.
  * [115] Euclid collaboration, _Euclid: Constraining linearly scale-independent modifications of gravity with the spectroscopic and photometric primary probes_ , 2306.12368.
  * [116] G. Aricò, R.E. Angulo and M. Zennaro, _Accelerating Large-Scale-Structure data analyses by emulating Boltzmann solvers and Lagrangian Perturbation Theory_ , 2104.14568.
  * [117] H. Hildebrandt et al., _KiDS+VIKING-450: Cosmic shear tomography with optical and infrared data_ , _Astron. Astrophys._ 633 (2020) A69 [1812.06076].
  * [118] B. Joachimi et al., _KiDS-1000 methodology: Modelling and inference for joint weak gravitational lensing and spectroscopic galaxy clustering analysis_ , _Astron. Astrophys._ 646 (2021) A129 [2007.01844].

##  Appendix A Connection between Linder gamma and nDGP

For DGP models without an additional dark energy component, we have
$\gamma_{\rm DGP}\approx 11/16$ (e.g., Ref. [49]). Here we will derive the
value for when we include a dark energy contribution. Following the same steps
as in Ref. [49], one can arrive to the following expression for the growth
index by solving the linearised growth equation, taking a matter dominated era
limit of $\Omega_{w}(a)/\Omega_{\rm m}(a)\ll 1$, $\Omega_{w}$ being the dark
energy density fraction, and $(\mu_{\rm L}-1)\ll 1$ (Eq. 22 in Ref. [49]):

$\displaystyle\gamma=\frac{1}{2}+\frac{1}{4\Omega_{w}(a)}\int_{0}^{1}\frac{\mathrm{d}u}{u}u^{5/2}\Omega_{w}(au)-\frac{3}{2\Omega_{w}(a)}\int_{0}^{1}\frac{\mathrm{d}u}{u}[\mu_{\rm
L}(au)-1]u^{5/2}\,,$ (A.1)

with $\Omega_{w}(a)=1-\Omega_{\rm m}(a)$. For early times
$\Omega_{w}(a)\propto a^{-3w}$ and $\mu_{\rm L}-1=A\,\Omega_{w}(a)$, with $A$
being a parameter to be determined. We can integrate the equation for $\gamma$
to get

$\gamma=\frac{3(1-w-A)}{5-6w}\,,$ (A.2)

where for GR $w=-1$ and $A=0$, so one obtains the standard result of
$\gamma_{\rm GR}=6/11$.

Now, the Friedman equation for flat DGP models is given by (“+” for the normal
branch or nDGP, “-” for the self-accelerating branch or sDGP):

$H^{2}\pm\frac{H}{r_{c}}=\frac{8\pi G}{3}\rho\,.$ (A.3)

We have 2 options: assume that a) $\rho=\rho_{\rm m}$, hence the additional
$H/r_{c}$ acts like an effective dark energy component
$\rho_{w}=\mp\frac{3}{8\pi G}\frac{H}{r_{c}}$; or b) $\rho=\rho_{\rm
m}+\rho_{\rm DE}$, hence the effective dark energy has two contributions
$\rho_{w}=\rho_{\rm DE}\mp\frac{3}{8\pi G}\frac{H}{r_{c}}$.

For the first case without dark energy we find from the Friedman equation:
$\Omega_{\rm m}(a)=1\pm 2\sqrt{\Omega_{rc}}/E(a)$ or
$E(a)=\mp\sqrt{\Omega_{rc}}+\sqrt{\Omega_{rc}+\Omega_{\rm m}a^{-3}}$. The
former implies that from $E(a=1)=1$ follows $\mp
2\sqrt{\Omega_{rc}}=1-\Omega_{\rm m}$. Hence, for all realistic $\Omega_{\rm
m}<1$ only the self-accelerating branch is relevant. From
$\dot{\rho}_{w}=-3H(1+w)\rho_{w}$ and a derivative of Eq. A.3 we find
$w=-1/(1+\Omega_{\rm m}(a))$. From Eq. 2.14 and 2.13 we compute
$\beta=-\frac{1+\Omega^{2}_{\rm m}(a)}{1-\Omega^{2}_{\rm m}(a)}$ and $\mu_{\rm
L}-1=-\frac{1}{3}\frac{1-\Omega^{2}_{\rm m}(a)}{1+\Omega^{2}_{\rm
m}(a)}=A\Omega_{w}(a)$ with $A=-\frac{1}{3}\frac{1+\Omega_{\rm
m}(a)}{1+\Omega^{2}_{\rm m}(a)}$. Combining these findings together with the
limit of $a\rightarrow 0$, i.e. $A\rightarrow-1/3$ and $w\rightarrow-1/2$, one
gets $\gamma=11/16$ from Eq. A.2 exactly like in Ref. [49] (recall Eq. A.1 is
derived for early times).

Overall, for small values of $\Omega_{rc}$, $\gamma=11/16$ is a good
approximation at all redshifts for this particular case of the flat DGP or
sDGP model, where the value of the extended parameter is directly related to
the matter density. While for larger values of $\Omega_{rc}$ (hence smaller
values of $\Omega_{\rm m}$) as $a\rightarrow 1$, $A\rightarrow-1/3$ and
$w\rightarrow-1$ (mimicking a cosmological constant today). This decreases the
value of the growth index to $\gamma\rightarrow 7/11$, when computed according
to Eq. A.2. Even the most extreme values of $\Omega_{rc}$ show at most $2\%$
deviation at lower redshifts when compared against solutions from numerically
equating $\mu^{\gamma}_{\rm L}=\mu^{\rm DGP}_{\rm L}$ at each redshift
individually. However, note that in this scenario the growth index value is
always larger than its GR limit. In other words, the structure growth is
suppressed relative to the standard cosmology.

For the second case with the cosmological constant, $\Omega_{rc}$ becomes an
independent parameter. We discuss only the case of only including a
cosmological constant, but the same logic can be applied to any other
parametrisation of the dark energy component. We fix the expansion to the
$\Lambda$CDM cosmology: from Eq. A.3, we have
$E(a)=\mp\sqrt{\Omega_{rc}}+\sqrt{\Omega_{rc}+\Omega_{\rm m}a^{-3}+\Omega_{\rm
DE}(a)}$, and from the equality $E=E_{\Lambda\rm CDM}$ we get $\Omega_{\rm
DE}(a)=(1-\Omega_{\rm m})\pm 2\sqrt{\Omega_{rc}}\sqrt{\Omega_{\rm
m}a^{-3}+(1-\Omega_{\rm m})}$. From Eq. 2.14 and 2.13: $\mu_{\rm
L}-1=\frac{2\sqrt{\Omega_{rc}}}{3}\frac{1}{2\sqrt{\Omega_{rc}}\pm
E(a)[2-\Omega_{\rm m}(a)]}$ with “+” for nDGP and “-” for sDGP. The
corresponding limits are: $\mu_{\rm
L}-1\rightarrow\pm\frac{2\sqrt{\Omega_{rc}}}{3\sqrt{\Omega_{\rm m}}}a^{3/2}$
for early times, while the late-time limit tends to $\mu_{\rm
L}-1\rightarrow\frac{2\sqrt{\Omega_{rc}}}{3}\frac{1}{2\sqrt{\Omega_{rc}}\pm[2-\Omega_{\rm
m}]}$. This means that at high redshift the assumption of $(\mu_{\rm
L}-1)/\Omega_{w}(a)\approx\mathrm{const}$ is not valid in this case as
$\Omega_{w}(a)\propto a^{3}$, while $(\mu_{\rm L}-1)\propto a^{3/2}$. As a
consequence the integration in Eq. A.2 for the early-time limit yields a
diverging value of $\gamma=\gamma_{\rm
GR}\mp\frac{\sqrt{\Omega_{rc}\Omega_{\rm m}}}{4(1-\Omega_{\rm m})a^{3/2}}$,
with “-” for nDGP, “+” for sDGP. We can repeat this calculation for
$a\rightarrow 1$ and get $\gamma=\gamma_{\rm
GR}-\frac{2}{5}\frac{\sqrt{\Omega_{rc}}}{(2\sqrt{\Omega_{rc}}\pm[2-\Omega_{\rm
m}])(1-\Omega_{\rm m})}$ , with “+” for nDGP, “-” for sDGP. Therefore, the
addition of the dark energy component in the nDGP model leads to $\gamma\leq
0.55$ for any $\Omega_{rc}\geq 0$, resulting in enhanced growth at linear
scales and in agreement with our findings. We can insert our values for
$\Omega_{rc}$ and $\Omega_{\rm m}$ for the late-time limit and obtain
$\gamma\approx\gamma_{\rm GR}-0.11\approx 0.44$, which is in excellent
agreement with our findings in Sections 4.2 and 4.5.

##  Appendix B Validation with GR: full posteriors and discussion

Figure 10: Full posterior distribution for validation with GR mock data for a
Stage-IV cosmic shear setup with $\ell_{\rm max}=3000$. Different colours
correspond to the growth index model with screening (green), the growth index
model with the unmodified nonlinear growth (orange), and the model used to
produce the mock data (purple). The dashed lines denote models with an
additional baryonic feedback parameter.

In Fig. 10 several degeneracies are apparent due to the fact that
$C_{\ell}\propto\sigma_{8}^{2}\Omega_{\rm m}$: the negative correlations of
$\sigma_{8}-\Omega_{\rm m}$ and $\gamma-\sigma_{8}$, and the positive
correlation of $\gamma-\Omega_{\rm m}$. We overlay contours from the
$\Lambda$CDM modelling to showcase the strong positive correlation between
$\log{10^{10}~{}A_{s}}$ and $\sigma_{8}$, which is weakened when an additional
parameter controlling the amplitude of the power spectrum on large scales,
$\gamma$, is added. However, in all scenarios $\log{10^{10}~{}A_{s}}$ is
strongly anti-correlated with $h$ and $\omega_{c}$. We also note that the zNLA
parameters $\mathcal{A}^{\rm IA}$ and $\eta^{\rm IA}$ are not strongly
degenerate with any other parameters but are anti-correlated with one another;
their constraints are model-independent. From the full posterior we also see
that the screening scale is not detected, which is to be expect for
$\mu_{L}^{\gamma\sim 0.55}\approx 1$ (see Eq. 2.15). Additionally, the growth
index is anti-correlated with the expansion rate $h$, which is due to the
inclusion of large scales, $10<\ell<100$, and tight informative priors on
$\omega_{b}$ and $n_{s}$ (see discussion in Appendix C).

Figure 11: Varying only $\Omega_{\rm m},~{}S_{8},~{}\gamma$ with re-scaled
covariance and using linear scales. Left panel: Constraints from the first
redshift bin. Time evolution of $S_{8}(z)$ in the standard cosmology (purple
crosses) and $\gamma$-parametrisation ($\gamma$-pseudo, orange crosses). Green
solid line denotes the combination of the lensing and IA kernels (re-scaled
and shifted for visualisation purposes). Black solid line denotes the maximum
of the kernel and same constraints on $S_{8}(z)$ in both models. Middle panel:
marginalised posterior distributions for $\log{10^{10}~{}A_{s}}$ from the
first bin (upper plot) and combined first two bins (lower plot) in the
standard cosmology (purple lines) and $\gamma$-parametrisation (orange lines).
Right panel: Marginalised posterior distribution for $\gamma$ and
$\log{10^{10}~{}A_{s}}$ with the first photometric bin (orange contour) and
combined with the second photometric bin (pink contour). Dashed lines denote
the fiducial values.

From the same full posterior distribution we also notice that the
$\gamma+q_{1}$ (solid green lines) and $\Lambda$CDM (solid purple lines)
models constrain the background cosmological parameters, $\Omega_{\rm
m},~{}h$, to the same level of uncertainty. This is well understood since
$\gamma$ only impacts the amplitude of the power spectrum and not the
background. However, we notice drastically different constraints on
$\sigma_{8}$ and $S_{8}$: the contours are much broader in the $\gamma$ models
than in the $\Lambda$CDM case. This difference arises due to the strong
sensitivity of the linear growth factor to the deviation of the growth index
from its GR value at lower redshifts (see Fig. 1). To demonstrate this we fix
all parameters in the model to their fiducial values and vary $\Omega_{\rm
m}$, $S_{8}$ and $\gamma$ in the first redshift bin with a re-scaled
covariance by a factor of 250 (corresponding to smaller error bars), and using
linear scales only, $\ell_{\rm max}=500$. In the left panel of Fig. 11 we show
that, as expected, both models measure the same value of
$S_{8}(z)=\sigma_{8}(z)\sqrt{\Omega_{\rm m}/0.3}$ around the peak of the first
lensing kernel at $z=0.125$. Therefore, it is the variation in $\gamma$-values
that affects the inferred constraints of $S_{8}$ at redshift $z=0$. For the
first redshift bin, the constraints of the primordial amplitude demonstrate
the same trend (see the upper middle panel of Fig. 11). This changes as soon
as we add additional redshift bins (in the lower middle panel). In Section 4.2
we discuss that the cross-correlated bins provide the most information on the
nonlinear scales as well. Overall, in Fig. 10 with 10 redshift bins the
constraints on the primordial amplitude are identical between both models
since we match the growth at high redshift. Contrary to Stage-III surveys
(see, for example, Ref. [118]), a Stage IV-like setup constrains the
primordial amplitude well due to its wide redshift range and large number of
redshift bins.

##  Appendix C Expansion rate: note on degeneracies and priors

Figure 12: Left panel: the parameter dependence of the shear angular power
spectrum in the 3-3 redshift bin without shape noise. Right panel: change in
the orientation of degeneracies between parameters depending on the scale-cuts
and priors used.

In Fig. 10 we see strong degeneracies between all cosmological parameters and
the expansion rate $h$. Fig. 12 shows the derivatives of $C_{\ell}$ in the
$3-3$ photo-$z$ bin with respect to the parameters listed in the legend. From
this figure, it is clear that $C_{\ell}$ is insensitive to $h$ in the region
of $10^{2}<\ell<10^{3}$. We also note the change in the dependence of
$C_{\ell}$ on $h$ in the $\ell<10^{2}$(prominent negative) and $\ell>10^{3}$
(weak positive) regimes. From Eq. 2.3, we see that the impact of $h$ on
$C_{\ell}$ is coming purely from the matter power spectrum. In the linear
regime, i.e., low $\ell$ and $k$, the matter power spectrum is an approximate
power law with its slope depending on $h$, $n_{s}$ and $\omega_{b}$ (for more
detailed discussion see Ref. [108]). We impose a Planck prior on $n_{s}$ and a
BBN prior on $\omega_{b}$, which breaks this degeneracy. In the right panel of
Fig. 12, we demonstrate the rotation of degeneracy between $\sigma_{8}$ and
$h$, when these informative priors are imposed in the analysis with a GR model
on GR mock data. For $\ell_{\rm max}=700$,
$\mathrm{dln}C_{\ell}/\mathrm{dln}h$ and
$\mathrm{dln}C_{\ell}/\mathrm{dln}\sigma_{8}$ have opposite signs, hence the
orange contour demonstrates a positive correlation. The situation changes to
the negative correlation for $\ell\in[200,2000]$ (the purple contour). When
all scales are combined (the pink contour), the orientation of the degeneracy
still slightly prefers the positive correlation characteristic for lower
$\ell$-values. However, this is no longer the case, when the priors are
omitted (the dashed green contour). Similar arguments are applicable to all
other cosmological parameters and their degeneracies with $h$. Also note that
due to the choice of a diagonal Gaussian covariance, our constraints can be
considered optimistic.

##  Appendix D Test with nDGP: full posteriors and discussion

Figure 13: Full posterior distribution for showcasing the screening impact
with nDGP mock data for a Stage-IV cosmic shear setup. Different colours
correspond to the screened growth index model with $\ell_{\rm max}=3000$
(green) and the pseudo model with the modified linear growth and standard
nonlinear structure formation with $\ell_{\rm max}=1000$ (orange). Solid grey
lines mark the true values of the synthetic data, dashed grey lines mark the
parameter values in the standard cosmology.

In Fig. 13 we show the two-dimensional marginalised posteriors for the full
parameter space in $\gamma+q_{1}$ (green) and $\gamma$-pseudo (orange) models
with with $\ell_{\rm max}=3000$ and $\ell_{\rm max}=1000$, respectively. The
only significantly biased parameter is $\log{10^{10}~{}A_{s}}$, which is
explained by a different time evolution of the linear growth factor in nDGP
and the Linder gamma parametrisation. In Fig. 1 we see that for values of
$\gamma\sim 0.3-0.4$ we obtain an offset of $D_{\gamma}/D_{\rm nDGP}\sim
0.97-0.98$. From the posterior maxima and the input fiducial values
$A_{s}^{\rm mean}/A_{s}^{\rm fid}=1.06$ ($\Delta A_{s}=6\%$). Therefore, lower
(than in nDGP) values of the growth factor in the growth index model are
compensated by a higher value of the primordial amplitude. We note that while
the posterior-maxima of $\omega_{c}$ and $h$ are unbiased, their 2-dimensional
contour demonstrates a $1\sigma$ bias that later propagates to a $1\sigma$
lower posterior maximum for $\Omega_{\rm m}$. We also found this bias when
sampling in $\Omega_{\rm m}$ directly.

We notice that even for $\ell_{\rm max}=1000$, there is a significant bias in
$h$ as well as in $\gamma$. Surprisingly, the pseudo $\gamma$ model finds a
value of $\sigma_{8}$ at its GR value (a lower than expected value compared to
nDGP), which is compensated by high values of $\Omega_{\rm m}$ ($1\sigma$ bias
towards lower values in $\omega_{c}$ and $3\sigma$ bias towards lower values
in $h$).

Figure 14: Full posterior distribution for tests with nDGP mock data for a
Stage-IV cosmic shear setup with $\ell_{\rm max}=3000$. Different colours
correspond to the time-dependent growth index model (green) and the same model
but with the second coefficient fixed to $\gamma_{1}=-0.19$(orange). Grey
lines mark the true values of the synthetic data.

In Fig. 14 we show the two-dimensional marginalised posteriors for the full
parameter space in the $\gamma(z)+q_{1}$ model. Clearly, the bias towards
higher values in $\log{10^{10}~{}A_{s}}$ vanishes due to a more accurate
representation of the growth evolution. However, we notice significantly
weaker constraints on $h$. This is expected from the degeneracy between the
expansion rate and $\gamma$ (in this case $\gamma_{0},\,\gamma_{1}$). We show
that this degeneracy is broken when $\gamma_{1}$ is fixed.

##  Appendix E Other full posteriors

In Figs. 15 and 16 we show the two-dimensional marginalised posteriors for the
full parameter space in the analyses mentioned in the main text.

Figure 15: Full posterior distribution for investigating the degeneracy
between the screening and baryonic parameters on $\gamma+q_{1}$ mock data for
a Stage-IV cosmic shear setup with $\ell_{\rm max}=3000$. Different colours
correspond to the flat broad priors on the baryonic feedback parameter
$\log_{10}{M_{c}}$, (green) and a flat tight prior $\mathcal{U}(13.12,13.52)$
on it (orange). Grey lines mark the true values of the synthetic data. Figure
16: Full posterior distribution for the nDGP model with $\Omega_{rc}=0.25$ for
a Stage-IV cosmic shear setup with $\ell_{\rm max}=3000$. Different colours
correspond to the $\gamma+q_{1}$ model (green), the $\gamma(z)+q_{1}$ model
with $\gamma_{1}=-0.2$ (orange), and the exact nDGP model (purple). Grey lines
mark the true values of the synthetic data.