Constraining f(R) gravity with cross-correlation of galaxies and cosmic microwave background lensing

We look for signatures of the Hu-Sawicki f(R) modified gravity theory, proposed to explain the observed accelerated expansion of the universe; in observations of the galaxy distribution, the cosmic microwave background (CMB), and gravitational lensing of the CMB. We study constraints obtained by using observations of only the CMB primary anisotropies, before adding the galaxy power spectrum and its cross-correlation with CMB lensing. We show that cross-correlation of the galaxy distribution with lensing measurements is crucial to breaking parameter degeneracies, placing tighter constraints on the model. In particular, we set a strong upper limit on $\log{\lvert f_{R_0}\lvert }<-4.61$ at 95% confidence level. This means that while the model may explain the accelerated expansion, its impact on large-scale structure closely resembles General Relativity. Studies of this kind with future data sets will probe smaller potential deviations from General Relativity.


Introduction
The cause of the late-time accelerated expansion of the Universe is one of the most profound problems facing modern cosmology (Riess et al. 1998;Perlmutter et al. 1999).Many theories have been proposed to explain this phenomena, the most popular due to its apparent simplicity is Λ, the cosmological constant.However, the interpretation of the cosmological constant as the energy of the vacuum results in theoretical predictions that are at least 55 orders of magnitude too large (e.g., Carroll et al. 2004;Solà 2013), motivating the study of alternative explanations.One possibility is new gravitational physics (Carroll et al. 2004).In addition to causing the accelerated expansion, modifications to gravity may alter structure formation in the Universe.We constrain a particular model of f (R) gravity, which can explain the accelerated expansion of the Universe, using observations of galaxy clustering, weak gravitational lensing of the cosmic microwave background (CMB), and temperature and polarization information from the CMB.
Deviations from General Relativity (GR) are tightly constrained on the scales of our Solar System (Everitt et al. 2011;Will 2014).Therefore, models of modified gravity must at the same time satisfy these constraints on small-scales, whilst modifying gravity on large scales to explain the cosmic acceleration.Carroll et al. (2004) presented a general class of models that can drive cosmic acceleration, by replacing the linear dependence of the Einstein-Hilbert action on the Ricci scalar R with a nonlinear function of R (R → R + f (R)).Hu & Sawicki (2007) (in the following, HS) presented a class of f (R) models capable of explaining the cosmic acceleration and evading the strong Solar System constraints through a chameleon mechanism (Khoury & Weltman 2004;Navarro & Van Acoleyen 2007;Faulkner et al. 2007).
Constraints have been placed on HS f (R) gravity using many different complementary observations.Such observations con-strain f R 0 , the value of the cosmological field today, that we introduce in more detail in Sect. 2. In particular, on cosmological scales HS f (R) has been constrained by Cataneo et al. (2015), who obtained the constraint log| f R 0 |< −4.79 at the 95% confidence level, using cluster number counts in addition to CMB, supernovae and BAO data.Hu et al. (2016) also found log| f R 0 |< −4.5 using the CMB (temperature, polarisation and lensing), supernovae, BAO and galaxy weak lensing measurements.Hojjati et al. (2016) obtained the upper bound log| f R 0 |< −4.15 at the 95% confidence level using similar observations.
The strongest constraints come from galactic scales.Naik et al. (2019) were able to exclude log| f R 0 |> −6.1 using galaxy rotation curves, and Desmond & Ferreira (2020) excluded log| f R 0 |> −7.85 based on the analysis of galaxy morphology.Astrophysical and cosmological constraints on HS f (R) gravity can also be found in the review by Lombriser (2014).Finally, Casas et al. (2023) forecasts the constraints that will be achievable using observations from Euclid.Despite the strong constraints on HS f (R) gravity from galactic studies, it is still a useful model to explore deviations from GR on cosmological scales.
Many different tools have been developed to predict the matter power spectrum in HS f (R) gravity.Boltzmann codes that calculate the linear matter power spectrum are mgcamb (Zhao et al. 2009;Hojjati et al. 2011;Zucca et al. 2019;Wang et al. 2023) and MGCLASS (Sakr & Martinelli 2022).There are several simulation based emulators of the matter-power spectrum into the mildy non-linear regime (Winther et al. 2019;Ramachandra et al. 2021;Arnold et al. 2022;Sáez-Casares et al. 2023), and ReACT (Bose et al. 2020, 2023), which uses a halo model reaction framework validated on N-body simulations.
In the following section we review HS f (R) gravity.In Sect. 3 we introduce the observations used in our analysis, and then in Sect. 4 we overview our methodology; the estimation of the angular power spectrum, our covariance matrix estimation, and our likelihood.Results are presented in Sect.5, before the conclusions in Sect.6.

Hu-Sawicki f(R) gravity
In f (R) theories of gravity, the Einstein-Hilbert action is modified such that R → R + f (R); therefore, the action becomes, where R is the Ricci scalar, κ = 8πG with G the gravitational constant (and the speed of light set to 1), g is the determinant of the spacetime metric, L m is the matter Lagrangian, and f is a function of the Ricci scalar.In HS (Hu & Sawicki 2007), f (R) follows a broken power law, where m is a mass scale given by m2 = κ 2 ρm /3 with ρm the mean matter density of the Universe, and c 1 , c 2 and n are three dimensionless constants.The derivative of f with respect to the Ricci scalar R is denoted, and can be interpreted as a new scalar field.Hu & Sawicki (2007) showed that a background close to ΛCDM can be recovered by imposing, where Ω Λ and Ω m are the present-day dark energy and matter densities (divided by the critical density) in the ΛCDM cosmology.Imposing this relation, there remain only two free parameters in Eq. ( 2): n and either c 1 or c 2 .In the high curvature regime (R ≫ m 2 ), which Oyaizu (2008) showed to be the appropriate regime, Eq. ( 3) can be written as, which, evaluated at the present-day background, leads to, The f R 0 parameter denotes the background value of f R at the present time, which we choose as our free parameter to constrain the model of HS f (R).Additionally, we fix n = 1.
For the small values of f R 0 probed in this work, the background expansion of the Universe is indistinguishable between f (R) and that of a cosmological constant.Instead, we constrain f (R) through its impact on the growth of structure.This can be seen by looking at the modified Poisson equation in f (R), where a is the cosmological scale factor and δρ m ≡ ρ m − ρm .We see directly that f R /2 can be seen as the potential of the modified gravity force.As mentioned in the introduction, this modified Poisson equation approaches the GR expression within the Solar System through the chameleon mechanism (Khoury & Weltman 2004;Hu & Sawicki 2007).
It is worth noting that unlike other theories of modified gravity, HS f (R) has little effect on the propagation of light in the weak-field limit (for example Hojjati et al. 2016).

BOSS galaxies
We use the DR12 data release of the BOSS survey from the SDSS collaboration (Alam et al. 2015).This large-scale spectroscopic survey was divided into two subsamples, LOWZ and CMASS.LOWZ contains galaxies at low redshift, up to approximately z ≃ 0.45, while CMASS contains higher redshift galaxies, (roughly up to z ≃ 0.8) and was constructed to create a sample of galaxies with approximately constant stellar mass.As in Loureiro et al. (2019), we restrict these samples to 0.15 < z < 0.45 and 0.45 < z < 0.8 for LOWZ and CMASS, respectively, such that the two samples do not overlap in redshift.This choice allows us to neglect the covariance between galaxies of LOWZ and CMASS.Using these redshift ranges, and after masking regions of low completeness, the two samples contain 366, 576 and 751, 067 galaxies.The redshift distribution of the two samples is shown in Fig. 1.The mask and map making is identical to that within Kou & Bartlett (2023), which follows Reid et al. (2016) and Loureiro et al. (2019).We transform the MANGLE1 (Swanson et al. 2008) acceptance and veto masks, provided with the galaxy catalogs, into high resolution binary masks in HEALPix 2 (Górski et al. 2005;Zonca et al. 2019) format with N SIDE = 8192.The acceptance mask represents the completeness of the observations, while the veto mask excludes regions that could not be observed.A first cut is made to exclude regions with completeness below 0.7, before degrading the resolution of the mask to N SIDE = 4096.A second cut is then applied, such that pixels with completeness below 0.8 are rejected.Finally, the galaxy maps are computed by summing the weighted number of galaxies in each pixel, divided by the completeness of the pixel.The weight w tot that is applied to each galaxy takes into account a number of observational effects, including fibre collisions, redshift failures, stellar density and seeing conditions (more details can be found in Ross et al. 2012).The galaxy overdensity map is, where, where C p pix is the completeness in pixel p.

CMB temperature and polarization observations
We use observations of the CMB temperature and polarization anisotropies from the Planck satellite, which observed the CMB for about 29 months and covered the full sky.In this work, we directly make use of the likelihood code provided (Planck Collaboration et al. 2020b) and whose cosmological results were analyzed in Planck Collaboration et al. (2020a).

CMB lensing convergence map
The observed CMB fluctuations are distorted as the CMB photons traverse the Universe because of gravitational lensing.An observational consequence of this is the correlation between different multipoles in both the temperature and polarization anisotropies, which would not be present in the unlensed CMB.
The CMB lensing potential can therefore be reconstructed from such correlations (see Lewis & Challinor (2006) for a comprehensive review).
We use the CMB lensing convergence map released by Planck Collaboration et al. (2020c).This map was obtained using a minimum variance quadratic estimator based on temperature and polarization maps.This map covers about 67% of the sky and led to the detection of lensing at 40σ.The map is provided with resolution N SIDE = 4096, together with the associated mask with N SIDE = 2048.

Theoretical angular power spectra
We calculate the matter power spectrum in f (R) using two different codes, MGCLASS and ReACT.
MGCLASS (Sakr & Martinelli 2022) is a modified version of the Boltzmann code CLASS (Blas et al. 2011) in which the equations of the linear perturbation theory are changed to take into account modifications to gravity.It can therefore be used to predict the linear matter power spectrum.
ReACT (Bose et al. 2020(Bose et al. , 2023) ) gives predictions for the nonlinear matter power spectrum in beyond ΛCDM cosmologies, including wCDM, f (R) and DGP gravity.ReACT uses a halo model based approach described in Cataneo et al. (2019) such that, where P NL is the non-linear matter power spectrum in modified gravity, and is the so-called non-linear "pseudopower spectrum".This pseudo-power spectrum is defined as a ΛCDM power spectrum with initial conditions chosen such that the ΛCDM linear matter power spectrum matches the modified gravity linear matter power spectrum at a given redshift.This choice was made in order to ensure that the halo mass function in ΛCDM and in the modified gravity theory are similar (which is anticipated since they have been defined to have exactly the same linear matter power spectrum).
The remaining term in Eq. ( 10), R is called the reaction and describes how the ΛCDM matter power spectrum changes due to the modifications to gravity.The reaction is calculated using the halo model and 1-loop perturbation theory.More details can be found in Cataneo et al. (2019) and Bose et al. (2023).When using ReACT to predict the non-linear modified gravity matter Effect of changing the value of log| f R 0 | on the auto-power spectrum of CMASS, with all other parameters fixed.The bottom panel shows the difference relative to the non-linear power spectrum in ΛCDM, together with the 1σ uncertainties of CMASS.Plain lines are obtained using ReACT and the dashed lines with MGCLASS.The grey dotted line shows the limit between the linear and non-linear regimes.We limit our analysis to the linear regime, but show the theoretical predictions with ReACT.We do not show the predictions with MGCLASS, which fail in this regime.power spectrum, it is required to provide a reliable non-linear ΛCDM matter power spectrum, for which we use the halo model based HMCode (Mead et al. 2015).
For our observations we compute the galaxy auto power spectrum C gg ℓ , the CMB lensing convergence auto power spectrum C κκ ℓ and the cross-correlation between the two C κg ℓ .We also compute the CMB temperature and polarization power spectra, which are sensitive to f (R) gravity through the integrated Sachs-Wolfe effect (ISW) and gravitational lensing.The CMB temperature, polarization and convergence power spectra are predicted by MGCLASS.
For the galaxy auto and cross power spectra, we use the matter power spectrum prediction from either ReACT or MGCLASS.The angular power spectra are then modeled using the Limber approximation (Limber 1953), where H(z) is the Hubble parameter at redshift z, c is the speed of light, χ is the comoving distance, P m is the matter power spectrum, and k is the comoving wavenumber.Finally, W g and W κ are the galaxy and CMB lensing kernels, Here, b g is the galaxy bias, (1/n tot )(dn/dz) is the normalized galaxy redshift distribution, H 0 is the present value of H, Ω m is the matter density parameter, and z * is the redshift of the surface of last scattering.
In Fig. 2 we show the effect of changing log| f R 0 | on the galaxy angular power spectrum, using the galaxy redshift distribution of CMASS (see Sect. 3 for more details).It can be seen that HS f (R) gravity increases the formation of structure on small scales, leading to more power in the angular power spectrum at larger multipoles.For a given value of log| f R 0 |, MG-CLASS predicts slightly more power than ReACT, except at the highest multipoles, where power might be missing in the prediction of MGCLASS as MGCLASS only predicts the linear matter power spectrum.This is also the reason why we do not show the predictions using MGCLASS after the dotted grey line marking the transition into the non-linear regime.The non-linear regime is not used in our analysis as we do not use a theoretical model that can reliably model the galaxy bias in this regime.

Angular power spectra estimation
The angular cross-correlation power spectrum of two fields A and B is defined as, where a ℓm and b ℓ ′ m ′ are respectively the spherical harmonic coefficients of fields A and B, b * denotes the complex conjugate of b, δ is the Kronecker symbol, and ⟨ ⟩ is the ensemble average.Given full sky coverage, this power spectrum can be estimated using, In practice, however, fields A and B are not observed on the full sky, but on a limited sky fraction defined by their masks W A and W B , such that what we really observe is Ã It can then be shown (Hivon et al. 2002;Brown et al. 2005) that taking the ensemble average of Eq. ( 16) with the observed fields Ã and B gives, where M ℓℓ ′ is a coupling matrix that depends on the masks W A and W B .It is then possible to recover an unbiased estimate of C AB ℓ by inverting the coupling matrix.
In our analysis we use the public code NaMaster 3 (Alonso et al. 2019) to estimate the power spectra of CMASS and LOWZ, as well as their cross-correlation with the CMB lensing convergence map from Planck.For the CMB lensing auto-correlation, we take directly the Planck lensing likelihood.We also apodize the CMB lensing mask with a scale of 10 arcmin, and we have verified that the estimated power spectra do not depend on the apodization scale.Table 1.Scale cuts for each power spectrum such that the relative difference between the linear and non-linear power spectra is inferior to 5%.

Noise removal
The measured galaxy angular power spectrum is biased by the shot-noise contribution.Therefore, this contribution is subtracted from the estimated power spectrum.The galaxy shot noise is given by, where N is the weighted number of galaxies in each sample.

Scale cuts
As mentioned in Sect.4.1, for our cosmological constraints we only use power spectra in the linear regime.This limitation mainly comes from the fact that we use a linear galaxy bias and that this simple modeling is not reliable in the non-linear regime.
To determine the maximum multipole that can be used, we follow the approach of Loureiro et al. (2019).Namely, we use the fiducial cosmology of Planck Collaboration et al. ( 2020a) and predict the theoretical linear and non-linear power spectra in the ΛCDM model.The non-linear power spectra uses halofit (Smith et al. 2003;Takahashi et al. 2012) to model the non-linear matter power spectrum and we keep using a linear galaxy bias.We then determine the transition between the linear and non-linear regimes to correspond to the largest multipole ℓ max such that the relative difference between the linear and non-linear power spectra is inferior to 5%.The determined scales are presented in Tab. 1.In practice we use ℓ max = 200 for CMASS, ℓ max = 100 for LOWZ, and ℓ max = 400 for CMB lensing.
Finally, since we use the Limber approximation, which is not valid on large scales, we limit our study to multipoles above ℓ min = 20.

Covariance matrix
We use a Gaussian covariance matrix following Saraf et al. (2022).This allows us to incorporate non-overlapping regions of the sky within our cross-correlation analysis, where A,B,C,D label one of the two galaxy density fields g or the CMB lensing convergence field κ, and f AB sky is the sky fraction common to fields A and B.
We estimate an initial covariance matrix using Eq. ( 19) with the observed power spectra.This covariance matrix is used to fit our theoretical model to the observations, and we then determine a second covariance matrix using the best fit theoretical power spectra from the first analysis.This procedure reduces our sensitivity to the noise in the estimated power spectra.

Likelihood
Our log-likelihood is the sum of the log-likelihood of the galaxy power spectra and the galaxy -CMB lensing cross-correlations, which we denote here as ln L 2×2 pt , the Planck log-likelihood for temperature and polarization (Planck Collaboration et al. (20) The name 2 × 2 pt comes from the fact that it is the combination of two different types of two point statistics (C gg ℓ and C κg ℓ ).The third two point statistic is the Planck lensing power spectrum, C κκ ℓ , making our analysis a 3 × 2 pt analysis in combination with the Planck temperature and polarisation likelihood.
We adopt a Gaussian for the 2 × 2 pt likelihood, ln where X T denotes the transpose of vector X, θ is the parameter vector, C is the covariance matrix, and X is a concatenation of power spectra such that, The likelihood L Planck contains the likelihood for the TT, TE and EE power spectra.We have neglected the covariance between C κκ ℓ and C κg ℓ (as is done, for instance, in Abbott et al. 2023), which is motivated by the fact that the C κκ ℓ is estimated on a much larger sky fraction than C κg ℓ , and that CMB lensing power comes mainly from redshifts greater than 0.8, where we do not have any galaxy in either of the two samples.

Priors
We use flat priors for all of the cosmological parameters, namely ω m , ω b , h, τ reio , log 10 10 A s , n s , as well as for the two galaxy bias parameters, for LOWZ and CMASS, and log| f R 0 |, for which we impose −7 < log| f R 0 |< 0. We use the recommended priors for the calibration and nuisance parameters required by Planck's likelihood.In total, we end up with 30 free parameters (the 6 parameters of ΛCDM, log| f R 0 |, 2 galaxy bias parameters and 21 calibration and nuisance parameters).We then use the MCMC (Monte Carlo Markov Chain) sampler emcee4 (Foreman-Mackey et al. 2013) to sample the resulting posterior distributions.

Computing prior-independent constraints
As will be seen in Sect.5, in many cases we are only able to put upper bounds on log| f R 0 |.It is then non-trivial to estimate the 95% constraint on log| f R 0 |, since the percentile depends strongly on the lower value of the prior.This arises because the posterior on log| f R 0 | at low values is non zero (and almost flat).Therefore, without a lower bound on the prior, very low values of log| f R 0 | would be explored, rather than more interesting regions of the  posterior, which would in turn be poorly sampled.This choice of lower limit for log| f R 0 | changes the estimated 95 th percentile.
In order to resolve this issue, we follow the approach of Piga et al. ( 2023), who rely on Gordon & Trotta (2007).We consider the ratio of the marginalized posterior and our prior, where x is the parameter we wish to constrain (in our case log| f R 0 |), d is the data, p the prior, and P the posterior.Then, for two different values of x, say x 1 and x 2 , we apply Bayes theorem to obtain, where B is the Bayes factor and L(d|x) is the marginalized likelihood of the data for x.The Bayes factor, which quantifies the support for the model with x = x 1 over the model with x = x 2 , is therefore prior independent.Gordon & Trotta (2007) showed that B(x 1 , x 2 ) = 2.5 means that the model with x = x 1 is favored compared to the model with x = x 2 at 95%.We then fix x 1 to its lower bound of log| f R 0 |= −7 and find the value of x 2 such that B(x 1 , x 2 ) = 2.5.In this way, we are able to compute 95% confidence intervals which do not depend on the prior.We check this is indeed the case in Sect.5.2.4,where we vary the prior used in the analysis.

Results
In this section we present our results for different combinations of observations, with different choices of binnings, theoretical power spectra and nuisance parameters.The upper limits we place on log| f R 0 | are shown in Tab. 2 and, the corner plot for the fiducial 3 × 2 pt + CMB analysis, obtained thanks to the use of GetDist5 (Lewis 2019), is shown in Fig. 6.

TTTEEE + A lens TTTEEE + lensing + A lens TTTEEE + lensing + A lens in TTTEEE only
Fig. 4. Degeneracy between A lens and log| f R 0 | for CMB temperature and polarization only (blue); CMB temperature, polarization and lensing (green); and when A lens is considered as a systematic effect and is applied only to the primary CMB anisotropies (magenta).The contours correspond to the 68 th and 95 th percentiles of the posterior samples.

Results from CMB only
We first look at the constraints obtained on log| f R 0 | when using only observations of the CMB.All these constraints are obtained using MGCLASS, as ReACT can only make predictions for the matter power spectrum at redshifts z < 2.5.
Firstly, we use the CMB temperature and polarization power spectra.They are sensitive to f (R) gravity through the imprint of the integrated Sachs-Wolfe effect and gravitational lensing.Lensing is the dominant effect and causes a smoothing of the acoustic peaks (for more detail, see Lewis & Challinor 2006).This is distinct from the CMB lensing convergence reconstruction described in Sect.3.3.Since f (R) increases the power in the lensing potential, f (R) models will have a larger smoothing effect on the CMB.
With just the CMB temperature and polarization observations, we find a strong preference for a non-zero value of f R 0 at more than 3σ: a best-fit value of log| f R 0 |= −2.35, with the bounds −3.09 < log| f R 0 |< −1.88 at 95% and −5.76 < log| f R 0 |< −1.64 at 99.7% confidence levels.The marginalized posterior for log| f R 0 | is shown in Fig. 3, where we see a clear peak at log| f R 0 |= −2.35.This is not a new result and has been shown by other authors, for example Dossett et al. (2014) and Hojjati et al. (2016).
However, when we add the CMB lensing power spectrum, which is estimated from the mode mixing in the primary CMB (as described in Sect.3.3), the large values of log| f R 0 | are excluded, and instead we set an upper limit on log| f R 0 |: log| f R 0 |< −2.31 at 95% confidence (computed following the approach described in Sect.4.8).
We clearly see that the CMB lensing convergence power spectrum is consistent with GR and low log| f R 0 | values, while the smoothing of the acoustic peaks in the CMB temperature and polarization anisotropies prefers a higher value of log| f R 0 |.This issue (or tension) is closely related to the Planck A lens tension (Planck Collaboration et al. 2020a).
The A lens parameter was introduced in Calabrese et al. ( 2008) as a phenomenological parameter scaling the CMB lensing potential amplitude as, It therefore changes the amplitude of the CMB lensing convergence power spectrum, and also the smoothing in the temperature and polarization power spectra.It was shown in Planck Collaboration et al. (2020a) that the CMB lensing convergence power spectrum is perfectly compatible with A lens = 1, which is not the case of the temperature and polarization power spectra.For instance, the 1σ constraint that Planck Collaboration et al. (2020a) find using the TT,TE,EE+low E likelihood is A lens = 1.180 ± 0.065, which is in tension with A lens = 1 at 2.8σ.Interestingly, the log| f R 0 | tension is larger than the A lens tension.This suggests that the modification of the lensing potential by HS f (R) better describes the smoothing of the CMB than the simple rescaling of the potential with A lens .
We further elucidate the relation between A lens and log| f R 0 | by running the Planck CMB TTTEEE likelihood with both parameters.Figure 3 shows the result as the blue curve.The constraints on log| f R 0 | are broadened, with low values of log| f R 0 | now acceptable, and compatible with GR.The posterior, however, remains consistent with higher values of log| f R 0 |. Figure 4 illustrates the degeneracy between A lens and log| f R 0 |.
The green and magenta curves in Fig. 3 show the posteriors obtained when adding the CMB lensing convergence power spectrum likelihood.The difference between the two curves is that for the magenta curve, the effect of A lens is only applied to the temperature and polarization power spectra; in this case, we are modeling A lens > 1 not as a physical effect, but rather as an unknown systematic in the observations.In both cases, the posteriors agree with GR, and we can put an upper limit on the value of log| f R 0 |.
We show the degeneracy between log| f R 0 | and A lens in these two cases in Fig. 4. When A lens is treated as a systematic (magenta contours), we recover high A lens values because the convergence power spectrum constrains log| f R 0 | to low values where it cannot reproduce the smoothing observed in the acoustic peaks.These results clearly show that HS f (R) gravity cannot explain the A lens tension, in agreement with the findings of Hojjati et al. (2016).

Results from combined CMB and 3 × 2 pt observations
We now present our fiducial analysis, the combination of the CMB (temperature and polarization power spectra) and the 3 × 2 pt analysis (CMB convergence power spectrum, galaxy power spectrum and the cross-correlation between the two).In this section, the value of A lens is fixed to 1 if not stated otherwise.
Figure 5 gives the marginalized posterior of log| f R 0 | obtained when using either MGCLASS (black curve) or ReACT (red curve) when modeling the 2 × 2 pt observables (C gg ℓ and C κg ℓ ) and we recall that the CMB lensing convergence auto power spectrum is always modeled using MGCLASS.The 95% confidence level constraints are log| f R 0 |< −4.12 and log| f R 0 |< −4.61 with MG-CLASS and ReACT, respectively.These constraints are consistent with GR and are much tighter than when using only the CMB observations.This also excludes HS f (R) as an explanation of the A lens tension.
The difference between MGCLASS and ReACT is also evident in Fig. 2. MGCLASS only predicts the linear matter power 6 4 2 0 5. Marginalized posterior of log| f R 0 | using observations of CMB temperature and polarization, and the 3×2 pt observables (CMB lensing, galaxy distribution and their cross-correlation).The CMB observations including the CMB lensing auto power spectrum are always modeled using MGCLASS, while the galaxy power spectra and the galaxy -CMB lensing cross-correlations are modeled with ReACT (red) or MGCLASS (black).spectrum, whereas ReACT predicts the non-linear power spectrum.Although our scale cuts were chosen to minimise the effects of non-linear structure formation, it seems likely that the origin of this difference here arises from the mildly non-linear regime where the MGCLASS predictions have less power than those of ReACT.Therefore, our MGCLASS constraint is likely too conservative and hence not as strong as it should be.
Figure 6 shows the marginalized two-dimensional posteriors on all parameters when using MGCLASS (blue) or ReACT (red).We see that the posteriors agree well for most parameters.We also see that log| f R 0 | exhibits a degeneracy with the galaxy bias parameters.This arises because both parameters change the amplitude of the galaxy power spectra.The parameters are not completely degenerate, however, as log| f R 0 | also changes the shape of the power spectra, and the cross-correlation separates the two effects to a certain extent, as explored in the following subsection.
We show in Fig. 7 the measured angular power spectra (orange) and the theoretical best fit (blue) using MGCLASS.The model fits well the data except at the largest scales for the crosscorrelation between the CMB lensing of Planck and the galaxies of CMASS, where the amplitude of the theoretical power spectrum is larger than the amplitude of the observation.This result has been observed in previous studies (Pullen et al. 2016;Singh et al. 2017;Kou & Bartlett 2023).We see that HS f (R) is unable to resolve this problem.
The constraints on log| f R 0 | are consistent and competitive with previous studies using galaxy clustering observations, such as Hu et al. (2016).Their tightest constraint is log| f R 0 |< −4.5 when combining observations of the CMB (temperature, polarization, lensing), supernovae, baryon acoustic oscillation (BAO) measurements (including, but not limited to, BAO measurements of LOWZ and CMASS) and galaxy weak lensing shear correlation functions estimated from the Canada-France-Hawaii Telescope Lensing Survey (Heymans et al. 2013).The crosscorrelation of galaxy -CMB lensing spectra enables us to obtain competitive constraints with a reduced data-set.

Benefit of the cross-correlation
In order to isolate the advantage of the cross-correlation, we performed the analysis without the cross-correlation power spectra.The cross-correlation is primarily useful in reducing the degeneracy between log| f R 0 | and the galaxy bias parameters.The degeneracy is reduced as C gg ℓ is proportional to b 2 g while C κg ℓ is proportional to b g (see Eq. 13 and Eq. 14).This effect is seen in Fig. 8 where the red contours do not contain the cross-correlation and the blue contours do.As a result, the 95% constraint obtained without the cross-correlation is log| f R 0 |< −2.95, compared to log| f R 0 |< −4.12 with the cross-correlation, i.e., more than an order of magnitude improvement.

Result including A lens as a systematic
As was shown in Sect.5.1, HS f (R) cannot explain the A lens tension.When the 2 × 2 pt observations are added, log| f R 0 | is constrained to even smaller values, ruling out further this kind of modified gravity resolution to the A lens tension.Since it has been suggested that A lens could be due to a systematic error (Planck Collaboration et al. 2020a), rather than a physical effect, we also perform our analysis with A lens on the CMB temperature and polarization power spectra only.We expect this to give tighter constraints, because A lens can explain the excess of smoothing in the CMB temperature and polarization power spectra without the need for a large value of log| f R 0 |.The resulting marginalized contours on A lens and log| f R 0 | are shown in Fig. 9.As expected, A lens > 1 is preferred.Unsurprisingly, the 95% constraint of log| f R 0 |< −4.24 is tighter than the fiducial constraint (log| f R 0 |< −4.12).

Effect of the binning scheme
We examine how changing the binning scheme impacts our constraints.In particular, we defined three binning schemes, namely ∆ℓ = 10, ∆ℓ = 20 (the fiducial case) and a unequal binning scheme with average ∆ℓ = 35.The marginalized posteriors on log| f R 0 | obtained with the three different binning schemes are presented in Fig. 10.We see that the constraints become tighter as ∆ℓ decreases.The use of narrow multipole bin-widths makes it possible to better use the shape of the power spectra to constrain f (R) gravity.The 95% confidence constraints are log| f R 0 |< −4.18, log| f R 0 |< −4.12, and log| f R 0 |< −3.84 for ∆ℓ = 10, ∆ℓ = 20, and ∆ℓ = 35, respectively.All of our constraints are summarized in Tab. 2.

Effect of the priors
Our 95% confidence intervals rely on the approach described in Sect.4.8, which aims at building prior-independent confidence intervals.In order to check that our constraints indeed are prior independent, we ran a second MCMC analysis in the fiducial setup, with a flat prior, −9 < log| f R 0 |< 0 (instead of the previous −7 < log| f R 0 |< 0).The posteriors (normalized by the average value of the plateau) are shown in Fig. 11.It can be seen that the 95% constraint is largely insensitive to the lower bound on the prior.More precisely, the constraint using the new prior is log| f R 0 |< −4.17, which is close to the limit log| f R 0 |< −4.12 that we found with the previous prior.
The approach described in Sect.4.8 indeed produces constraints independent of the adopted prior.In contrast, if we directly use the 95 th percentile of the samples, the values associ-

ReACT mgclass
Fig. 6.Constraints and degeneracies on all parameters when using observations of the CMB and 3 × 2 pt observables, modeled with MGCLASS (blue) or ReACT (red).The contours correspond to the 68 th and 95 th percentiles of the posterior samples.
ated to each prior (lower bound at −7 or −9) would be, respectively, −3.93 and −4.27.This dependence on the choice of prior demonstrates the danger in determining the constraints directly from the sample percentiles; for instance, we extract a tighter constraint when the lower bound of the prior is lower.This finding is intuitive, since by decreasing the lower bound of the prior, we allow the MCMC walkers to explore lower values in parameter space; the posterior distribution thus shifts towards lower parameter values and, consequently, the 95 th percentile becomes lower as well.

Conclusion
Modified gravity is a possible explanation for the observed accelerated expansion of the Universe (Carroll et al. 2004).Hu-Sawicki (HS) f (R) gravity (Hu & Sawicki 2007) is an attractive example, motivating the search for other possible observational signatures of the model.We searched for such signatures as deviations from the predictions of General Relativity for large-scale structure observations in a combined analysis of CMB, galaxy, and CMB lensing measurements.
If the HS f (R) model is to explain the accelerated expansion, the key parameter is log| f R 0 |.Primary CMB observations alone constrain this parameter through the ISW effect and smoothing of the temperature and polarization anisotropies by gravitational lensing.In agreement with previous analyses, we find that measurements by Planck prefer high values of log| f R 0 |, which would imply a remarkable deviation from General Relativity (see Sect. 5.1).However, this preference disappears when the CMB lensing convergence power spectrum is added, reflecting tension between the effects of lensing on the primary anisotropies and the reconstructed lensing power spectrum.This C g, LOWZ Fig. 7. Measured angular power spectra (orange) and best fit (blue) obtained using MGCLASS, as a function of multipole.The grey dashed line represents the limit between the linear and non-linear regimes.We only used multipoles above this limit in this analysis.Note the different multipole ranges for CMASS and LOWZ.

Case
Upper tension is closely related to the problem known as A lens .We illustrate this by exhibiting the degeneracy between log| f R 0 | and A lens .This analysis also demonstrates that HS f (R) cannot by itself resolve this tension.
Setting A lens = 1 and adding galaxy power spectra from BOSS and their cross-correlation with CMB lensing, we then constrain log| f R 0 |< −4.61 at 95% confidence (Sect.5.2).This is our central result.It means that while HS f (R) may still explain the accelerated expansion, there is no signature of the model in current observations of large-scale structure; the model predictions do not substantially deviate from those of General Relativity.
We also showed that the cross-correlation of galaxy and lensing measurements is essential in breaking the degeneracy be-tween galaxy bias and log| f R 0 |.Its addition improves the constraint on | f R 0 | by more than an order of magnitude (Sect.5.2.1).This paper is the first to make use of the cross-correlation between CMB lensing and galaxy measurements to constrain HS f (R) gravity.It paves the way for future large-scale galaxy surveys to place more stringent constraints, benefiting from lower noise and using galaxy lensing in addition to CMB lensing.

Fig. 3 .
Fig.3.Marginalized posterior of log| f R 0 | using observations of Planck modeled with MGCLASS.TTTEEE refers to the temperature and polarization power spectra of Planck, "lensing" to the addition of the CMB lensing power spectrum, and A lens to the inclusion of the lensing amplitude (see the text for details).

Fig. 8 .
Fig.8.Marginalized constraints on log| f R 0 | and the galaxy bias of CMASS, when using CMB observations and the full 3×2 pt observables (blue) or only the galaxy and CMB lensing auto power spectra (red).It can be seen that adding the galaxy -CMB lensing cross-correlation reduces the degeneracy between the two parameters.Those contours correspond to the 68 th and 95 th percentiles of the posterior samples.

Fig. 9 .
Fig.9.Marginalized constraints on log| f R 0 | and A lens , where A lens is considered as a systematic effect and is only applied to the CMB temperature and polarization anisotropies.Here, the 3 × 2 pt observables are included, but are not affected by A lens .A value of A lens > 1 is still preferred.The contours correspond to the 68 th and 95 th percentiles of the posterior samples.

Table 2 .
Constraints on log| f R 0 | obtained in the different cases.CMB here refers to the temperature and polarization power spectra.The CMB alone, which is not included within this table, gives an incongruous constraint, reflecting the A lens tension, and is discussed in the text in Sect.5.1.