# Combining gravitational lensing and gravitational redshift to measure the
anisotropic stress with future galaxy surveys

Isaac Tutusaus<EMAIL_ADDRESS>Institut de Recherche en
Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES,
14 Av. Edouard Belin, F-31400 Toulouse, France  Départment de Physique
Théorique and Center for Astroparticle Physics, Université de Genève, Quai E.
Ansermet 24, CH-1211 Genève 4, Switzerland  Daniel Sobral Blanco
<EMAIL_ADDRESS>Départment de Physique Théorique and Center for
Astroparticle Physics, Université de Genève, Quai E. Ansermet 24, CH-1211
Genève 4, Switzerland  Camille Bonvin<EMAIL_ADDRESS>Départment de
Physique Théorique and Center for Astroparticle Physics, Université de Genève,
Quai E. Ansermet 24, CH-1211 Genève 4, Switzerland

###### Abstract

Galaxy surveys provide one of the best ways to constrain the theory of gravity
at cosmological scales. They can be used to constrain the two gravitational
potentials encoding time, $\Psi$, and spatial, $\Phi$, distortions, which are
exactly equal at late time within general relativity. Hence, any small
variation leading to a nonzero anisotropic stress, i.e. a difference between
these potentials, would be an indication for modified gravity. Current
analyses usually consider gravitational lensing and redshift-space distortions
to constrain the anisotropic stress, but these rely on certain assumptions
like the validity of the weak equivalence principle, and a specific time
evolution of the functions encoding deviations from general relativity. In
this work, we propose a reparametrization of the gravitational lensing
observable, together with the use of the relativistic dipole of the
correlation function of galaxies to directly measure the anisotropic stress
with a minimum amount of assumptions. We consider the future Legacy Survey of
Space and Time of the Vera C. Rubin Observatory and the future Square
Kilometer Array, and show that combining gravitational lensing and
gravitational redshift with the proposed approach we will achieve model-
independent constraints on the anisotropic stress at the level of $\sim
20\,\%$.

###### pacs:

Valid PACS appear here

##  I Introduction

Testing the laws of gravity at cosmological scales is one of the main goals of
current and future large-scale structure surveys. In particular, one key test
of gravity consists in comparing the two gravitational potentials encoding
distortions in the geometry of our Universe 111We use the perturbed Friedmann
metric:
$\mathrm{d}s^{2}=a^{2}[-(1+2\Psi)\mathrm{d}\tau^{2}+(1-2\Phi)\mathrm{d}\mathbf{x}^{2}]$,
where $\tau$ denotes conformal time.: namely the time distortion, $\Psi$, and
the spatial distortion, $\Phi$. In the $\Lambda$CDM model, these two
potentials are equal at late times 222Note that massive neutrinos generate a
difference between these two potentials at late time, but the amplitude of
this difference has been shown to be very small: $(\Phi+\Psi)/(2\Psi)-1\simeq
10^{-5}$ Adamek  _et al._ (2017).. Most dark energy models do preserve this
equality, whereas theories of gravity beyond general relativity (GR)
generically induce a difference between $\Psi$ and $\Phi$, called anisotropic
stress Motta  _et al._ (2013); Amendola  _et al._ (2018). Measuring the
anisotropic stress provides therefore a stringent way of testing dark energy
and modified gravity theories Song  _et al._ (2011); Amendola  _et al._
(2013a, 2018). A detection of a nonzero anisotropic stress is often referred
to as a _smoking gun_ for modified gravity. In particular, it would allow us
to distinguish between interactions in the dark sector (that preserve the
equality between $\Phi$ and $\Psi$) and a modification of gravity.

At present, two methods have been used to measure the anisotropic stress. The
first one consists in combining measurements of galaxy peculiar velocities
obtained from redshift-space distortions (RSD), with measurements of the Weyl
potential, $(\Phi+\Psi)/2$, obtained from gravitational lensing. Assuming that
galaxies obey Euler’s equation, RSD measurements can be translated into a
measurement of the time distortion, $\Psi$. Comparing this with the Weyl
potential allows one to measure the anisotropic stress, which, at the
precision of current surveys, is consistent with zero Abbott  _et al._ (2019);
Alam  _et al._ (2021); Abbott  _et al._ (2022). This method is very neat and
powerful but it clearly fails if galaxies do not obey Euler’s equation, which
is the case if dark matter is affected by a fifth force (like for example in
coupled quintessence models Amendola (2000)). As shown in Bonvin and Pogosian
(2022), this breaking of the equivalence principle invalidates the method. It
would indeed lead to a measurement of an effective non-zero anisotropic stress
from RSD and lensing, even if $\Phi=\Psi$.

The second method to constrain the anisotropic stress consists in looking at
the speed of propagation of gravitational waves (GWs). In scalar-tensor
theories of gravity (Horndeski theories Horndeski (1974)), one of the
parameters that generates a non-zero anisotropic stress also governs the speed
of propagation of GWs. Combining the constraints on this speed obtained
through the GW and electromagnetic observations of the binary neutron stars
system GW170817 Abbott  _et al._ (2017), with the consideration that dark
energy perturbations in Horndeski may become unstable in the presence of GWs
Creminelli  _et al._ (2020) (which may further limit the parameter space),
leads to constraints on the anisotropic stress of the order of a few percent
Noller (2020). This elegant method, which combines GWs propagation with large-
scale structure, is however applicable only to a specific class of modified
gravity theories and is therefore not model-independent. As before, it relies
on the validity of the weak equivalence principle in Horndeski theories.
Moreover, it assumes that the speed of GWs is independent of frequency and of
time, something that has not been validated observationally.

In this paper, we propose an alternative method to measure the anisotropic
stress, which has the key advantage of being model-independent. Our method
combines direct measurements of the time distortion, $\Psi$, which can be
detected via the effect of gravitational redshift, with measurements of the
Weyl potential from gravitational lensing, to _directly_ probe the relation
between $\Psi$ and $\Phi$. This method does not rely on a specific model of
gravity, nor on a particular behavior of dark matter: it uses the data in a
completely agnostic way to extract measurable quantities and compare them. In
particular, contrary to measurements from RSD and weak lensing, this method is
also valid if the weak equivalence principle is violated, i.e. if dark matter
obeys a fifth force.

Our analysis follows the method presented in Sobral-Blanco and Bonvin (2021),
where it was shown that combining multipoles of the galaxy power spectrum with
the galaxy-galaxy lensing power spectrum provides a direct model-independent
determination of $\eta\equiv\Phi/\Psi$. Here we apply this method to the
coming generation of surveys to determine the precision with which $\eta$ can
be measured. Note that we use here the correlation function and the angular
power spectrum, instead of the power spectrum, to properly account for wide-
angle effects Bonvin  _et al._ (2014); Tansella  _et al._ (2018). Our
forecasts show that combining measurements of gravitational redshift from the
Square Kilometer Array with measurements of the Weyl potential from the Legacy
Survey of Space and Time of the Vera C. Rubin Observatory (LSST) will allow us
to constrain the anisotropic stress, through $2/(1+\eta)=2/(1+\Phi/\Psi)$,
with a precision of $\sim 20\,\%$.

The rest of the paper is organized as follow: in Sec. II we describe the
modeling considered for both the gravitational lensing and gravitational
redshift observables, as well as how they are combined to constrain the
anisotropic stress. In Sec. III we present the specifications of the future
galaxy surveys considered and the methodology used to forecast the constraints
on the different parameters under study. We then present the main results of
the analysis in Sec. IV for both observables and their combination, and we
conclude in Sec. V.

##  II Modeling

We start by deriving the modeling used to express the predictions for our
observables in a model-independent way, and thus directly test the anisotropic
stress with a minimum amount of assumptions. We first consider the
gravitational lensing observable related to the Weyl potential; then the
gravitational redshift related to $\Psi$; and finally the ratio related to the
anisotropic stress.

###  II.1 Gravitational lensing: $\Phi+\Psi$

The standard method to model the evolution of the two gravitational potentials
consists in relating them to the matter density, allowing for a non-zero
anisotropic stress and a modification to Poisson’s equation. Here we develop a
different method, which, as we will show, allows us to measure directly the
evolution of $\Phi+\Psi$ with redshift.

We start by reviewing how $\Phi+\Psi$ evolves with redshift in $\Lambda$CDM.
We define the transfer function of any field, $F$, as

$\displaystyle F(\mathbf{k},z)=T_{F}(k,z)\Psi_{\rm in}(\mathbf{k})\,,$ (1)

where $\Psi_{\rm in}$ denotes the primordial gravitational potential generated
by inflation. In $\Lambda$CDM the transfer function of $\Phi+\Psi$ can be
related to that of the matter density through

$T_{\Phi+\Psi}(k,z)=2T_{\Phi}(k,z)=-3\Omega_{\rm
m}(z)\left[\frac{\mathcal{H}(z)}{k}\right]^{2}T_{\delta}(k,z)\,,$ (2)

where $\Omega_{\rm m}(z)$ stands for the matter density of the Universe at a
given redshift $z$, $\mathcal{H}$ denotes the Hubble rate in conformal time,
and $\delta$ is the matter density fluctuations
$\delta=\delta\rho/\bar{\rho}$. The matter density transfer function is
usually split into a linear part and a boost factor, $B$, which accounts for
the nonlinear evolution of matter fluctuations at small scales

$T_{\delta}(k,z)=T_{\delta}^{\rm lin}(k,z)\sqrt{B(k,z)}\,.$ (3)

In $\Lambda$CDM, the growth of matter density fluctuations in the linear
regime is scale-independent at late times (once radiation is negligible)
333The only scale-dependence within general relativity (GR) is indeed due to
massive neutrinos and it is very small (see, e.g., Blanchard  _et al._
(2020))., and we can therefore relate the transfer function at redshift $z$ to
a redshift $z_{*}$ through the growth function $D_{1}(z)$

$T_{\delta}^{\rm lin}(k,z)=\frac{D_{1}(z)}{D_{1}(z_{*})}T_{\delta}^{\rm
lin}(k,z_{*})\,.$ (4)

Inserting Eqs. (3) and (4) into (2) we find

$\displaystyle T_{\Phi+\Psi}(k,z)=$
$\displaystyle\frac{\mathcal{H}^{2}(z)\Omega_{\rm
m}(z)D_{1}(z)}{\mathcal{H}^{2}(z_{*})\Omega_{\rm m}(z_{*})D_{1}(z_{*})}$ (5)
$\displaystyle\times\sqrt{B(k,z)}T_{\Phi+\Psi}(k,z_{*})\,.$

Equation (5) tells us that, in $\Lambda$CDM, the evolution of $\Phi+\Psi$ is
governed by $\Omega_{\rm m}(z)D_{1}(z)$, i.e. that $\Phi+\Psi$ follows
directly the evolution of the density.

In the case of modified gravity, however, the two equalities present in Eq.
(2) are generically modified: first $\Phi$ and $\Psi$ can be different from
each other, and second $\Phi$ may not be related to the density, $\delta$, via
Poisson’s equation. Therefore, generically, the evolution of $\Phi+\Psi$ will
differ from the evolution of the density, governed by $D_{1}$. To account for
this, without restricting ourselves to any particular model of gravity, we
replace $\Omega_{\rm m}(z)D_{1}(z)$ in Eq. (5) by an agnostic function
$J(k,z)$. This function encodes the growth of $\Phi+\Psi$, and by treating it
as a _new degree of freedom_ , that can be directly measured from the data and
that is independent from the growth function $D_{1}$, we effectively allow for
any deviations from GR.

We define the redshift $z_{*}$ to be well in the matter era, before the
accelerated expansion of the Universe started. We assume that, at that
redshift, GR is recovered. In other words, we expect structures to grow as in
GR, when the background evolution of the Universe behaves as in GR. We
therefore set $J(k,z_{*})=\Omega_{\rm m}(z_{*})D_{1}(z_{*})=D_{1}(z_{*})$,
since $\Omega_{\rm m}(z_{*})=1$ in the matter era. With this, the transfer
function for $\Phi+\Psi$ evolves as

$\displaystyle T_{\Phi+\Psi}(k,z)=$
$\displaystyle\frac{\mathcal{H}^{2}(z)J(k,z)}{\mathcal{H}^{2}(z_{*})D_{1}(z_{*})}\sqrt{B(k,z)}T_{\Phi+\Psi}(k,z_{*})\,.$
(6)

Combining Eq. (6) with Eq. (3) and (4), and using that at $z_{*}$ GR is
recovered, i.e.

$T_{\delta}^{\rm
lin}(k,z_{*})=-\frac{1}{3}\left[\frac{k}{\mathcal{H}(z_{*})}\right]^{2}T_{\Phi+\Psi}(k,z_{*})\,,$
(7)

we can relate the transfer function of $\Phi+\Psi$ to the density transfer
function at $z_{*}$ via

$\displaystyle T_{\Phi+\Psi}(k,z)$
$\displaystyle=-3\left[\frac{\mathcal{H}(z)}{k}\right]^{2}J(k,z)\sqrt{B(k,z)}\frac{T_{\delta}^{\rm
lin}(k,z_{*})}{D_{1}(z_{*})}\,.$ (8)

Let us now determine how the function $J$ can be measured with gravitational
lensing. Gravitational lensing can be measured either through shear-shear
correlations or through density-shear correlations (the so-called galaxy-
galaxy lensing). In this paper we focus on the latter, since, as we will see,
it allows us to probe directly the function $J$ at the redshift of the lenses.
In addition, we include in our analysis the galaxy-galaxy correlations since
it allows us to break the degeneracy between the galaxy bias and the function
$J$.

The harmonic power spectra of the galaxy-galaxy lensing can be computed as

$\displaystyle C_{\ell}$
${}^{\\!\\!\Delta\kappa}(z_{i},z_{j})=\int\text{d}z\,n_{i}(z)b_{i}(z)\int\text{d}z^{\prime}n_{j}(z^{\prime})C_{\ell}^{\delta\kappa}(z,z^{\prime})$
$\displaystyle=$
$\displaystyle-\frac{A}{\pi}\int\text{d}z\,n_{i}(z)b_{i}(z)\int\text{d}z^{\prime}n_{j}(z^{\prime})\int\frac{\text{d}k}{k}\left(\frac{k}{k_{*}}\right)^{n_{\rm
s}-1}$ $\displaystyle\times
T_{\delta}(k,z)j_{\ell}(k\chi)\frac{1}{\chi^{\prime}}\int_{0}^{\chi^{\prime}}\text{d}\chi^{\prime\prime}\frac{\chi^{\prime}-\chi^{\prime\prime}}{\chi^{\prime\prime}}\ell(\ell+1)$
$\displaystyle\times
T_{\Phi+\Psi}(k,\chi^{\prime\prime})j_{\ell}(k\chi^{\prime\prime})\,.$ (9)

Here $\Delta$ denotes the galaxy over-density, evaluated at the effective
redshift of the tomographic bin $i$

$\displaystyle\Delta(z_{i},\mathbf{n}_{i})=b_{i}\delta(z_{i},\mathbf{n}_{i})+\frac{1}{\mathcal{H}}\partial_{\chi}(\mathbf{V}\cdot\mathbf{n}_{i})\,,$
(10)

where $\mathbf{n}_{i}$ is the direction of the pixel $i$, $\mathbf{V}$ is the
galaxy peculiar velocity, $b_{i}$ is the linear bias, and $\partial_{\chi}$
denotes a derivative with respect to the comoving distance $\chi$. $\kappa$
represents the convergence, evaluated in pixel $j$ with effective redshift
$z_{j}$. It is related to the Weyl potential by

$\displaystyle\kappa(z_{j},\mathbf{n}_{j})$
$\displaystyle=\int_{0}^{\chi_{j}}\\!d\chi\,\frac{\chi_{j}-\chi}{2\chi_{j}\chi}\,\Delta_{\Omega}(\Phi+\Psi)\big{(}\chi,\mathbf{n}_{j}\big{)}\,,$
(11)

where $\chi_{j}\equiv\chi(z_{j})$ and $\Delta_{\Omega}$ is the Laplace
operator on the sphere. In Eq. (9), $n_{i}$ and $n_{j}$ denote the galaxy
distribution function of the lenses and sources, respectively, and $j_{\ell}$
stands for the spherical Bessel function of order $\ell$. The parameters $A$,
$n_{\rm s}$, and $k_{*}$ denote the amplitude, spectral index, and pivot scale
of the primordial power spectrum defined through444Note that $A$ is related to
the amplitude $A_{\rm s}$ defined in Planck through $A=8\pi^{2}A_{\rm s}/9$.
$k^{3}\langle\Psi_{\rm in}(\mathbf{k})\Psi_{\rm
in}(\mathbf{k}^{\prime})\rangle=(2\pi)^{3}A(k/k_{*})^{n_{\rm
s}-1}\delta(\mathbf{k}+\mathbf{k}^{\prime})$. Note that in Eq. (9), we have
neglected the correlation between the convergence and RSD [the second term in
Eq. (10)], as is done, e.g., in Blanchard  _et al._ (2020), since those are
subdominant for thick tomographic bins.

Using the Limber approximation and inserting Eqs. (8), (3), and (4) we can
express the harmonic power spectra as

$\displaystyle C_{\ell}^{\Delta\kappa}(z_{i}$
$\displaystyle,z_{j})=\int\text{d}z\,n_{i}(z)b_{i}(z)\int\text{d}z^{\prime}n_{j}(z^{\prime})\frac{\chi^{\prime}-\chi}{\chi\chi^{\prime}}$
$\displaystyle\times\frac{3\ell(\ell+1)}{2(\ell+1/2)^{2}}\mathcal{H}^{2}(z)\frac{J(z)}{D_{1}(z)}P_{\delta\delta}\left(k_{\ell},\chi\right)\,,$
(12)

where $k_{\ell}\equiv(\ell+1/2)/\chi$. Here we have neglected the
$k$-dependence of the functions $D_{1}$ and $J$. This is a common
approximation in large-scale structure analyses, see, e.g. Alam  _et al._
(2021), which is motivated by the fact that in the quasistatic approximation,
many models of modified gravity lead to a scale-independent growth of
structure Gleyzes  _et al._ (2015, 2016); Raveri  _et al._ (2023). However,
the methodology can be generalized to include the scale-dependence in case of
need. From Eq. (II.1), we see that $J$ is evaluated at the redshift of the
lenses, whose distribution is given by $n_{i}(z)$.

Writing now the density power spectrum as

$P_{\delta\delta}\left(k_{\ell},\chi\right)=\left[\frac{D_{1}(z)}{D_{1}(z_{*})}\right]^{2}P_{\delta\delta}^{\rm
lin}\left(k_{\ell},\chi_{*}\right)B\left(k_{\ell},\chi\right)\,,$ (13)

we can rewrite Eq. (II.1) as

$\displaystyle C_{\ell}^{\Delta\kappa}(z_{i},z_{j})=$
$\displaystyle\frac{3}{2}\int\text{d}z\,n_{i}(z)\mathcal{H}^{2}(z)\hat{b}_{i}(z)\hat{J}(z)$
$\displaystyle\times B\left(k_{\ell},\chi\right)\frac{P_{\delta\delta}^{\rm
lin}\left(k_{\ell},\chi_{*}\right)}{\sigma_{8}^{2}(z_{*})}$
$\displaystyle\times\int\text{d}z^{\prime}n_{j}(z^{\prime})\frac{\chi^{\prime}(z^{\prime})-\chi(z)}{\chi(z)\chi^{\prime}(z^{\prime})}\,,$
(14)

where we have defined

$\hat{J}(z)\equiv\frac{J(z)\sigma_{8}(z)}{D_{1}(z)}=\frac{J(z)\sigma_{8}(z_{*})}{D_{1}(z_{*})}\,,$
(15)

and

$\hat{b}_{i}(z)\equiv b_{i}(z)\sigma_{8}(z)\,.$ (16)

From Eq. (14), we see that galaxy-galaxy lensing is affected by four distinct
ingredients:

  1. 1.

The background expansion of the Universe, through the Hubble function
$\mathcal{H}(z)$ and the comoving distance $\chi(z)$.

  2. 2.

The density fluctuations at redshift $z_{*}$, before acceleration started.

  3. 3.

The nonlinear boost factor $B$.

  4. 4.

The evolution of density and gravitational potential at late times, through
the two functions $\hat{b}$ and $\hat{J}$.

Ingredients 1 and 2 are tightly constrained by cosmic microwave background
(CMB) measurements. In the following we therefore fix them, for simplicity,
using the latest Planck $\Lambda$CDM values for $A_{\rm s}$, $n_{\rm s}$,
$\Omega_{\rm m,0}$, $\Omega_{\rm b,0}$, and $h$ Aghanim  _et al._ (2020)
555Note that here we assume that the background evolution is consistent with
$\Lambda$CDM predictions, since this has been so far confirmed by
observations. This assumption is common in large-scale structure analyses,
see, e.g., Alam  _et al._ (2021).. But we note that the methodology described
in this work would allow for a combination of galaxy survey data with CMB
observations to constrain all parameters together, including the cosmological
and the $\hat{J}$ and $\hat{b}$ parameters. Ingredient 3 depends in principle
on the theory of gravity. However, the standard method to infer the nonlinear
boost is by generating cosmological simulations, which require the choice of a
specific modified gravity model. Therefore, it is not possible to obtain a
general nonlinear boost factor for our parametrization. In this work, we
follow the approach considered in the Dark Enery Survey analysis Abbott  _et
al._ (2022) and keep the standard halofit nonlinear boost Takahashi  _et al._
(2012), while limiting the range of scales used in the analysis to avoid
entering deeply into the nonlinear regime. We also consider a more stringent
scale cut to assess the impact of this choice of boost in our results, showing
that it is subdominant and constitutes, therefore, an acceptable
approximation. Finally, ingredient 4 is what we want to measure in this work.
In practice, we assume $\hat{J}$ and $\hat{b}$ to be free parameters with a
constant value within each tomographic bin, and we focus on constraining them.

As can be seen from Eq. (14), once we consider $\hat{J}$ and $\hat{b}$ as free
parameters with a constant amplitude in each tomographic bin, their product is
fully degenerate. Therefore, as mentioned at the beginning of the section, we
do not consider galaxy-galaxy lensing measurements alone, but rather their
combination with galaxy clustering, using the same photometrically-selected
lenses. This is the so-called 2x2pt analysis.

The harmonic power spectra for galaxy clustering, using the Limber
approximation, is given by

$\displaystyle C_{\ell}^{\Delta\Delta}(z_{i},z_{j})=$
$\displaystyle\int\text{d}z\,n_{i}(z)n_{j}(z){\color[rgb]{0,0,0}\frac{\mathcal{H}(z)(1+z)}{\chi^{2}(z)}}\hat{b}_{i}(z)\hat{b}_{j}(z)$
$\displaystyle\times B\left(k_{\ell},\chi\right)\frac{P_{\delta\delta}^{\rm
lin}\left(k_{\ell},\chi_{*}\right)}{\sigma_{8}^{2}(z_{*})}\,.$ (17)

Combining this expression for the galaxy clustering observable with the
galaxy-galaxy lensing observable in Eq. (14), we can break the degeneracy
between $\hat{J}$ and $\hat{b}$ and constrain both sets of parameters at the
same time.

This procedure allows us to measure the evolution of the gravitational
potentials, $\Phi+\Psi$, in a model-independent way. The only assumption
underlying such an analysis is that at high redshift $z_{*}$, before
acceleration started, we recover GR. Note that this procedure is similar to
the one used in RSD measurements, where the density power spectrum at $z_{*}$
is constrained by CMB measurements, and RSD are used to measure $\hat{b}$ and
$\hat{f}\equiv f\sigma_{8}$ in a model-independent way in each redshift bin,
see, e.g., Alam  _et al._ (2021).

###  II.2 Comparison with the standard $\mu-\Sigma$ parametrization

Before moving to the next observable, let us compare our approach with the
standard parametrization used in weak lensing analyses. Modifications to GR
are usually encoded into two phenomenological functions, $\mu$ and $\eta$,
that modify Poisson’s equation and the relation between the two gravitational
potentials:

$\displaystyle k^{2}\Psi$ $\displaystyle=-4\pi Ga^{2}\mu(z,k)\delta\rho\,,$
(18) $\displaystyle\Phi$ $\displaystyle=\eta(z,k)\Psi\,.$ (19)

The sum of the gravitational potentials can then be written as

$k^{2}(\Phi+\Psi)=-3\mathcal{H}^{2}(z)\Omega_{\rm
m}(z)\Sigma(z)\delta(k,z)\,,$ (20)

where $\Sigma=\mu(1+\eta)/2$. With this the galaxy-galaxy lensing power
spectra and the galaxy clustering power spectra become

$\displaystyle C_{\ell}^{\Delta\kappa}(z_{i},z_{j})=$
$\displaystyle\frac{3}{2}\int\text{d}z\,n_{i}(z)\mathcal{H}^{2}(z)\Omega_{\rm
m}(z)b_{i}(z)\Sigma(z)$ $\displaystyle\times
B\left(k_{\ell},\chi\right)P_{\delta\delta}^{\rm
lin}\left(k_{\ell},\chi\right)$
$\displaystyle\times\int\text{d}z^{\prime}n_{j}(z^{\prime})\frac{\chi^{\prime}(z^{\prime})-\chi(z)}{\chi(z)\chi^{\prime}(z^{\prime})}\,,$
(21) $\displaystyle C_{\ell}^{\Delta\Delta}(z_{i},z_{j})=$
$\displaystyle\int\text{d}z\,n_{i}(z)n_{j}(z){\color[rgb]{0,0,0}\frac{\mathcal{H}(z)(1+z)}{\chi^{2}(z)}}b_{i}(z)b_{j}(z)$
$\displaystyle\times B\left(k_{\ell},\chi\right)P_{\delta\delta}^{\rm
lin}\left(k_{\ell},\chi\right)\,.$ (22)

The galaxy-galaxy lensing power spectra depend directly on the function
$\Sigma(z)$. In addition, both the galaxy-galaxy lensing power spectra and the
clustering power spectra depend on $\mu$, since the evolution of $\delta$, and
consequently the matter power spectrum at redshift $z$, are sensitive to
$\mu$.

Table 1: Comparison of the parameters used in the standard $\mu-\Sigma$ approach and in our approach. $\mu-\Sigma$ approach |  Our approach
---|---
$A_{s},\ n_{s},\ \Omega_{\rm b,0},\ \Omega_{\rm m,0},\ h$ |  $A_{s},\ n_{s},\ \Omega_{\rm b,0},\ \Omega_{\rm m,0},\ h$
$b(z)$ |  $\hat{b}(z)=b(z)\sigma_{8}(z)$
$\Sigma(z)$ |  $\hat{J}(z)=\Sigma(z)\Omega_{\rm m}(z)\sigma_{8}(z)$
$\mu(z)$ |  adding RSD: $\hat{f}(z)=f(z)\sigma_{8}(z)$

Comparing Eq. (21) with Eq. (14) we see that

$\displaystyle\hat{J}(z)=\Sigma(z)\Omega_{\rm m}(z)\sigma_{8}(z)\,.$ (23)

In Table 1, we list the set of free parameters in the standard $\mu-\Sigma$
parametrization and in our parametrization. The parameters in the first line
are the standard cosmological parameters that are best determined by the CMB,
and that are the same in both cases. In our forecasts, we will keep them fixed
for simplicity, but in practice, in both approaches, one can combine lensing
with CMB measurements to constrain these parameters. In addition to these
standard parameters, the $\mu-\Sigma$ parametrization has three free
parameters per redshift bin ($b,\Sigma,\mu$), whereas our parametrization has
two free parameters per redshift bin ($\hat{b},\hat{J}$). One could conclude
that the $\mu-\Sigma$ parametrization is more powerful since it allows us to
measure one more parameter at each redshift. This is however not the case,
since $b$, $\Sigma$, and $\mu$ are strongly degenerated in Eq. (21) and (22).
To break the degeneracy it is necessary to add new information, through RSD,
which are sensitive to $\mu$ and $b$. In our approach, adding RSD would add
one new free function $\hat{f}(z)=f(z)\sigma_{8}(z)$, that can directly be
measured from RSD and leads to the same number of free parameters in both
approaches.

From this we see that the first key property of our approach is that it
separates clearly the information that can be measured with the 2x2pt lensing
measurements, from the information that can be measured from RSD. The 2x2pt
data measure the evolution of $\Phi+\Psi$ and the evolution of the galaxy
density, whereas RSD measure the evolution of the velocity (and again the
evolution of the galaxy density). In the standard $\mu-\Sigma$ parametrization
this separation cannot be applied since the parameters $b,\Sigma$, and $\mu$
can only be measured by combining lensing with RSD.

The second specificity of our approach is to allow for a direct measurement of
the parameters $\hat{J}$ and $\hat{b}$ in each redshift bin. In contrast, in
the $\mu-\Sigma$ approach, it is not straightforward to measure $\mu(z)$ in
each redshift bin. The reason for this is that $\mu$ enters in Eqs. (21) and
(22) through its impact on the matter power spectrum, $P_{\delta\delta}(k,z)$,
which depends on the whole time evolution of $\mu(z)$. In other words, to
constrain $\mu(z)$ from the 2x2pt data and RSD, we need to solve the evolution
equation for the matter density $\delta$, which reads

$\displaystyle\delta^{\prime\prime}(k,a)+\left(1+\frac{\mathcal{H}^{\prime}(a)}{\mathcal{H}(a)}\right)\delta^{\prime}(k,a)$
$\displaystyle-\frac{3}{2}\frac{\Omega_{m,0}}{a}\left(\frac{\mathcal{H}_{0}}{\mathcal{H}(a)}\right)^{2}\mu(a)\,\delta(k,a)=0\,,$
(24)

where a prime denotes derivatives with respect to $\ln a$. This means that to
measure $\mu$ in the bin $z_{i}$, it is not enough to measure clustering at
$z_{i}$ and weak lensing at $(z_{i},z_{j})$. We need instead measurements of
these quantities in a variety of redshifts larger than $z_{i}$. Two methods
have been used to account for this fact. The simplest way is to assume a given
evolution for $\mu$ with $a$, for example, see Abbott  _et al._ (2019); Alam
_et al._ (2021)

$\displaystyle\mu(a)=1+\mu_{0}\frac{\Omega_{\Lambda}(a)}{\Omega_{\Lambda,0}}\,,$
(25)

where only $\mu_{0}$ is a free parameter, that is constrained from the data.
However, if $\mu$ does not evolve in this way, the constraints on $\mu_{0}$
are not valid. Another possibility is to parametrize $\mu$ and $\Sigma$ in
terms of their values in a number of redshift nodes, and then interpolate
between the nodes to obtain continuous functions, that can be used to solve
Eq. (24), see Raveri  _et al._ (2023). This method does not assume any time
evolution, but it depends on the chosen interpolation method, which introduces
arbitrary correlations between the nodes. For example, using a cubic spline
tends to suppress sharp changes in these functions Raveri  _et al._ (2023). In
contrast, in our approach, since no evolution equation needs to be solved, no
interpolation is needed: $\hat{J}$ can be measured directly in each redshift
bin. A theory prior on the evolution of $\hat{J}$ with redshift can be
introduced if we want, for example, to reduce the number of free parameters,
but this is not required.

Finally, our approach has the advantage to be fully model-independent: we
directly measure $\hat{J}$, i.e. the evolution of $\Phi+\Psi$, without any
assumption on the theory of gravity or on the behavior of dark matter. This is
not the case for the standard $\mu-\Sigma$ approach, which relies on the
validity of Eq. (24) to find $\delta$ for a given $\mu$. As shown in Castello
_et al._ (2022), if dark matter couples differently than baryonic matter to
gravity, or if dark matter is affected by a fifth (nongravitational) force
Bonvin and Pogosian (2022), this equation is modified and the constraints on
$\mu$ and $\Sigma$ are not valid.

Measuring $\hat{J}$ is therefore much more direct and robust than measuring
$\Sigma$ and $\mu$. It is actually completely equivalent to measuring
$\hat{f}=f\sigma_{8}$ from RSD. These measurements can indeed also be done in
a completely model-independent way, in each of the redshift bins of the
survey. In contrast, going from $f\sigma_{8}$ to $\mu$ requires to solve the
evolution equation (24) for $\delta$. In the next section, we will present an
example where measuring $\hat{J}$ is very useful and allows us to measure
$\eta$ redshift bin by reshift bin, without having to assume anything on the
evolution of $\mu$, and, as stated previously, without having to assume a
specific behavior for dark matter. Note that from Table 1, we expect the
relative constraints on $\hat{J}$ and $\hat{f}$ to be of similar amplitude as
the relative constraints on $\Sigma$ and $\mu$.

To finish this section of comparison with other parametrizations, it is worth
considering the recent analysis presented in Camera (2022). The author
considers a template-fitting approach to constrain the growth of matter
perturbations with cosmic shear analyses in a model-independent way, which is
similar to our goal. There is however an important difference with respect to
our approach, which is that in Camera (2022), GR is assumed, and therefore Eq.
(2) is used to relate the evolution of $\Phi+\Psi$ to the evolution of the
density. The main goal of this method is therefore to measure the growth of
density, without being affected by the bias. Hence, the free function that is
fitted from the data is $\Omega_{\rm m}(z)\sigma_{8}(z)$. In our case, on the
other hand, we use lensing to measure directly the evolution of $\Phi+\Psi$,
i.e. $\hat{J}(z)$, without assuming GR. Note that because of that, we choose
our reference power spectrum at $z=z_{*}$ (where we assume GR to be recovered)
instead of $z=0$, as is done in Camera (2022). Furthermore, we assume the
shape of the template power spectrum to be given by CMB measurements, instead
of assuming a fiducial cosmology and accounting for it with an Alcock-
Paczyński parameter. We also differ in the fact that we use mildly nonlinear
scales with a GR nonlinear boost factor. Finally, another relevant difference
is that in Camera (2022) the author considers the cosmic shear observable and,
because of this, he includes the BNT nulling technique proposed in Bernardeau
_et al._ (2014). Such technique allows the author to obtain localized cosmic
shear kernels and therefore constrain the growth at different redshift bins.
In our analysis, we consider the galaxy-galaxy lensing observable. Therefore,
the agnostic $\hat{J}$ function only appears at the level of the lenses, which
are already localized.

Other recent model-independent analyses to test modified gravity that are
worth mentioning are the Dark Energy Survey Year 1 analysis splitting between
growth and geometry Muir  _et al._ (2021) and the Dark Energy Survey Year 3
analysis binning $\sigma_{8}$ as a function of redshift Abbott  _et al._
(2022). The former considered a split of a subset of cosmological parameters,
such that one parameter was sensitive to the growth of perturbations and the
other one was sensitive to the geometry of the Universe. Within $\Lambda$CDM,
the two parameters should agree and provide the same value. In the Year 3
analysis, a successor of this method was considered by introducing a set of
amplitudes (one per redshift bin) that scale the linear power spectrum. These
translate into the value of $\sigma_{8}$ at redshift 0 based on the amplitude
of structure in a given redshift bin. These methods can be seen as consistency
tests of $\Lambda$CDM, where any deviation would consist in an indication for
beyond-$\Lambda$CDM physics. Instead, our method focuses on directly measuring
the evolution of the gravitational potentials and the anisotropic stress as a
function of redshift in a model-independent way.

###  II.3 Gravitational redshift: $\Psi$

The method to measure directly the time distortion $\Psi$ has been presented
in detail in Sobral-Blanco and Bonvin (2022). Here we simply summarize the
main points. Following the same steps as in Sec. II.1, we introduce a new
function $I(k,z)$ to encode the evolution of the gravitational potential
$\Psi$ such that

$\displaystyle
T_{\Psi}(k,z)=\frac{\mathcal{H}^{2}(z)I(k,z)}{\mathcal{H}^{2}(z_{*})D_{1}(z_{*})}{\color[rgb]{0,0,0}\sqrt{B(k,z)}}T_{\Psi}(k,z_{*})\,,$
(26)

where $B$ is the boost defined in Eq. (3). In our forecasts, we restrict the
correlation function to separations larger than $d_{\rm min}$, chosen such
that nonlinearities are negligible, i.e., that the boost plays no role, see
Sec. III.2.

As shown in Sobral-Blanco and Bonvin (2022), the evolution of $\Psi$ can be
measured by cross-correlating the over-density, $\Delta$, of two populations
of galaxies, e.g., a bright and a faint population. $\Psi$ contributes to the
galaxy over-density through the effect of gravitational redshift, which shifts
to the red the spectrum of galaxies situated in a gravitational potential
well. This effect adds to the density and RSD in Eq. (10) and generates a
contribution of the form $\partial_{r}\Psi/\mathcal{H}$ in $\Delta$ Bonvin
(2014). This term has the particularity to produce asymmetries in the
distribution of galaxies Bonvin  _et al._ (2014). Hence it was proposed to
measure it by fitting for a dipole in the cross-correlation of bright and
faint galaxies McDonald (2009); Croft (2013); Bonvin  _et al._ (2014).
However, since Doppler effects also contribute to such a dipole, one needs a
method to disentangle the two types of contributions. More precisely, the
dipole is sensitive to the following contributions in $\Delta$ (that are often
called relativistic effects) 666Note that $\Delta^{\rm rel}$ contains other
relativistic contributions, like Shapiro time-delay, integrated Sachs-Wolfe,
and gravitational lensing Bonvin and Durrer (2011); Yoo  _et al._ (2009);
Challinor and Lewis (2011), but these contributions are strongly subdominant
and negligible at the scales and redshifts used in our analysis Jelic-Cizmek
_et al._ (2021); Lepori  _et al._ (2022).

$\displaystyle\Delta^{\rm rel}(z,\mathbf{n})=$
$\displaystyle\frac{1}{\mathcal{H}}\partial_{r}\Psi+\frac{1}{\mathcal{H}}{\dot{\mathbf{V}}}\cdot\mathbf{n}$
(27)
$\displaystyle+\left(1-5s+\frac{5s-2}{\mathcal{H}r}-\frac{{\dot{\mathcal{H}}}}{\mathcal{H}^{2}}+f^{\rm
evol}\right)\mathbf{V}\cdot\mathbf{n}\,,$

where a dot denotes derivative with respect to conformal time, $s$ is the
magnification bias and $f^{\rm evol}$ is the evolution bias.

In Sobral-Blanco and Bonvin (2022), it was shown that by combining
measurements of the monopole, quadrupole, and hexadecapole of the correlation
function, with a measurement of the dipole, one can measure separately the
evolution of $\Psi$ and the evolution of the velocity, in a model-independent
way. More precisely it was shown that the quantity $\hat{I}$ given by

$\displaystyle\hat{I}(z)\equiv\frac{I(z)\sigma_{8}(z)}{D_{1}(z)}=\frac{I(z)\sigma_{8}(z_{*})}{D_{1}(z_{*})}$
(28)

can be measured directly, redshift bin by redshift bin. Note that as before,
we neglect here the $k$-dependence of $\hat{I}$, since it is a good
approximation in the quasi-static limit.

The parameter $\hat{I}$ is a new parameter, that adds to the ones defined in
Table 1. If galaxies obey Euler’s equation, which is the case if the weak
equivalence principle for dark matter and baryonic matter is valid, then
$\hat{I}$ can be related to $\hat{f}$ through

$\displaystyle\hat{I}=\frac{2\hat{f}}{3}\left[\frac{\mathcal{H}^{\prime}}{\mathcal{H}}+\frac{\hat{f}^{\prime}}{\hat{f}}+1\right]\,.$
(29)

It is therefore possible to reconstruct the evolution of $\Psi$ from RSD, as
has been shown in Motta  _et al._ (2013). On the other hand, if Euler’s
equation is not valid, for example if dark matter is sensitive to a fifth
force, or if baryons and dark matter are not coupled in the same way to
gravity, then Eq. (29) is not valid Bonvin and Pogosian (2022). In this case,
$\hat{I}$ has to be considered as an independent function, that cannot be
inferred from RSD measurements. The dipole is therefore an important new
observable, since it will allow us to measure $\hat{I}$ directly, without
having to assume anything on the behavior of dark matter. With current
surveys, the dipole is unfortunately not detectable Gaztanaga  _et al._
(2017). Forecasts show however that the coming generation of surveys, like
DESI and SKA2, will be able to measure it robustly Bonvin and Fleury (2018);
Beutler and Di Dio (2020); Saga  _et al._ (2022).

###  II.4 Anisotropic stress

Now that we have derived the parametrization of the gravitational lensing and
gravitational redshifit observables, we can combine them to measure
$\eta=\Phi/\Psi$. We have

$\frac{\hat{I}(z)}{\hat{J}(z)}=\frac{T_{\Phi+\Psi}(k,z_{*})}{T_{\Psi}(k,z_{*})}\frac{T_{\Psi}(k,z)}{T_{\Phi+\Psi}(k,z)}=\frac{2\Psi}{\Phi+\Psi}=\frac{2}{1+\eta}\,,$
(30)

where in the first equality we use Eqs. (6) and (26). The boost, $B(k,z)$,
cancels in the ratio since it affects in the same way the evolution of $\Psi$
and of $\Phi+\Psi$ (it encodes indeed the nonlinear evolution of matter
density – see Eq. (3)). Moreover, in the second equality we use that
$(\Phi+\Psi)(z_{*})=2\Psi(z_{*})$, since GR is recovered at $z_{*}$.

We note that we consider the ratio of $\hat{I}$ over $\hat{J}$, or
equivalently the ratio of $\Psi$ over the Weyl potential, instead of the
inverse, because our constraints on $\Psi$ are weaker, as we will see in the
following sections. The fact of having weaker constraints allows $\Psi$ to
become compatible with a null value, which would introduce numerical
instabilities in the inverse of Eq. (30).

Equation (30) is a novel estimator of the anisotropic stress, $\eta$, which is
fully model-independent. It is directly built from measurements of the
functions $\hat{J}$ and $\hat{I}$ in the bins of the surveys. If the ratio
$\hat{I}/\hat{J}$ differs from 1, we can then unambiguously conclude that
gravity is modified. In contrast, the standard $\mu-\Sigma$ approach does not
allow us to measure $\eta$ in a model-independent way. In particular, with the
standard approach, we could detect an apparent deviation from GR: $\eta\neq
1$, even if gravity is not modified and $\Psi=\Phi$. As has been discussed in
Bonvin and Pogosian (2022), this is due to the fact that if dark matter is
affected by a fifth force, RSD do not provide a measurement of the true $\mu$.
As a consequence, the observed $\eta$, which is inferred from a measurement of
$\Sigma$ and $\mu$, will not be the true $\eta$. Our method is therefore
crucial if we want to regard $\eta\neq 1$ as a smoking gun for modified
gravity.

##  III Methodology

In this section we present the methodology used to forecast the constraints on
the observables presented in Sec. II. We first describe the galaxy surveys
considered and their settings, and then show the Fisher matrix forecast used
in this work.

Figure 1: Number of galaxies as a function of redshift for the different
samples considered. The solid lines stand for the LSST lenses, the dashed
lines represent the LSST sources, and the bars represent the SKA spectroscopic
top-hat bins, that have been adapted to the LSST lenses. The effective
redshift for each tomographic bin is represented with a vertical dotted line.

###  III.1 Gravitational lensing with LSST

In the following years we will have access to a huge amount of exquisite
photometric data from the future Stage-IV galaxy surveys. Some examples are
the Vera C. Rubin Observatory that will carry out the Legacy Survey of Space
and Time (LSST 777https://www.lsst.org, Ivezić  _et al._ (2019)), the Euclid
satellite 888https://www.euclid-ec.org Laureijs  _et al._ (2011), or the Nancy
G. Roman Space Telescope 999https://roman.gsfc.nasa.gov Spergel  _et al._
(2015). In this work we focus on LSST, whose main science objective from a
cosmological perspective is to probe dark energy and dark matter. In order to
achieve this, the baseline survey will cover 18 000 square degrees of the
Southern sky during 10 years to obtain accurate photometry in multiple optical
bands. Such a survey will allow us to obtain accurate photometric redshifts
and weak lensing shear measurements for about 27 galaxies per arcmin2. This
will enable precise photometric galaxy clustering, cosmic shear, and galaxy-
galaxy lensing cosmological analyses, also called 3x2pt analyses.

In more detail, we consider the survey settings provided with the public
version of CosmoSIS Zuntz  _et al._ (2015). These consist in five
equipopulated tomographic bins for the sources and ten equipopulated
tomographic bins for the lenses with a total of 27 galaxies per arcmin2, for
both sources and lenses. We represent the galaxy distributions in Fig. 1.
Additionally, we consider a linear galaxy bias model with a constant fiducial
set to $b=2$, but treated as a nuisance parameter in each tomographic bin and
marginalized over 101010We note that a realistic galaxy bias might have some
redshift dependence. However, since we marginalize over its value at each
redshift bin, we have verified that a different fiducial value for the galaxy
bias does not change the size of the final constraints. We therefore keep
$b=2$ as fiducial, for simplicity.. We also consider the nonlinear alignment
model for intrinsic alignments Bridle and King (2007); Kirk  _et al._ (2012)
with a fiducial amplitude set to ${\color[rgb]{0,0,0}A_{\rm IA}=1}$ but
allowed to vary. We note that we limit our analysis to multipoles between
$\ell=20$ and $\ell=2627$ in the optimistic scenario and between $\ell=20$ and
$\ell=750$ in the pessimistic scenario. In both cases we do not go too much
into the nonlinear regime, and therefore the linear modeling for galaxy bias
and intrinsic alignment is still a good approximation for our forecasts.
Finally, we consider an ellipticity total dispersion of
$\sigma_{\epsilon}=0.3$.

It is important to mention that we assume GR is valid at small scales where
intrinsic alignments are important, that is, we include the intrinsic
alignment contribution using a standard GR modeling, even when constraining
$\hat{J}$. The main reason for this choice is the lack of intrinsic alignment
models for modified gravity theories. However, we marginalize over the
amplitude of this effect to account for its impact. We neglect some
observational systematic uncertainties, like a shear calibration bias, or
biases in the mean of our galaxy distributions. Moreover, we neglect some
contributions to the signal like the impact of magnification or RSD for
photometric galaxy clustering and galaxy-galaxy lensing observables. This is
also the approach followed in Blanchard  _et al._ (2020), for example. Note
that we do not expect magnification to alter the measurement of $\hat{b}$ and
$\hat{J}$ at small redshift, where its contribution to the signal is strongly
subdominant. Since, as we will see, $\eta$ is only well measured at small
redshift with our estimator, our analysis should be independent of this
contamination. Also RSD are partially washed out in thick tomographic bins and
therefore a subdominant contribution to the total signal. In any case, this
analysis focuses on determining whether a model-independent measurement of the
anisotropic stress can be performed with future observations. Therefore, our
first aim is to check its feasibility in an optimistic scenario where the main
systematic uncertainties are under control. Obviously this analysis will need
to account for observational uncertainties once real observations are
available and to obtain more precise predictions.

###  III.2 Gravitational redshift with SKA

Future Stage-IV galaxy surveys will also provide us with precise spectroscopic
data probing the galaxy clustering in our Universe. An example of such a
survey is the phase 2 of the future Square Kilometer Array (SKA)
111111https://www.skatelescope.org, with which we will be able to observe
close to a billion of galaxies between redshifts $z=0.1$ and $z=2.0$. It is
expected that this new data will drastically improve our measurements of the
growth history, achieving subpercent measurements of $f(z)\sigma_{8}(z)$ Bull
(2016) using the even multipoles of the galaxy clustering correlation
function. Furthermore, by combining the measurement of the even multipoles
with the dipole we will be able to measure the gravitational redshift with a
precision of $10-30\%$ at late times (see Sobral-Blanco and Bonvin (2022)). In
this work, we use this measurement to build an estimator for the anisotropic
stress by comparing the functions $\hat{I}(z)$ and $\hat{J}(z)$ of Eq. (28)
and Eq. (15), respectively. Importantly, we adapt the redshift bins considered
in Sobral-Blanco and Bonvin (2022) for $\hat{I}(z)$, in order to obtain
constraints at the same effective redshifts than those for $\hat{J}(z)$. This
can easily be done for SKA2, thanks to the very precise determination of
redshift. The new top-hat galaxy distributions and the effective redshifts are
shown in Fig. 1.

We consider the number density and volume specifications presented in Bull
(2016) and follow the approach used in Sobral-Blanco and Bonvin (2022) to
split the populations of galaxies between faint and bright such that we have
the same number of each luminosity type per redshift bin. For the fiducial of
the galaxy biases we use the exponential fitting functions

$\displaystyle b_{\rm B}$ $\displaystyle=ce^{dz}+\frac{\Delta b}{2}\,,$ (31)
$\displaystyle b_{\rm F}$ $\displaystyle=ce^{dz}-\frac{\Delta b}{2}\,,$ (32)

for the bright and faint populations, respectively, where $c=0.554$ and
$d=0.783$, following Bull (2016). As in Sobral-Blanco and Bonvin (2022), we
assume a difference between the two galaxy biases of $\Delta b=1$. This is
consistent with the $\mathcal{O}(1)$ difference that has been measured for
BOSS in Gaztanaga  _et al._ (2017). In our forecasts, we then treat $b_{\rm
B}$ and $b_{\rm F}$ as free parameters in each redshift bin and marginalize
over them. For the magnification bias, we use the model developed in Castello
_et al._ (2022), and we neglect the evolution bias for the two populations.
Once data will be available, both the evolution bias and the magnification
bias of the two populations will be directly measurable from the average
galaxy distribution. Therefore we keep these values fixed in our forecasts.

Finally, we shall perform two forecasts by choosing different values for the
minimum separations between galaxies. This $d_{\rm min}$ must be set at the
scale on which the linear regime is a good approximation. Following Bonvin
_et al._ (2020), for an optimistic scenario, we choose $d_{\rm min}=20\,{\rm
Mpc}/h$, while for the pessimistic case we set $d_{\rm min}=32\,{\rm Mpc}/h$.
In this regime, the boost in Eq. (26) plays a negligible role and can
therefore be set to 1 in our forecasts. In addition, we set the maximum
separation such that it is consistent with the size of our smallest redshift
bin, $d_{\rm max}=120\,{\rm Mpc}/h$, in all of our bins.

###  III.3 Fisher matrix forecasts

Once we have our theoretical predictions for the observables derived in Sec.
II with the specifications provided above, we can forecast the uncertainties
on the different parameters using a Fisher matrix formalism. We summarize the
parameters considered in our analyis, together with their fiducial values, in
Tables 2 and 3 in the Appendix, but for simplicity our vector of free
parameters is given by

$\theta=\\{\hat{J}_{i},\hat{b}_{i},A_{\rm
IA},\hat{I}_{j},\hat{f}_{j},\hat{b}_{\text{B},j},\hat{b}_{\text{F},j}\\}\,,$
(33)

where $i$ runs over the ten tomographic bins for LSST and $j$ runs over the
seven bins for SKA2.

We recall that the Fisher matrix is defined as the expectation value of the
second derivative of the logarithm of the likelihood with respect to the
parameters of the model:

$F_{\alpha\beta}=\Braket{-\frac{\partial^{2}\,\text{ln}(\mathcal{L})}{\partial\theta_{\alpha}\partial\theta_{\beta}}}\,.$
(34)

For a Gaussian likelihood, and neglecting any dependence of the covariance of
the observables on the model parameters, the Fisher matrix can be expressed as

$F_{\alpha\beta}=\sum_{pq}\frac{\partial\mu_{p}}{\partial\theta_{\alpha}}\left(\text{C}^{-1}\right)_{pq}\frac{\partial\mu_{q}}{\partial\theta_{\beta}}\,,$
(35)

where $\mu$ is the mean of the data vector and C is the covariance matrix of
the data. These will correspond to the galaxy-galaxy lensing and galaxy
clustering harmonic spectra for LSST, and to the correlation function
multipoles for SKA2. We consider a Gaussian covariance for the former, meaning
that we account for the cosmic variance and shape/shot noise, but neglect non-
Gaussian terms like the supersample covariance, as was done in Blanchard  _et
al._ (2020). For SKA2, we include shot noise and cosmic variance in the
variance of the multipoles (see Appendix C of Bonvin and Fleury (2018)) and
account for cross-correlations between the different combinations of
luminosity pairs. We note that, because of our $\Lambda$CDM fiducial, we use
GR to compute these covariances.

Once the Fisher matrix is computed, we estimate the covariance matrix of the
model parameters as the inverse of the Fisher matrix:

$C_{\alpha\beta}=\left(\text{F}^{-1}\right)_{\alpha\beta}\,.$ (36)

In practice, we use CosmoSIS to call the CAMB Boltzmann solver Lewis  _et al._
(2000); Howlett  _et al._ (2012) and compute the Fisher matrix for the
$\hat{J}$ parameters. We then use the same Boltzmann solver to build the
Fisher matrix for the $\hat{I}$ parameters from the correlation function
multipoles. We note that our gravitational lensing observable is also
sensitive to $\hat{b}$ and the amplitude of intrinsic alignments, while our
gravitational redshift observable depends on the generic growth function
$\hat{f}=f\sigma_{8}$ 121212Note that we treat $\hat{f}$ as a free function,
with unknown time evolution. In particular, we do not assume that
$f=\frac{d\ln D_{1}}{d\ln a}$, since this is only true if the continuity
equation for dark matter is valid. Since we want to remain agnostic about the
behavior of dark matter, we do not assume that this equation is valid. and
galaxy biases $\hat{b}_{\rm B}=b_{\rm B}\sigma_{8}$ and $\hat{b}_{\rm
F}=b_{\rm F}\sigma_{8}$. These additional parameters are considered nuisance
parameters and we marginalize over them when providing constraints on
$\hat{J}$ or $\hat{I}$. We also remind the reader that these observables
depend on the cosmological parameters providing the spectrum of matter
perturbations at $z=z_{*}$ but we consider the cosmology fixed.

In addition to the forecast constraints on $\hat{J}$ and $\hat{I}$, sensitive
to $\Phi+\Psi$ and $\Psi$, respectively, we want to combine them to constrain
the anisotropic stress. Instead of building a Jacobian transformation to move
from one set of parameters to another one, and to keep the nonlinearities that
may arise in the transformation, we generate synthetic chains from the
individual Fisher matrices. In more detail, given a Fisher matrix for the
$\hat{J}$ parameters and a Fisher matrix for the $\hat{I}$ parameters, we
invert them to obtain the covariances of the parameters. Then, for each one of
them, we generate a mock chain centered at our fiducial and with random points
drawn from a multidimensional Gaussian distribution with the corresponding
covariance. We finally generate a third chain from the ratio of the other two
at each point, which provides the posterior on $\hat{I}/\hat{J}=2/(1+\eta)$.
It is important to mention that we only consider seven redshift bins for
$\hat{I}$, while we use ten redshift bins for $\hat{J}$, as can be seen in
Fig. 1. The main reason for this choice is the lack of constraining power on
$\hat{I}$ at higher redshifts. Therefore, we further marginalize over the last
three tomographic bins for the gravitational lensing observable when it is
combined with the gravitational redshift.

##  IV Results

In this section we present the main results of the analysis. We first focus on
the gravitational lensing observable and the constraints on $\Phi+\Psi$. We
then present the results on $\Psi$ coming from the dipole of the correlation
function. Finally, we derive constraints on the anisotropic stress from the
combination of the two observables.

###  IV.1 Gravitational lensing

Figure 2: Forecast $1\,\sigma$ uncertainties for $\hat{J}$ (top panel) and
$\hat{I}$ (bottom panel) with respect to their fiducial value for each
tomographic bin. The red (light gray) error bars correspond to the optimistic
scenario, while the blue (dark gray) error bars represent the pessimistic
settings described in the text.

In the top panel of Fig. 2 we present the $1\,\sigma$ forecast uncertainties
for $\hat{J}$ with respect to their fiducial value as a function of redshift.
We show the optimistic scenario in red, corresponding to scale cuts $\ell_{\rm
max}=2627$, and the pessimistic scenario in blue, corresponding to $\ell_{\rm
max}=750$. A given offset has been added in the x axis for illustrative
purposes.

As can be seen in Fig. 2, there is a degradation of the constraints on
$\hat{J}$ as a function of redshift. This behavior is essentially due to two
effects. First, at high redshift, the tomographic bins are wider, which
implies that there is a more significant smoothing of the galaxy clustering
distribution along the line-of-sight. Because of this, the amplitude of the
galaxy-galaxy lensing spectra decreases, leading to worse constraints on
$\hat{J}$. Second, at high redshift, the lenses are necessarily closer to the
sources, which decreases the lensing kernel, and in addition, the lenses are
correlated with a smaller number of bins, which decreases the number of
independent measurements of $\hat{J}$ at that redshift. Let us for example
consider the tomographic bin number two for the lenses in Fig. 1, which is
centered at redshift $z\sim 0.4$. The lensing efficiency will peak at roughly
the double of this redshift, implying that the sources in the third
tomographic bin, which is centered at redshift $z\sim 0.9$, will provide a
high signal-to-noise measurement. In addition, this second bin is also
correlated (albeit less strongly) with the bins number two, four, and five of
the sources, providing four independent measurements of $\hat{J}$ at $z\sim
0.4$. On the other hand, if we consider the high-redshift lenses, like
tomographic bins number seven or higher (effective redshifts higher than
$z\sim 1.2$), there are no sources at the double of these effective redshifts,
where the signal would be the strongest. This means that only the very low-
redshift tail of the lenses distributions will be close to half the effective
redshift of the sources. Therefore, this will lead to a decrease of the
lensing kernel and thus a decrease of the galaxy-galaxy lensing signal.
Furthermore, since the lenses are not correlated with the sources at lower
redshift, the lenses bin number seven is only correlated with the sources bin
number four and five, which provides only two independent measurements of
$\hat{J}$. All together, these effects lead to worse constraints on $\hat{J}$
at high redshift. We note that all tomographic bins have the same number of
galaxies with the same ellipticity dispersion. Therefore, the shape and shot
noise are the same for all redshifts. Cosmic variance is instead a bit smaller
at high redshift, given the larger volume, but not enough to compensate the
effects mentioned above. However, even accounting for the degradation as a
function of redshift, LSST will be able to constrain $\hat{J}$ as a function
of redshift at less than percent level, and therefore probe the Weyl potential
with very high precision.

Let us mention that we have validated our method to measure $\hat{J}$, by
comparing it with the standard method, assuming $\Lambda$CDM. For this we have
proceeded in the following way. First, we have verified that replacing the
standard $\Lambda$CDM harmonic power spectra by the ones provided in Eqs.
(14,II.1) with the $\hat{J}$ and $\hat{b}$ values in Eqs. (15-16), we recover
the same constraints on the cosmological parameters at the level of $0.1\,\%$.
We note that in this test we have fixed the values of $\hat{J}$ and $\hat{b}$
and constrained the cosmological parameters, with the goal of testing the
implementation of the new harmonic spectra against the standard method. In a
second step, we have assumed $\hat{J}$ and $\hat{b}$ to be constant within
each redshift bin and moved them out of the integral. Their values have then
been fixed according to Eqs. (15-16) at the effective redshift of the
corresponding bin. Under this assumption, the recovered constraints on the
cosmological parameters degrade by a factor between 1 and 3.5 compared to the
constraints obtained with the standard approach. Such a discrepancy is
expected, since by fixing the value of $\hat{J}$ and $\hat{b}$ in each
redshift bin, we remove the information coming from the evolution of these
quantities inside the bins. Performing the same test with twice the number of
redshift bins, we have found that the degradation reduces to a factor between
1 and 1.8, showing that as we increase the number of bins and the
approximation of constant $\hat{J}$ and $\hat{b}$ is more valid, we recover
the standard constraints with our new implementation.

Besides validating our methodology, this test also shows the level of
degradation on the cosmological constraints due to our requirement of model-
independence. Since we want to measure the evolution of the Weyl potential
without assuming a specific theory of gravity, we are limited by the size of
the tomographic bins. Contrary to standard methods, where the evolution within
a redshift bin is given by the model, here we can only measure the value of
the Weyl potential at the effective redshift of the bins, therefore loosing
part of the information.

###  IV.2 Gravitational redshift

Let us now focus on the spectroscopic side of the analysis with the dipole
measurements. We present the main results in the bottom panel of Fig. 2. As
for $\hat{J}$ in the top panel, we present the $1\,\sigma$ forecast
uncertainty on $\hat{I}$ with respect to their fiducial value as a function of
redshift. We present both the optimistic (red) and pessimistic (blue)
scenarios, which correspond to $d_{\rm min}=20\,\text{Mpc}/h$ and $d_{\rm
min}=32\,\text{Mpc}/h$, respectively. As for the gravitational lensing
observable, the uncertainties on $\hat{I}$ increase as a function of redshift.
This is due to the fact that the signal decreases with increasing redshift,
given the decreasing growth at higher redshifts. Moreover, shot-noise
increases quickly as a function of redshift, as can be seen in the drop of the
number of galaxies in Fig. 1.

Overall, SKA2 will be able to constrain $\hat{I}$ as a function of redshift at
the level of $\sim 20\%$. These results are consistent with the forecasts
presented in Sobral-Blanco and Bonvin (2022), although different redshift bins
have been considered in this work. It is important to note that the
constraining power on $\hat{J}$ is much stronger than on $\hat{I}$, but this
is an expected result. The signal-to-noise ratio (SNR) of the galaxy-galaxy
lensing observable is already of 148 with current observations Prat  _et al._
(2022), while the expected SNR for gravitational redshift is much lower. Using
the bins defined in this analysis, the SNR can be computed as

$\displaystyle{\rm SNR}^{\Psi}=\sum_{z}\sum_{ij}\xi^{\Psi}_{1}(d_{i},z){\rm
cov}^{-1}(d_{i},d_{j},z)\xi^{\Psi}_{1}(d_{j},z)\,,$ (37)

where $\xi^{\Psi}_{1}(d_{i},z)$ is the gravitational redshift contribution to
the dipole, i.e. the contribution from the first term in Eq. (27), evaluated
at separation $d_{i}$ and in the redshift bin $z$, and cov$(d_{i},d_{j},z)$
stands for its covariance at separations $d_{i}$ and $d_{j}$ and redshift bin
$z$. Summing from 20 to 120 Mpc$/h$, we find a total SNR of 8. Note that the
SNR of the dipole is significantly larger, due to the Doppler effects, but
these terms do not contribute to the constraints on $\hat{I}$. Therefore, it
is not surprising that the constraints obtained with gravitational lensing
will be much more stringent. However, the particularity of the dipole is that
it allows us to directly constrain $\hat{I}$ and therefore directly probe
$\Psi$ in a model-independent way, something that cannot be done with any
other observable at cosmological scales.

Finally, let us stress that the constraining power on $\hat{I}$ directly
depends on the magnitude of the bias difference between the bright and the
faint population. In this work, we assumed a bias difference of 1, consistent
with what has been measured for BOSS in Gaztanaga  _et al._ (2017). If the
bias difference turns out to be smaller (see, e.g., the analysis of Alam  _et
al._ (2017), which finds a bias difference of order $\mathcal{O}(0.5)$), the
constraints presented in this work worsen by a factor of $\sim 2$. However,
the final results still provide stringent constraints.

###  IV.3 Anisotropic stress with the combination of probes

Figure 3: Reconstruction of $2/(1+\eta)$ as a function of redshift. The dashed
horizontal line represents the fiducial, while the vertical error bars show
the optimistic (red, light gray) and pessimistic (blue, dark gray) forecast
uncertainties.

After computing the forecast constraints for both $\hat{J}$ and $\hat{I}$, we
can combine them and place constraints on the anisotropic stress $\eta$.
Following the methodology described in Sec. III and after marginalizing over
the last three tomographic bins for $\hat{J}$, we generate a mock chain of 50
000 points for $\hat{J}$ and another one for $\hat{I}$. In each point we
compute the ratio presented in Eq. (30) and compute the posterior at each
redshift bin. We present in Fig. 3 the $1\,\sigma$ uncertainties on
$2/(1+\eta)$ as a function of redshift for both the optimistic (red) and
pessimistic (blue) scenarios. Again, an offset in the $x$ axis has been
included for illustrative purposes.

As it was the case for $\hat{J}$ and $\hat{I}$, there is a degradation of the
constraining power as a function of redshift, which is given by the same
physical effects since we are computing the ratio of the two quantities. We
can appreciate a more significant degradation compared to the $\hat{J}$ case,
but this is due to the fact that the uncertainties are driven by the
uncertainties on $\hat{I}$. Overall, the combination of LSST and SKA2 will
allow us to constrain the anisotropic stress, through the relation
$2/(1+\eta)$, in a model-independent way as a function of redshift and at the
level of $\sim 20\,\%$ in the optimistic scenario ($\sim 30\,\%$ in the
pessimistic case). Note again that if the bias difference is of order
$\mathcal{O}(0.5)$, the final constraints on $\eta$ worsen by a factor of
$\sim 2$. Finally, let us mention that since the uncertainties are essentially
dominated by the dipole, we do not expect a significant gain on the
constraints if we add cosmic shear into the analysis.

##  V Conclusions

In this work we have addressed the question of how well will future surveys
constrain the laws of gravity. Given how open this question is, one way to
answer it is by comparing the two gravitational potentials encoding the time
distortion and the spatial distortion, that is $\Psi$ and $\Phi$,
respectively. This comparison is a key test of the validity of GR, since it
compares directly the two independent degrees of freedom of the Universe’s
geometry. In contrast, the growth of structure is less direct: a deviation
from the $\Lambda$CDM prediction can indeed be generated not only by modified
gravity, but also by a fifth force acting on dark matter, or by a dark energy
component that is clustering. Because of this, measuring $\eta\neq 1$ is often
regarded as the smoking gun of modified gravity.

Many analyses exist in the literature constraining the anisotropic stress,
$\eta$, with current observations Abbott  _et al._ (2019); Alam  _et al._
(2021); Abbott  _et al._ (2022); Noller (2020), or forecasting the expected
constraints with future surveys Amendola  _et al._ (2013b); Motta  _et al._
(2013); Pinho  _et al._ (2018). However, these analyses all rely on the
validity of the weak equivalence principle. Hence, they are simply not valid
if dark matter obeys a fifth force. In this work we have built a new
estimator, which does not rely on the behavior (or even on the existence) of
dark matter, and we have forecast the expected constraints from future galaxy
surveys.

Our estimator improves the way gravity is tested in two ways. First, it relies
on a novel parametrization that we have developed for the galaxy-galaxy
lensing observable, which allows us to directly measure the evolution of the
Weyl potential, encoded in a free function $\hat{J}$, in each tomographic bin.
We have shown that combining the galaxy-galaxy lensing with the galaxy
clustering of the photometrically-detected galaxies allows us to break the
degeneracy between this new function and the galaxy bias. We have then
forecast the constraining capability of LSST with a Fisher matrix formalism
and determined the uncertainties on $\hat{J}$ as a function of redshift (see
Fig. 2). We have observed an expected degradation of the constraining power
toward higher redshifts, but the overall constraints are at less than percent
level, showing that LSST will be able to provide very precise model-
independent constraints of the late-time growth of the Weyl potential. This
novel parametrization has the strong advantage of being fully model-
independent, contrary to the standard $\mu-\Sigma$ parametrization, which
relies on the validity of the weak equivalence principle to measure $\mu$ and
$\Sigma$ from weak lensing observables Castello  _et al._ (2022); Bonvin and
Pogosian (2022).

The second improvement with respect to standard analyses is that our estimator
is built from a direct measurement of the evolution of $\Psi$ (encoded in a
new function $\hat{I}$) from the dipole of galaxy clustering. Using the Fisher
matrix formalism, we have forecast the constraining power of a future survey
like SKA2 on $\hat{I}$ as a function of redshift (see Fig. 2). We have
observed an expected degradation as a function of redshift and, overall,
significantly weaker constraints compared to $\hat{J}$. This is also expected
given the difficulty to measure the relativistic dipole of the correlation
function. Nevertheless, SKA2 will be able to put constraints of the order of
$20\,\%$ on $\hat{I}$, therefore directly constraining the gravitational
potential $\Psi$. These results are consistent with those presented in Sobral-
Blanco and Bonvin (2022), although a different redshift binning has been
considered in this work. This methodology has, again, the strong advantage of
being model-independent. Current measurements of $\Psi$ are instead performed
through RSD, assuming that the equivalence principle is valid. More precisely,
$\Psi$ is inferred from the velocity, using Euler’s equation.

Finally, we have combined the constraints on $\hat{J}$ from LSST and the
constraints on $\hat{I}$ from SKA2 to constrain the anisotropic stress. More
precisely, we have considered the combination $2/(1+\eta)$ and shown that
these two future galaxy surveys will be able to constrain it at the $\sim
20\,\%$ level in a direct and model-independent way. This level of precision
is similar to current constraints on $\eta$ obtained through the $\mu-\Sigma$
parametrization. The latest DES analysis does indeed constrain $\mu_{0}$ with
a precision of 20% and $\Sigma_{0}$ with a precision of 5% Abbott  _et al._
(2022), leading to a derived constraint on $\eta_{0}$ of the order of 25%.
However, as explained before, these constraints are not valid if dark matter
obeys a fifth force. Moreover, they assume a given time evolution for $\mu$
and $\Sigma$, meaning that only $\eta_{0}$, i.e. the value of $\eta$ today can
be constrained with this method. If this time evolution is not correct, then
the constraints that are obtained are invalid. In contrast, in our case,
$\eta$ is measured independently in each of the redshift bins of the surveys.

To conclude, let us mention that our analysis has considered a simple Fisher
matrix formalism, accounting for some astrophysical systematic uncertainties,
like galaxy bias or intrinsic alignments. However, an analysis using real
observations should go beyond this proof-of-concept and account for
observational systematic uncertainties, that we have assumed here to be under
control.

#### Acknowledgements.

This project has received funding from the Swiss National Science Foundation
and from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (Grant agreement No. 863929;
project title “Testing the law of gravity with novel large-scale structure
observables”).

## References

  * Adamek  _et al._ (2017) J. Adamek, R. Durrer, and M. Kunz, JCAP 11, 004 (2017), arXiv:1707.06938 [astro-ph.CO] .
  * Motta  _et al._ (2013) M. Motta, I. Sawicki, I. D. Saltas, L. Amendola, and M. Kunz, Phys. Rev. D 88, 124035 (2013), arXiv:1305.0008 [astro-ph.CO] .
  * Amendola  _et al._ (2018) L. Amendola _et al._ , Living Rev. Rel. 21, 2 (2018), arXiv:1606.00180 [astro-ph.CO] .
  * Song  _et al._ (2011) Y.-S. Song, G.-B. Zhao, D. Bacon, K. Koyama, R. C. Nichol, and L. Pogosian, Phys. Rev. D84, 083523 (2011), arXiv:1011.2106 [astro-ph.CO] .
  * Amendola  _et al._ (2013a) L. Amendola _et al._ (Euclid Theory Working Group), Living Rev. Rel. 16, 6 (2013a), arXiv:1206.1225 [astro-ph.CO] .
  * Abbott  _et al._ (2019) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 99, 123505 (2019), arXiv:1810.02499 [astro-ph.CO] .
  * Alam  _et al._ (2021) S. Alam _et al._ (eBOSS), Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Abbott  _et al._ (2022) T. M. C. Abbott _et al._ (DES), arXiv e-prints (2022), arXiv:2207.05766 [astro-ph.CO] .
  * Amendola (2000) L. Amendola, Phys. Rev. D 62, 043511 (2000), arXiv:astro-ph/9908023 .
  * Bonvin and Pogosian (2022) C. Bonvin and L. Pogosian, arXiv e-prints (2022), arXiv:2209.03614 [astro-ph.CO] .
  * Horndeski (1974) G. W. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Creminelli  _et al._ (2020) P. Creminelli, G. Tambalo, F. Vernizzi, and V. Yingcharoenrat, JCAP 05, 002 (2020), arXiv:1910.14035 [gr-qc] .
  * Noller (2020) J. Noller, Phys. Rev. D 101, 063524 (2020), arXiv:2001.05469 [astro-ph.CO] .
  * Sobral-Blanco and Bonvin (2021) D. Sobral-Blanco and C. Bonvin, Phys. Rev. D 104, 063516 (2021), arXiv:2102.05086 [astro-ph.CO] .
  * Bonvin  _et al._ (2014) C. Bonvin, L. Hui, and E. Gaztanaga, Phys. Rev. D 89, 083535 (2014), arXiv:1309.1321 [astro-ph.CO] .
  * Tansella  _et al._ (2018) V. Tansella, C. Bonvin, R. Durrer, B. Ghosh, and E. Sellentin, JCAP 03, 019 (2018), arXiv:1708.00492 [astro-ph.CO] .
  * Blanchard  _et al._ (2020) A. Blanchard _et al._ (Euclid), Astron. Astrophys. 642, A191 (2020), arXiv:1910.09273 [astro-ph.CO] .
  * Gleyzes  _et al._ (2015) J. Gleyzes, D. Langlois, M. Mancarella, and F. Vernizzi, JCAP 1508, 054 (2015), arXiv:1504.05481 [astro-ph.CO] .
  * Gleyzes  _et al._ (2016) J. Gleyzes, D. Langlois, M. Mancarella, and F. Vernizzi, JCAP 02, 056 (2016), arXiv:1509.02191 [astro-ph.CO] .
  * Raveri  _et al._ (2023) M. Raveri, L. Pogosian, M. Martinelli, K. Koyama, A. Silvestri, G.-B. Zhao, J. Li, S. Peirone, and A. Zucca, JCAP 02, 061 (2023), arXiv:2107.12990 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Takahashi  _et al._ (2012) R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Astrophys. J. 761, 152 (2012), arXiv:1208.2701 [astro-ph.CO] .
  * Castello  _et al._ (2022) S. Castello, N. Grimm, and C. Bonvin, Phys. Rev. D 106, 083511 (2022), arXiv:2204.11507 [astro-ph.CO] .
  * Camera (2022) S. Camera, arXiv e-prints (2022), arXiv:2206.03499 [astro-ph.CO] .
  * Bernardeau  _et al._ (2014) F. Bernardeau, T. Nishimichi, and A. Taruya, MNRAS 445, 1526 (2014), arXiv:1312.0430 [astro-ph.CO] .
  * Muir  _et al._ (2021) J. Muir _et al._ (DES), Phys. Rev. D 103, 023528 (2021), arXiv:2010.05924 [astro-ph.CO] .
  * Sobral-Blanco and Bonvin (2022) D. Sobral-Blanco and C. Bonvin, MNRAS 519, L39 (2022), arXiv:2205.02567 [astro-ph.CO] .
  * Bonvin (2014) C. Bonvin, Class. Quant. Grav. 31, 234002 (2014), arXiv:1409.2224 [astro-ph.CO] .
  * McDonald (2009) P. McDonald, JCAP 11, 026 (2009), arXiv:0907.5220 [astro-ph.CO] .
  * Croft (2013) R. A. C. Croft, MNRAS 434, 3008 (2013), arXiv:1304.4124 [astro-ph.CO] .
  * Bonvin and Durrer (2011) C. Bonvin and R. Durrer, Phys. Rev. D 84, 063505 (2011), arXiv:1105.5280 .
  * Yoo  _et al._ (2009) J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 083514 (2009), arXiv:0907.0707 [astro-ph.CO] .
  * Challinor and Lewis (2011) A. Challinor and A. Lewis, Phys. Rev. D 84, 043516 (2011), arXiv:1105.5292 [astro-ph.CO] .
  * Jelic-Cizmek  _et al._ (2021) G. Jelic-Cizmek, F. Lepori, C. Bonvin, and R. Durrer, JCAP 04, 055 (2021), arXiv:2004.12981 [astro-ph.CO] .
  * Lepori  _et al._ (2022) F. Lepori _et al._ (Euclid), Astron. Astrophys. 662, A93 (2022), arXiv:2110.05435 [astro-ph.CO] .
  * Gaztanaga  _et al._ (2017) E. Gaztanaga, C. Bonvin, and L. Hui, JCAP 01, 032 (2017), arXiv:1512.03918 [astro-ph.CO] .
  * Bonvin and Fleury (2018) C. Bonvin and P. Fleury, JCAP 05, 061 (2018), arXiv:1803.02771 [astro-ph.CO] .
  * Beutler and Di Dio (2020) F. Beutler and E. Di Dio, JCAP 07, 048 (2020), arXiv:2004.08014 [astro-ph.CO] .
  * Saga  _et al._ (2022) S. Saga, A. Taruya, M.-A. Breton, and Y. Rasera, MNRAS 511, 2732 (2022), arXiv:2109.06012 [astro-ph.CO] .
  * Ivezić  _et al._ (2019) Ž. Ivezić _et al._ (LSST), Astrophys. J. 873, 111 (2019), arXiv:0805.2366 [astro-ph] .
  * Laureijs  _et al._ (2011) R. Laureijs _et al._ , arXiv e-prints (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Spergel  _et al._ (2015) D. Spergel _et al._ , arXiv e-prints (2015), arXiv:1503.03757 [astro-ph.IM] .
  * Zuntz  _et al._ (2015) J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, and J. Kowalkowski, Astron. Comput. 12, 45 (2015), arXiv:1409.3409 [astro-ph.CO] .
  * Bridle and King (2007) S. Bridle and L. King, New J. Phys. 9, 444 (2007), arXiv:0705.0166 [astro-ph] .
  * Kirk  _et al._ (2012) D. Kirk, A. Rassat, O. Host, and S. Bridle, MNRAS 424, 1647 (2012), arXiv:1112.4752 [astro-ph.CO] .
  * Bull (2016) P. Bull, Astrophys. J. 817, 26 (2016), arXiv:1509.07562 [astro-ph.CO] .
  * Bonvin  _et al._ (2020) C. Bonvin, F. Oliveira Franco, and P. Fleury, JCAP 2020, 004 (2020), arXiv:2004.06457 [astro-ph.CO] .
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 [astro-ph] .
  * Howlett  _et al._ (2012) C. Howlett, A. Lewis, A. Hall, and A. Challinor, JCAP 2012, 027 (2012), arXiv:1201.3654 [astro-ph.CO] .
  * Prat  _et al._ (2022) J. Prat _et al._ , Phys. Rev. D 105, 083528 (2022), arXiv:2105.13541 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam, H. Zhu, R. A. C. Croft, S. Ho, E. Giusarma, and D. P. Schneider, MNRAS 470, 2822 (2017), https://academic.oup.com/mnras/article-pdf/470/3/2822/18314915/stx1421.pdf .
  * Amendola  _et al._ (2013b) L. Amendola, M. Kunz, M. Motta, I. D. Saltas, and I. Sawicki, Phys. Rev. D 87, 023501 (2013b), arXiv:1210.0439 [astro-ph.CO] .
  * Pinho  _et al._ (2018) A. M. Pinho, S. Casas, and L. Amendola, JCAP 11, 027 (2018), arXiv:1805.00027 [astro-ph.CO] .

## Appendix A

In this appendix we provide the fiducial values for all the parameters
considered in our analysis, including those that have been fixed. The latter
are presented in Table 2, while the former are shown in Table 3.

Table 2: List of fixed parameters considered in the Fisher analysis, together with their fiducial values in the $\Lambda$CDM model. Parameter |  Fiducial value
---|---
$\Omega_{\rm m,0}$ |  $0.3111$
$\Omega_{\rm b,0}$ |  $0.0490$
$h$ |  $0.677$
$n_{\rm s}$ |  $0.9665$
$A_{\rm s}$ |  $2.105\times 10^{-9}$
$s_{\text{B},1}$ |  $s_{\text{B}}(z=0.25)=0.3706$
$s_{\text{B},2}$ |  $s_{\text{B}}(z=0.38)=0.4665$
$s_{\text{B},3}$ |  $s_{\text{B}}(z=0.51)=0.5757$
$s_{\text{B},4}$ |  $s_{\text{B}}(z=0.65)=0.6817$
$s_{\text{B},5}$ |  $s_{\text{B}}(z=0.79)=0.7839$
$s_{\text{B},6}$ |  $s_{\text{B}}(z=0.95)=0.8974$
$s_{\text{B},7}$ |  $s_{\text{B}}(z=1.13)=1.0224$
$s_{\text{F},1}$ |  $s_{\text{F}}(z=0.25)=-0.1618$
$s_{\text{F},2}$ |  $s_{\text{F}}(z=0.38)=-0.1279$
$s_{\text{F},3}$ |  $s_{\text{F}}(z=0.51)=-0.1269$
$s_{\text{F},4}$ |  $s_{\text{F}}(z=0.65)=-0.1209$
$s_{\text{F},5}$ |  $s_{\text{F}}(z=0.79)=-0.1164$
$s_{\text{F},6}$ |  $s_{\text{F}}(z=0.95)=-0.1120$
$s_{\text{F},7}$ |  $s_{\text{F}}(z=1.13)=-0.1080$
Table 3: List of free parameters considered with their fiducial values in the $\Lambda$CDM model. All the parameters were varied in the Fisher analysis. Parameter |  Fiducial value
---|---
$\hat{J}_{1}$ |  $\hat{J}(z=0.25)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.25)=0.3388$
$\hat{J}_{2}$ |  $\hat{J}(z=0.38)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.38)=0.3666$
$\hat{J}_{3}$ |  $\hat{J}(z=0.51)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.51)=0.3846$
$\hat{J}_{4}$ |  $\hat{J}(z=0.65)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.65)=0.3949$
$\hat{J}_{5}$ |  $\hat{J}(z=0.79)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.79)=0.3979$
$\hat{J}_{6}$ |  $\hat{J}(z=0.95)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.95)=0.3948$
$\hat{J}_{7}$ |  $\hat{J}(z=1.13)=\Omega_{\rm m}\cdot\sigma_{8}(z=1.13)=0.3859$
$\hat{J}_{8}$ |  $\hat{J}(z=1.35)=\Omega_{\rm m}\cdot\sigma_{8}(z=1.35)=0.3707$
$\hat{J}_{9}$ |  $\hat{J}(z=1.7)=\Omega_{\rm m}\cdot\sigma_{8}(z=1.7)=0.3427$
$\hat{J}_{10}$ |  $\hat{J}(z=2.1)=\Omega_{\rm m}\cdot\sigma_{8}(z=2.1)=0.3111$
$\hat{b}_{1}$ |  $\hat{b}(z=0.25)=b\cdot\sigma_{8}(z=0.25)=1.4462$
$\hat{b}_{2}$ |  $\hat{b}(z=0.38)=b\cdot\sigma_{8}(z=0.38)=1.3512$
$\hat{b}_{3}$ |  $\hat{b}(z=0.51)=b\cdot\sigma_{8}(z=0.51)=1.2644$
$\hat{b}_{4}$ |  $\hat{b}(z=0.65)=b\cdot\sigma_{8}(z=0.65)=1.1796$
$\hat{b}_{5}$ |  $\hat{b}(z=0.79)=b\cdot\sigma_{8}(z=0.79)=1.1035$
$\hat{b}_{6}$ |  $\hat{b}(z=0.95)=b\cdot\sigma_{8}(z=0.95)=1.0259$
$\hat{b}_{7}$ |  $\hat{b}(z=1.13)=b\cdot\sigma_{8}(z=1.13)=0.9492$
$\hat{b}_{8}$ |  $\hat{b}(z=1.35)=b\cdot\sigma_{8}(z=1.35)=0.8683$
$\hat{b}_{9}$ |  $\hat{b}(z=1.7)=b\cdot\sigma_{8}(z=1.7)=0.7631$
$\hat{b}_{10}$ |  $\hat{b}(z=2.1)=b\cdot\sigma_{8}(z=2.1)=0.6690$
$A_{\rm IA}$ |  1.0
$\hat{I}_{1}$ |  $\hat{I}(z=0.25)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.25)=0.3388$
$\hat{I}_{2}$ |  $\hat{I}(z=0.38)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.38)=0.3666$
$\hat{I}_{3}$ |  $\hat{I}(z=0.51)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.51)=0.3846$
$\hat{I}_{4}$ |  $\hat{I}(z=0.65)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.65)=0.3949$
$\hat{I}_{5}$ |  $\hat{I}(z=0.79)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.79)=0.3979$
$\hat{I}_{6}$ |  $\hat{I}(z=0.95)=\Omega_{\rm m}\cdot\sigma_{8}(z=0.95)=0.3948$
$\hat{I}_{7}$ |  $\hat{I}(z=1.13)=\Omega_{\rm m}\cdot\sigma_{8}(z=1.13)=0.3859$
$\hat{f}_{1}$ |  $\hat{f}(z=0.25)=f\cdot\sigma_{8}(z=0.25)=0.4761$
$\hat{f}_{2}$ |  $\hat{f}(z=0.38)=f\cdot\sigma_{8}(z=0.38)=0.4826$
$\hat{f}_{3}$ |  $\hat{f}(z=0.51)=f\cdot\sigma_{8}(z=0.51)=0.4811$
$\hat{f}_{4}$ |  $\hat{f}(z=0.65)=f\cdot\sigma_{8}(z=0.65)=0.4733$
$\hat{f}_{5}$ |  $\hat{f}(z=0.79)=f\cdot\sigma_{8}(z=0.79)=0.4612$
$\hat{f}_{6}$ |  $\hat{f}(z=0.95)=f\cdot\sigma_{8}(z=0.95)=0.4444$
$\hat{f}_{7}$ |  $\hat{f}(z=1.13)=f\cdot\sigma_{8}(z=1.13)=0.4238$
$\hat{b}_{\text{B},1}$ |  $\hat{b}_{\text{B}}(z=0.25)=b_{\text{B}}\cdot\sigma_{8}(z=0.25)=0.8486$
$\hat{b}_{\text{B},2}$ |  $\hat{b}_{\text{B}}(z=0.38)=b_{\text{B}}\cdot\sigma_{8}(z=0.38)=0.8416$
$\hat{b}_{\text{B},3}$ |  $\hat{b}_{\text{B}}(z=0.51)=b_{\text{B}}\cdot\sigma_{8}(z=0.51)=0.8380$
$\hat{b}_{\text{B},4}$ |  $\hat{b}_{\text{B}}(z=0.65)=b_{\text{B}}\cdot\sigma_{8}(z=0.65)=0.8382$
$\hat{b}_{\text{B},5}$ |  $\hat{b}_{\text{B}}(z=0.79)=b_{\text{B}}\cdot\sigma_{8}(z=0.79)=0.8430$
$\hat{b}_{\text{B},6}$ |  $\hat{b}_{\text{B}}(z=0.95)=b_{\text{B}}\cdot\sigma_{8}(z=0.95)=0.8541$
$\hat{b}_{\text{B},7}$ |  $\hat{b}_{\text{B}}(z=1.13)=b_{\text{B}}\cdot\sigma_{8}(z=1.13)=0.8739$
$\hat{b}_{\text{F},1}$ |  $\hat{b}_{\text{F}}(z=0.25)=b_{\text{F}}\cdot\sigma_{8}(z=0.25)=0.1256$
$\hat{b}_{\text{F},2}$ |  $\hat{b}_{\text{F}}(z=0.38)=b_{\text{F}}\cdot\sigma_{8}(z=0.38)=0.1661$
$\hat{b}_{\text{F},3}$ |  $\hat{b}_{\text{F}}(z=0.51)=b_{\text{F}}\cdot\sigma_{8}(z=0.51)=0.2059$
$\hat{b}_{\text{F},4}$ |  $\hat{b}_{\text{F}}(z=0.65)=b_{\text{F}}\cdot\sigma_{8}(z=0.65)=0.2486$
$\hat{b}_{\text{F},5}$ |  $\hat{b}_{\text{F}}(z=0.79)=b_{\text{F}}\cdot\sigma_{8}(z=0.79)=0.2914$
$\hat{b}_{\text{F},6}$ |  $\hat{b}_{\text{F}}(z=0.95)=b_{\text{F}}\cdot\sigma_{8}(z=0.95)=0.3413$
$\hat{b}_{\text{F},7}$ |  $\hat{b}_{\text{F}}(z=1.13)=b_{\text{F}}\cdot\sigma_{8}(z=1.13)=0.3995$