© The Authors 2025

1. Introduction

In 1998, astronomers made the surprising discovery that the expansion of the Universe is accelerating, not slowing down (Riess et al. 1998; Perlmutter et al. 1999). This late-time acceleration of the Universe has become the most challenging problem in theoretical physics. The main aim of ongoing and future cosmological surveys is to address the key questions about the origin of the late-time acceleration. The acceleration can be driven by a cosmological constant or dark energy that evolves with the expansion of the Universe, Alternatively, there could be no dark energy if general relativity (GR) itself is in error on cosmological scales. There has been significant progress in developing modified theories of gravity and these have been developed into tests of GR itself via cosmological observations (Koyama 2016; Ishak 2019); Ferreira:2019xrr. Testing gravity is one of the main objectives of stage-IV dark energy surveys (Albrecht et al. 2006).

Of particular importance in these surveys is the Euclid mission (Euclid Collaboration: Mellier et al. 2025). The Euclid satellite undertakes a spectroscopic survey of galaxies and an imaging survey (targeting weak lensing (WL) which can also be used to reconstruct galaxy clustering using photometric redshifts). The combination of these two probes is essential for cosmological tests of gravity (Euclid Collaboration: Blanchard et al. 2020, EC20 hereafter).

Modified gravity models typically include an additional scalar degree of freedom that gives rise to a fifth force. To satisfy the stringent constraints on deviations from GR in the Solar System, many modified gravity models utilise screening mechanisms to hide modifications of gravity on small scales (Joyce et al. 2015; Brax et al. 2021). This is accomplished by nonlinearity in the equation that governs the dynamics of the scalar degree of freedom. This significantly complicates the nonlinear modelling of matter clustering in these models as the nonlinear equation for the scalar mode coupled to nonlinear density fields needs to be solved. A systematic comparison of N-body simulations in f(R) gravity and normal-branch Dvali–Gabadadze–Porratti (nDGP) models was done in Winther et al. (2015). Since then, new simulations have been developed, including those using approximate methods to accelerate simulations (Valogiannis & Bean 2017; Winther et al. 2017). A fitting formula (Winther et al. 2019) and emulators for the nonlinear matter power spectrum have been developed (Arnold et al. 2019a; Ramachandra et al. 2021; Sáez-Casares et al. 2023; Fiorini et al. 2023). At the same time, a semi-analytic method based on the halo model to predict the nonlinear matter power spectrum for general dark energy and modified gravity models has been developed (Cataneo et al. 2019; Bose et al. 2021, 2023; Carrilho et al. 2022). These nonlinear predictions were used to study modifications of the WL observables (Schneider et al. 2020; Harnois-Déraps et al. 2023; Spurio Mancini & Bose 2023; Carrion et al. 2024; Tsedrik et al. 2024), cross-correlation of galaxies, and cosmic microwave background (Kou et al. 2024), for example in f(R) gravity.

In Casas et al. (2023), Fisher Matrix forecasts were performed to predict Euclid’s ability to constrain f(R) gravity models. In the Hu–Sawicki f(R) model, it was shown that in the optimistic setting defined in EC20, and for a fiducial value of |f_R0| = 5 × 10⁻⁶, Euclid alone will be able to constrain the additional parameter log₁₀|f_R0| at the 3% level, using spectroscopic galaxy clustering alone; at the 1.4% level, using the combination of photometric probes on their own; and at the 1% level, using the combination of spectroscopic and photometric observations. The forecast for photometric probes used the fitting formula for the nonlinear matter power spectrum obtained in Winther et al. (2019). Further Fisher Matrix forecasts have been done for other modified gravity models with scale-independent linear growth (Frusciante 2024).

For real data analysis, it is imperative to check the effect of the accuracy of the theoretical modelling. This article is part of a series that collectively explores simulations and nonlinearities beyond the Λ cold dark matter (ΛCDM) model:

1. Numerical methods and validation (Euclid Collaboration: Adamek et al. 2025, Paper 1 hereafter).

2. Results from non-standard simulations (Euclid Collaboration: Racz et al. 2025, Paper 2 hereafter).

3. Constraints on f(R) models from the photometric primary probes (this work).

4. Cosmological constraints on non-standard cosmologies from simulated Euclid probes (D’Amico et al. in prep.).

In Paper 1, the comparison of N-body simulations performed in Winther et al. (2015) in the Hu–Sawicki f(R) and nDGP models was extended to add more simulations. The comparison was done for the matter power spectrum and the halo mass function. The measurements of these quantities were done using the dedicated pipeline developed by Paper 2. In Paper 2, additional simulations have been analysed in addition to those used in Paper 1. In this paper, we utilise these developments and compare several predictions for the nonlinear matter power spectrum in f(R) gravity. We perform Markov chain Monte Carlo (MCMC) simulations using synthetic data for Euclid photometric probes, and assess the impact of using different nonlinear models for the matter power spectrum at the level of parameter constraints. In addition, we add baryonic effects using a baryonification method that was not included in the Fisher Matrix forecast. Although this paper focuses on f(R) gravity, for which multiple public codes are available to predict the nonlinear matter power spectrum, the methodology and code developed in this paper are readily applicable to other modified gravity models. The validation of the nonlinear models for spectroscopic probes has been done in Euclid Collaboration: Bose et al. (2024) and D’Amico et al. (Paper 4) will perform a similar analysis to this paper’s for the spectroscopic probes.

This paper is organised as follows. In Sect. 2, we summarise theoretical predictions for Euclid observables based on EC20 and Casas et al. (2023). In Sect. 3, we introduce the Hu–Sawicki f(R) gravity model and summarise the four nonlinear models for the matter power spectrum. We then discuss the implementation of these models in the MontePython code introduced in Casas et al. (2024). In Sect. 4, we compare predictions for the nonlinear matter power spectrum with N-body simulations and study their impact on the angular power spectra for Euclid photometric probes. Forecasts for errors from the combination of photometric probes considered in EC20 are presented based on the synthetic data created by four different nonlinear models. We also compare the result with the Fisher Matrix forecast. We then study the bias in the recovered f(R) gravity parameter when the synthetic data are created by a different nonlinear model. Section 5 is devoted to the study of baryonic effects using the BCemu baryonification method. We show how the bias is affected by baryonic effects and obtain the upper bound on |f_R0| using the ΛCDM synthetic data. In Sect. 6, we implement theoretical errors to take into account the uncertainties of theoretical predictions. We conclude in Sect. 7.

2. Theoretical predictions for Euclid observables

In this section, we discuss how moving away from the standard GR assumption impacts the predictions for the angular power spectra C(ℓ) that will be compared with the photometric data survey. The observables that need to be computed and compared with the data are the angular power spectra for weak lensing (WL), photometric galaxy clustering (GC_ph) and their cross-correlation (XC_ph). In EC20, these were calculated using the Limber and flat-sky approximations in a flat ΛCDM Universe, as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\mathit{XY}}(\ell )=c{\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}z\frac{W_i^X(z)W_j^Y(z)}{H(z)r^2(z)}{P}_{\delta \delta }(k_{\ell },z),\end{array}$$(1)

where W_i^X(z) is the window function in each tomographic redshift bin, i, with X = {L, G} corresponding to WL and GC_ph, respectively. Here, k_ℓ = (ℓ+1/2)/r(z), r(z) is the comoving distance to redshift, z, and H(z) is the Hubble function. The nonlinear power spectrum of matter density fluctuations, P_δδ(k_ℓ, z), is evaluated at wave number k_ℓ and redshift z, in the redshift range of the integral from z_min = 0.001 to z_max = 4.

However, when abandoning the assumption of GR, one has to account for changes in the evolution of both the homogeneous background and cosmological perturbations. For WL, it is the Weyl potential, ϕ_W, that determines the angular power spectrum. The power spectrum of the Weyl potential is related to P_δδ as (Casas et al. 2023)

$$\begin{array}{c}\hfill P_{{\varphi }_W}(k,z)={[-3{\mathrm{\Omega }}_{\mathrm{m}}{\left(\frac{H_0}{c}\right)}^2(1+z)\mathrm{\Sigma }(k,z)]}^2{P}_{\delta \delta }(k,z),\end{array}$$(2)

where Σ(k, z) is a phenomenological parameterisation to account for deviations from the standard ΛCDM lensing prediction. Ω_m is the density parameter of matter and H₀ is the Hubble parameter – both at the present time. Here, we assumed a standard background evolution of the matter component, that is, ρ_m(z) = ρ_m, 0(1 + z)³.

We can therefore use the recipe of Eq. (1) with H, r and P_δδ provided by a Boltzmann solver, but with the new window functions (Spurio Mancini et al. 2019)

$$\begin{array}{cc}\hfill W_i^{\mathrm{G}}(k,z)=& \frac{1}{c}{b}_i(k,z)n_i(z)H(z),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill W_i^{\mathrm{L}}(k,z)=& \frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m}}{\left(\frac{H_0}{c}\right)}^2(1+z)r(z)\mathrm{\Sigma }(k,z)\hfill \\ \hfill & \times {\displaystyle {\int }_z^{z_{\max }}}\mathrm{d}z\prime n_i(z)\frac{r(z\prime )-r(z)}{r(z\prime )}+W_i^{\mathrm{IA}}(k,z),\hfill \end{array}$$(4)

where n_i(z) is the normalised galaxy number-density distribution in a tomographic redshift bin, i, such that ∫_zmin^z_maxn_i(z)dz = 1, and W_i^IA(k, z) encodes the contribution of intrinsic alignments (IA) to the WL power spectrum. We follow EC20 in assuming an effective scale-independent galaxy bias, constant within each redshift bin, and its values, b_i, are introduced as nuisance parameters in our analysis, with their fiducial values determined by $b_i=\sqrt{1+{\overline{z}}_i}$, where ${\overline{z}}_i$ is the mean redshift of each redshift bin.

The IA contribution is computed following the nonlinear alignment model with a redshift dependent amplitude (EC20), in which

$$\begin{array}{c}\hfill W_i^{\mathrm{IA}}(k,z)=-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m}}{\mathcal{F}}_{\mathrm{IA}}(z)}{\delta (k,z)/\delta (k,0)}{n}_i(z)\frac{H(z)}{c},\end{array}$$(5)

where

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{IA}}(z)={(1+z)}^{{\eta }_{\mathrm{IA}}}.\end{array}$$(6)

The parameters 𝒜_IA and η_IA are the nuisance parameters of the model, and 𝒞_IA is a constant accounting for dimensional units. The galaxy distribution is binned into 10 equipopulated redshift bins with an overall distribution following

$$\begin{array}{c}\hfill n(z)\propto {\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{3/2}],\end{array}$$(7)

with $z_0=0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is ${\overline{n}}_{\mathrm{g}}=30{\mathrm{arcmin}}^{-2}$ (EC20).

Changes in the theory of gravity impact the IA contribution, introducing a scale dependence through the modified perturbations’ growth. This is explicitly taken into account in Eq. (5) through the matter perturbation δ(k, z), while the modifications on the clustering of matter in the GC_ph case are accounted for in the new P_δδ(k_ℓ, z).

Finally, we present the likelihood that we used. We followed the approach presented in Casas et al. (2024). We first constructed an (N^G + N^L)×(N^G + N^L) angular power spectrum matrix for each multipole, where the different N correspond to the number of redshift bins for each probe (WL and GC_ph):

$$\begin{array}{c}\hfill {\mathsf{C}}_{\ell }=\left[\begin{array}{cc}{C}_{\mathit{ij}}^{\mathrm{LL}}(\ell )& C_{\mathit{ij}}^{\mathrm{GL}}(\ell )\\ [3pt]C_{\mathit{ij}}^{\mathrm{LG}}(\ell )& C_{\mathit{ij}}^{\mathrm{GG}}(\ell )\end{array}\right],\end{array}$$(8)

where lower-case Latin indexes i, … run over all tomographic bins. Similarly, the noise contributions can also be condensed into one noise matrix, N:

$$\begin{array}{c}\hfill {\mathsf{N}}_{\ell }=\left[\begin{array}{cc}{N}_{\mathit{ij}}^{\mathrm{LL}}(\ell )& N_{\mathit{ij}}^{\mathrm{GL}}(\ell )\\ [3pt]N_{\mathit{ij}}^{\mathrm{LG}}(\ell )& N_{\mathit{ij}}^{\mathrm{GG}}(\ell )\end{array}\right],\end{array}$$(9)

where the noise terms, $N_{\mathit{ij}}^{\mathit{AB}}(\ell )$, are given by

$$\begin{array}{cc}\hfill N_{\mathit{ij}}^{\mathrm{LL}}(\ell )& =\frac{{\delta }_{\mathit{ij}}^{\mathrm{K}}}{{\overline{n}}_i}{\sigma }_{ϵ}^2,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill N_{\mathit{ij}}^{\mathrm{GG}}(\ell )& =\frac{{\delta }_{\mathit{ij}}^{\mathrm{K}}}{{\overline{n}}_i},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill N_{\mathit{ij}}^{\mathrm{GL}}(\ell )& =0,\hfill \end{array}$$(12)

where σ_ϵ² = 0.3² is the variance of observed ellipticities. We can further define ${\widehat{\mathsf{C}}}_{\ell }:={\mathsf{C}}_{\ell }+{\mathsf{N}}_{\ell }$. This is the covariance of the spherical multipole moments, a_lm.

In the Gaussian approximation, it can be shown that the description in EC20, using the covariance of the angular power spectra, is equivalent to the description found in Casas et al. (2024), using the covariance of the a_ℓm. The likelihood then can be expressed in terms of the observed covariance

$$\begin{array}{c}\hfill {\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{obs}}(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{\ell }}{\left(a_{\ell m}\right)}_i{\left(a_{\ell m}\right)}_j^{\ast },\end{array}$$(13)

and the theoretically predicted one, ${\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{th}}(\ell )$, as

$$\begin{array}{c}\hfill {\chi }^2=f^{\mathrm{sky}}{\displaystyle \sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}}(2\ell +1)[\frac{d^{\mathrm{mix}}}{d^{\mathrm{th}}}+ln\left(\frac{d^{\mathrm{th}}}{d^{\mathrm{obs}}}\right)-N],\end{array}$$(14)

where the determinants, d, are defined as

$$\begin{array}{cc}\hfill d^{\mathrm{th}}(\ell )& =det\left[{\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{th}}(\ell )\right],\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill d^{\mathrm{obs}}(\ell )& =det\left[{\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{obs}}(\ell )\right],\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill d^{\mathrm{mix}}(\ell )& ={\displaystyle \sum _k^N}det\left[\{\begin{array}{cc}{\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{obs}}(\ell )& \text{for}k=j\\ {\widehat{\mathrm{C}}}_{\mathit{ij}}^{\mathrm{th}}(\ell )& \text{for}k\ne j\end{array}\right].\hfill \end{array}$$(17)

Here, N is the number of bins, and thus either (N^G + N^L) for multipoles for which we include the cross correlation, or the respective N for multipoles for which we treat the two probes separately. The additional factor, f^sky, comes from an approximation to account for having fewer available independent ℓ modes due to partial sky coverage. In this paper, we set the observed covariance ${\widehat{\mathrm{C}}}^{\mathrm{obs}}$ to the theoretically predicted one computed at the fiducial cosmology. Following EC20 and Casas et al. (2024), we do not include super-sample covariance (SSC) in this work. The SSC impact was shown to be non-negligible and will need to be included in future analyses as shown in Euclid Collaboration: Sciotti et al. (2024).

We consider two different scenarios: an optimistic and a pessimistic case. In the optimistic case, we consider ℓ_max = 5000 for WL, and ℓ_max = 3000 for GC_ph and XC_ph. Instead, in the pessimistic scenario, we consider ℓ_max = 1500 for WL, and ℓ_max = 750 for GC_ph and XC_ph.

In the smallest photometric redshift bin, the galaxy number density distribution, n(z), peaks at around z = 0.25. Under the Limber approximation for our fiducial cosmology, the corresponding maximum values of k evaluated in the power spectrum corresponding to the pessimistic and optimistic scenario for GC_ph are $k_{\max }=[0.7,2.9]h{\mathrm{Mpc}}^{-1}$, respectively, while for the WL, maximum wavenumbers probed are $k_{\max }=[1.4,4.8]h{\mathrm{Mpc}}^{-1}$ at the peak redshift z = 0.25 (Casas et al. 2023). Here, h denotes the dimensionless Hubble parameter h ≔ H₀/(100 km s⁻¹ Mpc⁻¹). For smaller values of z, the values of k at a given ℓ increase, but the window functions in Eqs. (3) and (4) suppress the power spectrum and we set it to zero after a fixed $k_{\max }=30h{\mathrm{Mpc}}^{-1}$. We list the specific choices of scales and settings used for each observable in Table 1. Although the currently foreseen specification of the Euclid Wide Survey differs from that assumed in EC20, for example in terms of the survey area, we use their results to allow for comparison with earlier forecasts.

3. MontePython implementation

In this section, we describe the implementation of the likelihood discussed in Sect. 2 in the MontePython code developed in Casas et al. (2024), in the Hu–Sawicki f(R) gravity model.

3.1. Hu–Sawicki f(R) gravity

Modified gravity f(R) models (Buchdahl 1970) are models constructed by promoting the Ricci scalar, R, in the Einstein–Hilbert action to a generic function of R → R + f(R), that is,

$$\begin{array}{c}\hfill S=\frac{c^4}{16\pi G}{\displaystyle \int {\mathrm{d}}^{4}x\sqrt{-g}[R+f(R)]}+S_{\mathrm{m}}[g_{\mu \nu },{\mathrm{\Psi }}_{\mathrm{m}}],\end{array}$$(18)

where g_μν is the metric tensor, g its determinant, and S_m represents the matter sector with its matter fields, Ψ_m. For further discussions, we refer to Sotiriou & Faraoni (2010), Clifton et al. (2012), and Koyama (2016).

How this modification changes gravity is more easily seen by formulating the theory in the so-called Einstein frame by performing a conformal transformation $g_{\mu \nu }={\stackrel{\sim }{g}}_{\mu \nu }{A}^2(\varphi )$ where $A(\varphi )=\sqrt{1+f_R}=e^{\varphi /\sqrt{6}}$, to obtain

$$\begin{array}{cc}\hfill S=& \frac{c^4}{16\pi G}{\displaystyle \int {\mathrm{d}}^{4}x\sqrt{-\stackrel{\sim }{g}}[\stackrel{\sim }{R}+\frac{1}{2}{(\partial \varphi )}^2-V(\varphi )]}\hfill \\ \hfill & +S_{\mathrm{m}}[A^2(\varphi ){\stackrel{\sim }{g}}_{\mu \nu },{\mathrm{\Psi }}_{\mathrm{m}}],\hfill \end{array}$$(19)

with the potential V = [f_RR − f(R)]/2(1 + f_R)² and f_R ≔ df(R)/dR. This demonstrates that the theory reduces to standard GR together with an extra scalar field that is coupled to the matter sector giving rise to a fifth force. This fifth force has a finite range, λ_C, and within this range, it mediates a force that has a strength that, in the linear regime, corresponds to 1/3 of the usual gravitational force:

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{fifth}}=\frac{1}{3}\frac{\mathit{GMm}}{r^2}(1+r/{\lambda }_{\mathrm{C}}){\mathrm{e}}^{-r/{\lambda }_{\mathrm{C}}}.\end{array}$$(20)

This effectively changes $G\to \frac{4}{3}G$ on small scales (r ≪ λ_C) while keeping the usual G on large scales. Such a large modification would be ruled out by observations if it was not for the fact that the theory possesses a screening mechanism (Khoury & Weltman 2004; Brax et al. 2008) which hides the modifications in high-density regions. This screening mechanism effectively suppresses the fifth force by a factor ∝f_R/Φ_N, where Φ_N is the Newtonian gravitational potential.

Not all f(R) models one can write down possess such a screening mechanism, which places some restrictions on their functional form (see e.g., Sotiriou & Faraoni 2010). One concrete example of a model that has all the right ingredients is the model proposed by Hu & Sawicki (2007), which in the large-curvature limit is defined by

$$\begin{array}{c}\hfill f(R)=-6{\mathrm{\Omega }}_{\mathrm{DE}}\frac{H_0^2}{c^2}+|f_{R0}|\frac{{\overline{R}}_0^2}{R},\end{array}$$(21)

where

$$\begin{array}{c}\hfill {\overline{R}}_0=3{\mathrm{\Omega }}_{\mathrm{m}}\frac{H_0^2}{c^2}(1+4\frac{{\mathrm{\Omega }}_{\mathrm{DE}}}{{\mathrm{\Omega }}_{\mathrm{m}}})\end{array}$$(22)

is the Ricci scalar in the cosmological background and Ω_DE is the density parameter of dark energy at the present time. The first term in Eq. (21) corresponds to a cosmological constant and the only free parameter is |f_R0|. In the limit |f_R0|→0, we recover GR and the ΛCDM model. This parameter controls the range of the fifth force and, in the cosmological background, we have ${\lambda }_{\text{C}}\simeq 32\sqrt{|f_{R0}|/{10}^{-4}}\mathrm{Mpc}$ at the present time. Solar System constraints require |f_R0|≲10⁻⁶, cosmological constraints currently lie around 10⁻⁶–10⁻⁴ depending on the probe in question (see e.g., Fig. 28 in Koyama 2016, for a summary) while astrophysical constraints at the galaxy scale can be as tight as ≲10⁻⁸ (Desmond & Ferreira 2020), but see Burrage et al. (2024) for a recent note of caution on such galactic-scale constraints.

The energy density of the scalar field contributes in general to the expansion of the Universe; however, for viable models, like the one we consider here, this contribution is tiny (of the order |f_R0| Ω_DE, 0) apart from the constant part of the potential, which is indistinguishable from a cosmological constant. The background evolution of such models is therefore very close to ΛCDM. Since f(R) models have a conformal coupling, light deflection is weakly affected as follows

$$\begin{array}{c}\hfill \mathrm{\Sigma }(z)=\frac{1}{1+f_R(z)}·\end{array}$$(23)

Since the maximum value of |f_R(z)| is given by |f_R0|, for the values of |f_R0| we consider in this paper, we can ignore this effect. Thus gravitational lensing is also not modified in the sense that the lensing potential is sourced by matter in the same way as in GR (though the underlying density field will of course be different). The main cosmological signatures of the model therefore come from having a fifth force, acting only on small scales r ≲ λ_C, in the process of structure formation. The main effect of the screening mechanism is that the prediction for the amount of clustering will generally be much smaller than what naive linear perturbation theory predicts.

3.2. Nonlinear modelling

We implement Ξ(k, z), defined as

$$\begin{array}{c}\hfill \mathrm{\Xi }(k,z)\frac{P_{f(R)}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)},\end{array}$$(24)

to obtain the nonlinear f(R) matter power spectrum. For the ΛCDM power spectrum P_ΛCDM(k, z), we use the halofit ‘Takahashi’ prescription (Takahashi et al. 2012) following Casas et al. (2023). It includes the minimum mass for massive neutrinos in P_ΛCDM(k, z), but ignores the effect of massive neutrinos on Ξ(k, z). This approximation was shown to be well justified for the minimum mass of neutrinos in Winther et al. (2019) using data from Baldi et al. (2014). We describe below four models for Ξ(k, z) used in this paper (see Table 2). We note that we can use any ΛCDM power spectrum prediction in our approach such as EuclidEmulator2 (Euclid Collaboration: Knabenhans et al. 2021) and Bacco (Angulo et al. 2021). The exercise we perform in this paper is a comparison of the different prescriptions for Ξ(k, z). Given this, the ΛCDM nonlinear spectrum prescription does not matter and it is common to all nonlinear models.

When we apply our pipeline to observational data, it is important to control the accuracy of the ΛCDM nonlinear spectrum prescription. We need to implement scale cuts to account for the inaccuracy of the ΛCDM prediction. The validation of the ΛCDM nonlinear prescription is still underway in the Euclid consortium, and we leave the investigation of the effect of using different prescriptions in ΛCDM for future work.

3.2.1. Fitting formula

A fitting formula for Ξ(k, z) was developed in Winther et al. (2019) and describes the enhancement in the power spectrum compared to a ΛCDM nonlinear power spectrum as a function of the parameter |f_R0|. This fitting function has been calibrated using the DUSTGRAIN (Giocoli et al. 2018) N-body simulations run by MG-Gadget (Puchwein et al. 2013) and the ELEPHANTN-body simulations (Cautun et al. 2018) run by ECOSMOG (Li et al. 2012; Bose et al. 2017). The main approximation is that the cosmological parameter dependence of Ξ(k, z) is ignored. In Winther et al. (2019), this assumption was checked using simulations with different σ₈, Ω_m, as well as the mass of massive neutrinos. Winther et al. (2019) also corrected the fitting formula to account for additional dependence on these parameters. In this paper, we do not include these corrections as the previous forecast paper (Casas et al. 2023) used the fitting formula without these corrections.

The fitting function has in total 54 parameters for the full scale, redshift and |f_R0| dependence. This fitting formula is not defined outside the range 10⁻⁷ < |f_R0| < 10⁻⁴. The code is publicly available¹.

3.2.2. Halo model reaction

We further consider the halo model reaction approach of Cataneo et al. (2019) which combines the halo model and perturbation theory frameworks to model corrections coming from non-standard physics. The nonlinear power spectrum is given by

$$\begin{array}{c}\hfill P_{\mathrm{NL}}(k,z)=\mathcal{R}(k,z)P_{\mathrm{NL}}^{\mathrm{pseudo}}(k,z),\end{array}$$(25)

where $P_{\mathrm{NL}}^{\mathrm{pseudo}}(k,z)$ is called the pseudo-power spectrum and is defined as a nonlinear ΛCDM spectrum with initial conditions tuned such that the linear clustering at the target redshift z matches the beyond-ΛCDM case. This choice ensures the halo mass functions in the beyond-ΛCDM and pseudo-universes are similar, giving a smoother transition of the power spectrum over inter- and intra-halo scales. We model the pseudo-cosmology nonlinear power spectrum using HMCode2020 (Mead et al. 2015, 2016, 2021) by supplying the code with the linear f(R) power spectrum.

The halo model reaction, ℛ(k, z), models all the corrections to the pseudo spectrum coming from nonlinear beyond-ΛCDM physics. We refer the reader to Cataneo et al. (2019), Bose et al. (2021) and Frusciante (2024) for the exact expressions for this term. The halo model reaction can be computed efficiently using the publicly available ReACT² (Bose et al. 2020, 2021, 2023) code.

Despite ReACT being highly efficient, having been used in previous real cosmic shear analyses (Tröster et al. 2021), it is still too computationally expensive for the number of tests we wish to perform. To accelerate our inference pipeline, we created a neural network emulator using the Cosmopower package³ (Spurio Mancini et al. 2022) for the halo model reaction-based boost

$$\begin{array}{c}\hfill \mathrm{\Xi }(k,z)=\frac{\mathcal{R}(k,z)P_{\mathrm{HMCode}2020}^{\mathrm{pseudo}}(k,z)}{P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z)},\end{array}$$(26)

where P_ΛCDM(k, z) and $P_{\mathrm{HMCode}2020}^{\mathrm{pseudo}}(k,z)$ were calculated using HMCode2020⁴ (Mead et al. 2021) while ℛ(k, z) was calculated using ReACT. We chose HMCode2020 to model the pseudo-power spectrum as it has been shown to have improved accuracy and does not show suppression of power with respect to ΛCDM (Ξ < 1) at high redshifts, which is not expected.

To train the emulator, we followed the procedure of Spurio Mancini & Bose (2023) but widened the parameter priors. We produce ∼10⁵ boost predictions in the range k ∈ [0.01, 3] h Mpc⁻¹ and z ∈ [0, 2.5]. We take cosmologies from the Latin hypercube given by the ranges in Table 3 and Table 4, with the massive neutrino density parameter today’s range being Ω_ν ∈ [0.0, 0.00317]. Emulation of the boost speeds up the computation by 4 orders of magnitude and we find similar accuracy of our emulator as found in Spurio Mancini & Bose (2023). The emulator is publicly available⁵. Finally, we note another small difference between our emulator and that of Spurio Mancini & Bose (2023). In this work, we assume P_ΛCDM in Eq. (26) with Ω_m equal to the total of the true cosmology, whereas in Spurio Mancini & Bose (2023) they assume that Ω_m = Ω_b + Ω_c, Ω_b and Ω_c being, respectively the baryon and cold dark matter density parameters today. The emulators give the same output for Ω_ν = 0, which is what we assume in this work for Ξ.

3.2.3. FORGE

The FORGE emulator was introduced in Arnold et al. (2022). It is based on simulations for 50 combinations of |f_R0|, Ω_m, σ₈^ΛCDM and h with all other parameters fixed. We note that σ₈^ΛCDM is the σ₈ we would obtain in a ΛCDM model with the same cosmological parameters and initial amplitude A_s and not σ₈ in an f(R) gravity model. The emulator accuracy is better than 2.5% around the centre of the explored parameter space, up to scales of k = 10 h Mpc⁻¹. f(R) simulations are performed by a modified version of the Arepo code, MG-AREPO (Springel 2010; Arnold et al. 2019a; Weinberger et al. 2020) that solves the nonlinear scalar field equation using a relaxation method. (see Paper 1 for more details.)

The emulation was made for the ratio between the power spectrum in f(R) and the halofit prediction in ΛCDM. We note that this is different from Ξ(k, z) as the power spectrum in ΛCDM models in these simulations can have slight deviations from the halofit prediction. To account for this effect, the authors provided the ratio of the power spectrum in a reference ΛCDM model to the halofit prediction. This ΛCDM simulation uses Ω_m = 0.31315, h = 0.6737, σ₈^ΛCDM = 0.82172. This can be used to obtain Ξ(k, z). However, this correction is provided only in this ΛCDM model. Thus, the assumption here is that this correction is independent of cosmological parameters. The latest version of the FORGE simulations has ΛCDM counterparts to the f(R) simulations using the same seed. These simulations were analysed in Paper 2. Hence, it is in principle possible to emulate directly Ξ(k, z) from these simulations. However, as this is not publicly available, we opted for using the original FORGE emulator⁶ for this paper. We use one of these simulations for validation in Sect. 4.

3.2.4. e-Mantis

We also consider the predictions given by the e-Mantis emulator presented in Sáez-Casares et al. (2023), which can predict the ratio Ξ(k, z) between the nonlinear matter power spectrum in f(R) gravity and in ΛCDM. The predictions are calibrated from a large suite of N-body simulations run with the code ECOSMOG (Li et al. 2012; Bose et al. 2017), a modified gravity version of the Adaptive Mesh Refinement (AMR) N-body code RAMSES (Teyssier 2002; Guillet & Teyssier 2011). The e-Mantis simulation suite covers the $\{|f_{R_0}|,{\mathrm{\Omega }}_{\text{m}},{\sigma }_8^{\mathrm{\Lambda }\text{CDM}}\}$ parameter space with 110 cosmological models sampled from a Latin hypercube (McKay et al. 1979). The remaining ΛCDM parameters, h, n_s and, Ω_b, remain fixed, which means that their impact on Ξ(k, z) is ignored. In Sáez-Casares et al. (2023), it was shown that the error made by this approximation is at the sub-percent level. The quantity Ξ(k, z) is measured from pairs of f(R) and ΛCDM simulations, both run with the same initial conditions, which leads to a large cancellation of cosmic variance and numerical resolution errors. The emulation is done through a Gaussian Process regression (see e.g. Rasmussen & Williams 2005). The emulator can give predictions for redshifts z ∈ [0,2], wavenumbers k ∈ [0.03,10] h Mpc⁻¹ and cosmological parameters in the following range: |f_R0| ∈ [10⁻⁷,10⁻⁴], Ω_m ∈ [0.2365,0.3941] and σ₈^ΛCDM ∈ [0.6083,1.0140]. The estimated accuracy of e-Mantis, including emulation errors and systematic errors in the training data, is better than 3% for scales k ≲ 7 h Mpc⁻¹ and across the whole parameter space although in most cases, the accuracy is of order 1%. The code is publicly available⁷.

3.3. MontePython implementation

Our implementation of the likelihood is based on the MontePython implementation described in Casas et al. (2024). For a more detailed explanation of the numerical implementation, we refer to that work. For our purposes, we have modified this code to include our prediction for Ξ(k, z) to obtain the nonlinear power spectrum in the f(R) gravity model. As the different implementations have different ranges in cosmological parameters, wavenumbers, and redshifts, we have chosen an extrapolation scheme to unify the ranges. We note that the exact implementation does not affect the final results strongly. We checked that the contributions from high redshifts and high wavenumbers which used extrapolations were sub-dominant.

For the modified gravity parameter |f_R0|, we used flat priors in terms of log₁₀|f_R0| to stay within the validity range of each model. For the largest scales beyond the range of the emulators, we set Ξ = 1. This is because the effect of the fifth force has to vanish on these scales. For small scales, we did a power-law extrapolation. To obtain the spectral index for the extrapolation, we proceeded as follows: we calculated Ξ on a grid close to the edges of the emulators; we then fitted, for fixed z, a linear function in the log-log space onto this grid. This was done to average out any numerical noise at the edge. We did the same for high redshifts by fitting the power law for fixed wavenumbers. For regions of both high k and z, we did a constant extrapolation of the spectral index. We checked the dependence of the extrapolation methods (the power law and constant extrapolation as well as the ΛCDM-extrapolation to Ξ = 1) at the level of the angular power spectra and confirmed that the effect was smaller than the errors on the synthetic data. To keep the extrapolated function Ξ(k, z) from going to non-physical values, we set a hard lower bound of Ξ = 1 and an upper bound of Ξ = 2. The ranges of the different emulators are found in Table 3. We did a constant extrapolation in the ΛCDM parameters. This can be done because during a typical MCMC the majority of the suggested points are within the ranges of our emulators. In addition, Ξ is insensitive to cosmological parameters as shown by Winther et al. (2019). The validity range of ΛCDM cosmological parameters is summarised in Table 4.

We also added the effect of baryonic feedback in the form of an additional correction, Ξ^BFM, to our power spectrum prediction. The physics of the baryonic feedback effects is discussed in Sect. 5. We refer the reader to Schneider et al. (2020) and Mead et al. (2021) for further details. We obtained the correction from BCemu by Giri & Schneider (2021). This emulator is only trained for redshifts below z = 2 and wavenumbers k < 12.5 h Mpc⁻¹. In this case, we use constant extrapolation for both k and z. We checked that this had little effect on our results. This is because the Euclid main probes are mostly sensitive to redshifts around z ∼ 1. For these redshifts, the extrapolation only affects very high multipoles, ℓ ≳ 4500. Thus, the extrapolation has negligible effect on the angular power spectrum. We estimate the overall effect this has on our result to be at most at the percent level. In the absence of modified gravity, to obtain the nonlinear baryonic feedback power spectrum, we multiply the ΛCDM nonlinear power spectrum by Ξ^BFM. When adding the effect of modified gravity, we combine both boosts as

$$\begin{array}{c}\hfill P_{f(R)}^{\mathrm{BFM}}(k,z)={\mathrm{\Xi }}^{\mathrm{BFM}}(k,z)\mathrm{\Xi }(k,z)P_{\mathrm{\Lambda }\mathrm{CDM}}(k,z).\end{array}$$

We can make this approximation if both effects are independent. This was shown to be the case for small deviations from ΛCDM in f(R) gravity models considered in this paper using hydrodynamical simulations by Arnold et al. (2019a) and Arnold & Li (2019).

Our final addition to the code is the inclusion of theoretical errors. For this, we have adjusted the prescription of Audren et al. (2013). The theory and numerical implementation is further discussed in Sect. 6.

4. Comparison of different nonlinear models

In this section, we compare four different predictions for the nonlinear dark matter power spectrum in terms of Ξ(k, z) defined in Eq. (24) introduced in Sect. 3. We start from a comparison of the matter power spectrum with N-body simulations. We then compare the angular power spectra in the Euclid reference cosmology. We perform MCMC analysis to compare errors and investigate bias due to the difference in the prediction of the nonlinear matter power spectrum using the settings defined in EC20.

4.1. Comparison of predictions

4.1.1. Comparison with N-body simulations

Figure 1 shows a comparison of the ratio of the power spectra between f(R) gravity and ΛCDM measured from N-body simulations with the theoretical predictions for Ξ. These simulations use the same initial conditions and the ratio removes the cosmic variance and the effect of mass resolution. The left-hand side plot shows a comparison using the measurements from Paper 1. This is based on the comparison project in Winther et al. (2015). These simulations were run in a ΛCDM cosmology with Ω_m = 0.269, h = 0.704, n_s = 0.966 and σ₈ = 0.801. The simulations have N_p = 512³ particles of mass M_p ≃ 8.756 × 10⁹ h⁻¹ M_⊙ in a box of size B = 250 h⁻¹ Mpc and start at redshift z = 49. We picked a model called F5 with |f_R0| = 10⁻⁵ and showed the result at z = 0.667 that is presented in Paper 1. In the comparison, we included the measurements from MG-AREPO and ECOSMOG, as FORGE is based on MG-AREPO, while e-Mantis is based on ECOSMOG. The prediction of e-Mantis agrees with ECOSMOG very well. On the other hand, the FORGE prediction deviates from MG-AREPO as well as ECOSMOG. We note that the FORGE prediction is corrected using the ΛCDM simulation (Node 0) in the FORGE simulation suite with Ω_m = 0.31315, h = 0.6737, σ₈^ΛCDM = 0.82172 to obtain Ξ. We find that if we use Ω_m = 0.31315 and h = 0.6737 in the FORGE prediction, the agreement with MG-AREPO is much better.

Fig. 1. Left panel: Comparisons between N-body simulations, e-Mantis and FORGE for |fR0| = 10−5. N-body data from simulations run by MG-AREPO and ECOSMOG codes are taken from the latest comparison project presented in Paper 1. We also include the prediction of FORGE where the cosmological parameters in the fiducial ΛCDM simulation, Ωm = 0.31315 and h = 0.6737, are used to make the prediction (FORGE Node 0). Right panel: Comparison of four prescriptions with an N-body simulation run by MG-AREPO. This simulation is one of the simulations used for training to construct the FORGE emulator (Node 13). The measurement of the power spectrum was done in Paper 2.

To further investigate this issue with FORGE, we used the measurement of the power spectrum in one of the nodes in the FORGE simulation suite (Node 13) run by MG-AREPO that is closest to the Euclid reference cosmology that we will use in this paper with non-zero |f_R0|. This simulation has the following cosmological parameters: Ω_m = 0.34671, h = 0.70056 and σ₈ = 0.78191 and we show a comparison at z = 0.652. Both ΛCDM and f(R) simulations with |f_R0| = 10^−4.90056 are available so that we can measure Ξ directly. We note that the pipeline developed in Paper 2 measures the power spectrum only up k = 3 h Mpc⁻¹. In this case, both FORGE and e-Mantis agree with MG-AREPO within 1% up to k = 3 h Mpc⁻¹. The fitting formula and ReACT agree with MG-AREPO within 1% up to k = 1 h Mpc⁻¹. Given this result, the large discrepancy between FORGE and MG-AREPO is likely to be caused by emulation errors as well as calibrations using the halofitΛCDM power spectra predictions.

4.1.2. Comparison in the Euclid reference cosmology

In this paper, we consider the model called HS6 in Casas et al. (2023), which has the following parameters:

$$\begin{array}{cc}\hfill \mathbf{\Theta }& =\{{\mathrm{\Omega }}_{\mathrm{m}},{\mathrm{\Omega }}_{\mathrm{b}},h,n_{\mathrm{s}},{\sigma }_8,{\log }_{10}|f_{R0}|\},\hfill \\ \hfill \mathrm{HS}6:{\mathbf{\Theta }}_{\mathrm{fid}}& =\{0.32,0.05,0.67,0.96,0.853,-5.301\}.\hfill \end{array}$$(27)

The cosmological parameters are the same ones adopted in EC20. Casas et al. (2023) show that this value of |f_R0| = 5 × 10⁻⁶ can be well constrained by the Euclid photometric probes. Also the range of |f_R0| covered by the four models is wider than the errors predicted in Casas et al. (2023). Our fiducial cosmology includes massive neutrinos with a total mass of ∑m_ν = 0.06 eV, but we keep ∑m_ν fixed in the following analysis.

Figure 2 shows a comparison of the power spectrum and the angular power spectrum for WL, GC_ph and their cross-correlation XC_ph. In these plots, we show the ratio to the ΛCDM prediction and error bars from the diagonal part of the covariance matrix. For the power spectrum comparison at z = 0.5, we see that e-Mantis and fitting agree best. This is not surprising as these two predictions are based on N-body simulations ran by the same code ECOSMOG. On the other hand, FORGE overestimates Ξ at k = 0.1 h Mpc⁻¹. This is similar to the deviation we find in the comparison with the N-body simulation from Paper 1 although the deviation is smaller, at the 2% level. This is likely due to the fact that Ω_m in HS6 is closer to Ω_m in the fiducial ΛCDMFORGE simulation (Ω_m = 0.31315). On the other hand, ReACT underestimates Ξ at k < 0.1 h Mpc⁻¹ and even predicts Ξ < 1. As we discussed in the previous section, we enforce Ξ ≥ 1.

Fig. 2. Comparisons between all models considered for the HS6 fiducial cosmology for angular power spectra of WL (top left), galaxy clustering (top right) and their cross-correlation (bottom right) as well as for the matter power spectrum (bottom left). We show these for the WL redshift bin 6 and galaxy clustering bin 2, with the matter power spectrum plotted at z = 0.5, the redshift at which the kernels of both observables peak for the chosen bins. For all angular spectra, the error bands taken from the diagonal of the covariance are also shown.

4.2. Forecasting errors and biases

4.2.1. Forecasting errors for WL

We first compared errors obtained by running MCMC using the synthetic data created by one of the four nonlinear models and fitting it by the same model. In this case, by definition, we recover the input parameters that were used to create the synthetic data. Our interests are constraints on the |f_R0| parameter and cosmological parameters. We first considered the WL-only case. In this case, we imposed a tight Gaussian prior on the spectral index, n_s, taken from the Planck results (Planck Collaboration VI 2020) and on the baryon density parameter, using Big Bang nucleosynthesis constraints (Pisanti et al. 2021):

$$\begin{array}{c}\hfill n_{\mathrm{s}}=0.96\pm 0.004,{\mathrm{\Omega }}_{\mathrm{b}}{h}^2=0.022445\pm 0.00036,\end{array}$$(28)

as we do not expect to obtain strong constraints on these parameters from WL alone. The parameters that are used in the MCMC runs are summarised in Table 5, including their fiducial values and prior ranges. As a convergence criterion, we used a Gelman–Rubin (Gelman & Rubin 1992) value of R − 1 < 0.01 for each individual sampling parameter using MontePython. For post-processing chains, we used GetDist⁸ (Lewis 2019).

Figure 3 shows the 2D contours of the constraints on parameters, and Table 6 and Table 7 summarise constraints on log₁₀|f_R0| for the optimistic and pessimistic settings. The constraints on cosmological parameters are consistent among the four different nonlinear models, while we see some notable differences in the constraints on |f_R0|. In the case of the optimistic setting, e-Mantis gives the tightest constraints, which is also closer to Gaussian. The fitting formula agrees with e-Mantis for small log₁₀|f_R0| but has a weaker constraint for larger log₁₀|f_R0|. Constraints from FORGE agree with the fitting formula for large log₁₀|f_R0|, but give a weaker constraint for small log₁₀|f_R0|. The degeneracy between log₁₀|f_R0|, Ω_c h² and ln(10¹⁰A_s) also presents some notable differences. The degeneracy for larger log₁₀|f_R0| is different for the fitting formula when compared with e-Mantis and FORGE. This could be attributed to the fact that the fitting formula does not include any cosmological parameter dependence in the prediction for Ξ(k, z). We observe similar agreements and disagreements for the pessimistic setting, but the agreement among e-Mantis, FORGE and the fitting formula is better, particularly for small log₁₀|f_R0|. ReACT gives a weaker constraint on log₁₀|f_R0|, but the constraints on cosmological parameters are consistent among the four different models. This is due to the weak cosmology dependence of Ξ(k, z), as discussed in Winther et al. (2019), and the fact that the constraint on cosmological parameters is coming from the ΛCDM power spectrum, which is common in all these four models.

Fig. 3. Constraints on parameters where the same model is used to create the synthetic data and perform the fitting such that the input parameters are reproduced. Left panel: WL optimistic case. Right panel: WL pessimistic case.

We compared the constraints from the MCMC analysis with the Fisher Matrix forecast and show this in Fig. 3. The MontePython pipeline can be used as a Fisher Matrix forecast tool (Casas et al. 2024), which was shown to agree very well with previous Fisher Matrix forecasts given by EC20. We note that in Casas et al. (2023), no prior was imposed on n_s and Ω_bh², thus no direct comparison is possible. Instead, we used MontePython as a Fisher Matrix forecast tool and compared the result with the MCMC analysis to validate the Fisher Matrix forecast. To be consistent with Casas et al. (2023), we used the fitting formula as the nonlinear model. The 1σ error is very consistent: the Fisher Matrix forecast gives 0.111, while the MCMC analysis gives 0.116. The constraint from MCMC is non-Gaussian and the 1D posterior has a slightly longer tail for large log₁₀|f_R0|. The constraints on cosmological parameters agree very well between the Fisher Matrix forecast and the MCMC result.

4.2.2. Assessing biases for WL

Next, we created the synthetic data using FORGE and we fitted it by different nonlinear models to assess the biases in the recovered parameters due to the difference in the nonlinear modelling. We selected FORGE to create the data because it has a relatively narrow prior range for |f_R0| and we encountered a problem with that range when using FORGE as the model to fit when including baryonic effects. We should note that FORGE has a larger discrepancy with other nonlinear models. As we discussed before, this could be attributed to emulation errors and calibration with halofit. In particular, this choice is disadvantageous to ReACT as the difference of the prediction for Ξ from FORGE is the largest. Thus, the estimation of bias in recovered parameters presented in this section is conservative, particularly for ReACT.

Figure 4 shows the 2D contours of the constraints on parameters, and Table 8 and Table 9 summarise constraints on log₁₀|f_R0| for the optimistic and pessimistic settings. We first start from the optimistic setting. Since ReACT is valid only up to k = 3 h Mpc⁻¹, we did not include ReACT in this case. The fitting formula and e-Mantis recover the input parameters within 1σ. In the case of e-Mantis, cosmological parameters are well recovered, but there is a slight bias in the recovered log₁₀|f_R0|. On the other hand, a slight bias appears in h in the case of the fitting formula. This may be attributed to the fact that the fitting formula does not include cosmological parameter dependence in the prediction of Ξ. To quantify the bias, we define a 1D bias as

Fig. 4. Bias due to different nonlinear modelling. The synthetic data are created by FORGE and fitted by four different models. Left panel: WL optimistic case. Right panel: WL pessimistic case.

$$\begin{array}{c}\hfill B_{1\mathrm{D}}=\frac{\mu -{\mu }_{\mathtt{FORGE}}}{{\sigma }_{\mathtt{FORGE}}},\end{array}$$(29)

where μ and σ are the mean and 1σ error computed from the 1D marginalised posteriors. If μ > μ_FORGE (μ < μ_FORGE), we use the upper (lower) 68.3% confidence interval to obtain σ_FORGE. The 1D bias for log₁₀|f_R0| is B_1D = 0.273 and 0.602 for the fitting formula and e-Mantis, respectively.

We see a similar result in the pessimistic setting for e-Mantis and the fitting formula. In this case, the input log₁₀|f_R0| is well within 1σ although h is again slightly biased. The 1D bias for log₁₀|f_R0| is B_1D = 0.441 and 0.518 for the fitting formula and e-Mantis, respectively, while the 1D bias for h is 0.988 and 0.812 for the fitting formula and e-Mantis, respectively. On the other hand, we find a bias in the recovered log₁₀|f_R0| when ReACT is used. As mentioned above, we should bear in mind that the choice of FORGE as a fiducial is disadvantageous for ReACT. In the case of ReACT, the difference from FORGE leads to biases also in the cosmological parameters leading to larger Ω_ch² and smaller ln(10¹⁰A_s) and h. This is due to the different k dependence of Ξ between ReACT and FORGE and this is compensated by adjusting cosmological parameters as well as log₁₀|f_R0| leading to a stronger bias. The 1D bias for log₁₀|f_R0| reaches B_1D = 3.11. We will discuss how to mitigate this bias by correcting the prediction for the fiducial cosmology and including theoretical errors in Sect. 6.

Since the FORGE prediction deviates from the other three models, we also performed the same analysis using e-Mantis as fiducial data for the pessimistic setting. We show the results in Appendix A. Qualitatively, we obtain consistent results. As we can see from Fig. A.1, the 2D contours overlap in the same way as the case where FORGE is used as fiducial, while the mean of log₁₀|f_R0| obtained by ReACT is closer to the input value. However, we still find that the 1D bias is at 3 sigma level due to the smaller errors of e-Mantis compared with FORGE.

4.2.3. 3 × 2 pt analysis

We considered the constraints from 3 × 2 pt statistics by adding GC_ph and its cross-correlation XC_ph to WL. In this paper, we use a simple model of the scale-independent linear bias as in ΛCDM. This assumption needs to be reexamined in f(R) gravity as the scale-dependent growth will lead to a scale-dependent bias. The effect of f(R) gravity on the halo bias was studied in simulations (Arnold et al. 2019b) and the anlytic model was developed by Valogiannis & Bean (2019). However, the linear bias assumption also needs to be relaxed even in ΛCDM and it is beyond the scope of this paper to implement a more complete bias description. For this reason, we consider the pessimistic case only and focus our attention on a comparison with the Fisher Matrix forecast and check if the bias when FORGE is used to create the data becomes worse with the increasing statistical power. We also do not include ReACT as we already observe a significant bias with WL only.

We have ten scale-invariant bias parameters, one for each redshift bin. For the 3 × 2 pt analysis, we vary n_s, but we still impose a tight Gaussian prior on Ω_b h² as it is unlikely to get a tighter constraint than the Big Bang nucleosynthesis constraints.

Fig. 5. Constraints from 3 × 2 pt analysis in the pessimistic setting. In the left panel, the same nonlinear model is used for the data and the fitting while in the right panel, the data are generated by FORGE.

The parameters used in addition to the WL analysis are listed in Table 10. We did not change the fiducial bias parameters from EC20 as there is no prediction of the linear bias in f(R) gravity. The fiducial bias adopted in EC20 will need to be improved even in ΛCDM. Our prime focus is the study of the effect of nonlinear models and we vary the linear bias in the MCMC analysis resulting in different constraints depending on the nonlinear model used. The observed covariance was built from synthetic data. Thus the effect of f(R) gravity was included in the covariance. We independently checked that the constraints on parameters do not change by using the ΛCDM covariance. Thus we expect that the effect of f(R) on galaxy bias has also negligible effects on the covariance.

The left-hand side of Fig. 5 shows the 2D contours of the constraints on parameters, and Table 11 summarises constraints on log₁₀|f_R0| where the same nonlinear model is used for the data and the model. The right-hand side of Fig. 5 shows the 2D contours of the constraints on parameters, and Table 12 summarises constraints on log₁₀|f_R0| as well as 1D bias for log₁₀|f_R0| when FORGE is used to create the data.

In the case where the same nonlinear model is used for the data and the model, errors are consistent between the fitting formula and e-Mantis, although we still observe the same longer tail for larger log₁₀|f_R0| for the fitting formula as we observe for WL. FORGE gives slightly tighter constraints on log₁₀|f_R0|. Constraints on cosmological parameters are very consistent among the three different nonlinear models. The Fisher Matrix forecast using the fitting formula is also very consistent with the MCMC results. We note that errors from 3 × 2 pt analysis in the pessimistic setting are comparable to or better than WL alone in the optimistic setting. This demonstrates the strength of the 3 × 2 pt analysis although we should bear in mind the limitation of the bias model used in this analysis.

We find that the increased statistical power does not degrade the agreement between FORGE and e-Mantis. When FORGE is used to create the data, the input parameters are well recovered by e-Mantis, although h is slightly biased as in the WL-only case. The 1D bias for log₁₀|f_R0| is given by 0.150. In the case of the fitting formula, the bias in cosmological parameters becomes worse compared with WL. This could be attributed to the fact that the fitting formula does not take into account the cosmological parameter dependence of Ξ as mentioned before. The 1D bias for log₁₀|f_R0| also becomes slightly larger, with a value of B_1D = 0.667.

Fig. 6. Effect of log10|fR0|, two baryonic parameters log10Mc and θej, and h on the WL angular power spectrum in bin 10. The dotted lines show predictions of e-Mantis with varying these parameters individually. The plot shows the ratio to the ΛCDM power spectrum. The plot includes the prediction of FORGE and e-Mantis with the fiducial parameters and the e-Mantis prediction using the best-fit values (log10|fR0| = −5.765, 100h = 65.736, log10(Mc) = 12.999) with the data created by FORGE. We can see that the e-Mantis prediction can fit the FORGE by adjusting these parameters. We note that the cosmological parameters for the ΛCDM power spectrum, $C_{\ell }^{\mathrm{\Lambda }\text{CDM}}$ is fixed to be the fiducial ones when we vary 100h.

5. Baryonic effects

5.1. Adding baryonic effects using BCemu

In this section, we study the impact of baryonic effects on the constraints on the f(R) parameter. We use the seven-parameter emulator of baryonic effects called BCemu (Giri & Schneider 2021) and we assume that the baryonic effects and the modified gravity effects can be treated independently as discussed in Sect. 3 (Arnold & Li 2019). The BCemu parameters govern the gas profiles and stellar abundances in haloes. It is not our intention to study in detail the baryonic effects. We use the default values for the baryonic parameters provided with BCemu and use the full prior range. We choose not to include any redshift dependence. This is because we need to impose a tight prior range for these parameters at z = 0 to keep them within the prior range and we may miss important degeneracies between baryonic parameters and log₁₀|f_R0|. We found that two baryonic parameters, $M_{\mathrm{c}}^{\mathrm{Giri}\mathrm{et}\mathrm{al}.}$ and θ_ej, have the strongest effects on the matter power spectrum and they are well constrained in the presence of the f(R) parameter. Thus, we only vary these two parameters. In BCemu, the gas profile is modelled as a cored double power law. $M_{\mathrm{c}}^{\mathrm{Giri}\mathrm{et}\mathrm{al}.}$ controls the dark matter halo mass dependence of the logarithmic slope of the first-cored power law. It allows the profile to become less steep than the Navarro et al. (1996) one for $M<M_{\mathrm{c}}^{\mathrm{Giri}\mathrm{et}\mathrm{al}.}$. The parameter θ_ej determines the scale radius (with respect to the virial radius) of the second-cored power law. Table 13 summarises the fiducial values and priors for baryonic parameters. We will use the dimensionless parameter M_c that is related to the one defined in Giri & Schneider (2021) via $M_{\mathrm{c}}:=M_{\mathrm{c}}^{\mathrm{Giri}\mathrm{et}\mathrm{al}.}/M_{\odot }$ where M_⊙ is the solar mass.

5.2. Effects of adding baryons

Figure 6 shows the effect of changing log₁₀|f_R0|, h, and two baryonic parameters on the ratio of the WL angular power spectrum to the ΛCDM one. Baryonic effects introduce scale-dependent modifications to the angular power spectrum that are similar to the effect of f(R) at ℓ > 100. On the other hand, h changes the overall amplitude of the angular power spectrum as it changes Ω_m. We will see that the interplay between these parameters leads to interesting degeneracies.

Figure 7 shows the 2D contours of the WL constraints on parameters for the optimistic settings with baryons and Table 14 summarises constraints on log₁₀|f_R0|. In the case of the optimistic setting, when e-Mantis is used as the data and the model, the two baryonic parameters degrade the constraints on log₁₀|f_R0|: the 1σ error becomes 0.281 while it was 0.088 without baryons. However, it is important to stress that log₁₀|f_R0| and the two baryonic parameters are still constrained well within the prior range of these parameters. It means that the impact of modified gravity on the high-k tail of the matter power spectrum is not washed out by the baryonic feedback, and we can still distinguish between the effect of f(R) and baryons.

Fig. 7. Constraints on parameters from the WL optimistic setting with baryons. We show the case where e-Mantis is used to create both the data and fitting and the case where FORGE is used to create the data and it is fitted by e-Mantis. We also include the case where we impose a Gaussian prior on log10 Mc with the width of 0.2 in the case of FORGE fiducial data. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.

Nevertheless, if FORGE is used as the data, the difference between FORGE and e-Mantis in the matter power spectrum can be absorbed by the shift in log₁₀M_c and h, and the constraint on log₁₀|f_R0| is shifted along the degeneracy direction between log₁₀|f_R0| and log₁₀M_c as well as h. This can be understood from Fig. 6. The difference of the scale dependence between e-Mantis and FORGE at ℓ > 100 can be adjusted by decreasing log₁₀|f_R0| and decreasing log₁₀M_c. This leads to a lower amplitude that can be adjusted by decreasing h. Due to the combination of these effects, the best-fit log₁₀|f_R0| becomes smaller. We note that e-Mantis does not include the h dependence in the prediction for Ξ(k, z) explicitly, but it depends on h implicitly through σ₈^ΛCDM. However, this dependence of h on Ξ(k, z) is very weak. To compute the power spectrum in f(R) we use P_f(R) = Ξ(k, z)P_ΛCDM(k, z). The e-Mantis emulator provides us with the first factor, Ξ(k, z), and a ΛCDM emulator provides the last factor. The h dependence comes from the ΛCDM power spectrum not from Ξ(k, z).

As a result, we obtain a 95.5% upper limit of log₁₀|f_R0| as log₁₀|f_R0| < −5.477, which is incompatible with the input value of log₁₀|f_R0| = −5.301. This is partly due to prior volume effects shifting the contour to lower values of log₁₀|f_R0|, caused by the strong degeneracy between log₁₀|f_R0| and log₁₀M_c, as well as h.

To confirm this, we obtain a profile likelihood for log₁₀|f_R0| by fixing log₁₀|f_R0| and finding a minimum χ² by varying other parameters. The profile likelihood was obtained by PROSPECT⁹ (Holm et al. 2024). This is shown in Fig. 8. The global best-fit value for log₁₀|f_R0| is obtained as log₁₀|f_R0| = −5.765, which is larger than the mean of the 1D marginalised posterior log₁₀|f_R0| = −6.087. Figure 6 also shows the prediction of e-Mantis with the best-fit values, which agrees well with FORGE. The Δχ² for the input value of log₁₀|f_R0| = −5.301 is found as Δχ² = 2.657, thus it is still within 2σ. We also observe that the Δχ² curve becomes flat for smaller values of log₁₀|f_R0|.

Fig. 8. Profile likelihood for log10|fR0| where FORGE is used to create data fitted by e-Mantis with baryons in the optimistic setting of WL. The square points show Δχ2 from the best-fit indicated by a diamond. The 95.5% confidence interval from the 1D marginalised posterior is shown as the error bar. The vertical dotted line indicates the input value of −log10|fR0| = 5.301. We see that the 1D marginalised posterior is shifted to smaller values of log10|fR0|.

We also tested the fitting formula as a model to fit to the data generated by FORGE. In this case, the posterior is highly non-Gaussian and the chains did not converge well.

This implies that we can break the degeneracy between f(R) gravity and baryonic effects, but we need to have an accurate nonlinear model for the f(R) gravity model. Otherwise, we could obtain a significantly biased result in terms of the 1D marginalised constraint. We note that if we have a physical prior on log₁₀M_c from baryonic physics, we could break the degeneracy between log₁₀|f_R0| and log₁₀M_c. For example, in the case where FORGE was used as the data and they were modelled by e-Mantis, the mean of log₁₀M_c is significantly shifted to a smaller value as we can see in Fig. 7. If we had a prior on log₁₀M_c to prevent this, this bias could be avoided. To test this idea, we also ran an analysis imposing a prior on log₁₀M_c and showed the result in Fig. 7. We used a Gaussian prior with the width of 0.2, which was estimated from the WL informed gas and stellar mass fraction measurements of massive haloes by Grandis et al. (2024). This prior is consistent with the error on log₁₀M_c that we obtained by using e-Mantis both for data and model. We observe that the bias is relaxed slightly; however, the degeneracy between log₁₀M_c and log₁₀|f_R0| is quite strong and this prior is not enough to alleviate the bias in the marginalised constraints. An improved prior from external data is needed to fully break this degeneracy.

We assumed that the effect of f(R) gravity and baryonic effect can be treated independently. Recently, the coupling between baryonic feedback and cosmology has been studied showing that the combined effect of baryonic and non-baryonic suppression mechanisms is greater than the sum of its parts for decaying dark matter (Elbers 2025). The effect of this coupling on the degeneracy needs to be studied in the future.

Next, we consider the pessimistic setup. Figure 9 shows the 2D contours of the WL constraints on parameters for the optimistic settings with baryons and Table 15 summarises constraints on log₁₀|f_R0|. In this case, due to the weaker constraining power, the 99.7% confidence level upper bound on log₁₀|f_R0| reaches the prior boundary of log₁₀|f_R0| = −7. Also, the 95.5% confidence level upper bound of θ_ej is bounded by the prior. If the data are generated by FORGE, strong degeneracies appear among log₁₀|f_R0| and two baryonic parameters. The input log₁₀|f_R0| is consistent within 1σ for the fitting formula and e-Mantis although the prior bound log₁₀|f_R0| = −7 is reached at 95.5% confidence level lower bound. On the other hand, for ReACT, the input value is outside the 95.5% confidence level upper bound. As we already observed without baryons, h is biased to a lower value.

Fig. 9. Constraints on parameters from the WL pessimistic setting with baryons. As in Fig. 7, we show the case where e-Mantis is used to create both the data and fitting as well as the case where FORGE is used to create the data and it is fitted by e-Mantis, ReACT and the fitting formula. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.

We also checked the case where fiducial data is generated by e-Mantis. As we can see in Fig. A.2, unlike the case without baryons, the agreements between e-Mantis, fitting and ReACT are improved significantly compared with the case where FORGE is used for fiducial data. Particularly, for ReACT, the mean of log₁₀|f_R0| is now consistent with the input value within 1σ.

5.3. Constraints on |f_R0| with the ΛCDM data

Based on these results, we consider the pessimistic setting and use e-Mantis to obtain an upper limit of log₁₀|f_R0| in the presence of baryonic effects to be conservative. Figure 10 shows the 2D contours of the WL constraints on parameters for the pessimistic settings where ΛCDM is used to create the data. We find that the recovered cosmological parameters are consistent with the input values. The mean of 1D marginalised constraints of baryonic parameters are slightly biased due to the strong degeneracy between them, but it is still consistent within 1σ. We obtain the 95.5% confidence level upper bound on log₁₀|f_R0| as

Fig. 10. Constraints on parameters from the WL pessimistic setting with baryons where the data are created in ΛCDM and it is fitted by e-Mantis. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.

$$\begin{array}{c}\hfill {\log }_{10}|f_{R0}|<-5.21.\end{array}$$(30)

We note that this bound depends on the prior log₁₀|f_R0| = −7. To obtain the prior-independent bound, we follow the approach presented in Gordon & Trotta (2007), Piga et al. (2023) and Kou et al. (2024). The ratio of the marginalised posterior and prior is given by

$$\begin{array}{c}\hfill b(x\mid d,p)=\frac{\mathcal{P}(x\mid d,p)}{p(x)},\end{array}$$(31)

where x is the parameter we are interested in (i.e. log₁₀|f_R0|), d is the data, p is the prior, and 𝒫 is the posterior. The Bayes factor, B(x₁, x₂), which quantifies the support of the models with x = x₁ over the models with x = x₂, is given by

$$\begin{array}{c}\hfill B(x_1,x_2)=\frac{b(x_1\mid d,p)}{b(x_2\mid d,p)}=\frac{\mathcal{L}(d\mid x_1)}{\mathcal{L}(d\mid x_2)},\end{array}$$(32)

where ℒ(d ∣ x) is the marginalised likelihood of the data for parameter x. For our purpose, we chose x = log₁₀|f_R0| and fixed x₁ to be the upper bound of the prior, x₁ = −7. Following Gordon & Trotta (2007), we then used B(x₁, x₂) = 2.5 to find x₂ so that the model with x = x₁ is favoured compared to the model with x = x₂ at 95.5% confidence level. This gives the 95.5% confidence interval of log₁₀|f_R0| that does not depend on the prior. We note that this method applies only to a ℒ(d ∣ x) that is a monotonic function of x. Employing this technique, we obtain

$$\begin{array}{c}\hfill {\log }_{10}|f_{R0}|<-5.58.\end{array}$$(33)

Finally, in Fig. 11, we present the constraints on Ω_m and log₁₀|f_R0| so that we can compare these with those in the literature (Schneider et al. 2020; Harnois-Déraps et al. 2023; Spurio Mancini & Bose 2023). It is not possible to make a direct comparison due to various differences in the setting. Nonetheless, our constraint is comparable to the one presented in Schneider et al. (2020) [log₁₀|f_R0| < −5.3] where a similar analysis was done by combining BCemu and the fitting formula. We cannot also make a comparison with Harnois-Déraps et al. (2023) as baryonic effects are not included in their analysis, but again the constraint is comparable [log₁₀|f_R0| < −5.24]. On the other hand, our constraint is much weaker than the one found in Spurio Mancini & Bose (2023) that used ReACT with BCemu. We refer the readers to Spurio Mancini & Bose (2023) for possible explanations. In our analysis, ReACT tended to prefer a lower value of log₁₀|f_R0| when FORGE is used as data and this might be one of the reasons for this difference. Finally, a recent study by Tsedrik et al. (2024) used similar settings to forecast cosmic shear constraints on model-independent parametrisations of modified gravity theories with scale-independent linear growth. They found that a much better understanding of baryonic feedback is needed in order to detect a screening transition.

Fig. 11. Same as in Fig. 10, but we show the constraint on the derived quantity Ωm.

6. Theoretical error

Any theoretical prediction will in general come with an associated error and should be included in the likelihood if it is not subdominant to the statistical error. In this section, we discuss the implementation of theoretical errors in our pipeline. To be conservative, we estimated the theoretical errors using ReACT and FORGE as these two give the most discrepant results. We applied this method to ReACT for the pessimistic case of WL to check if we could remove the parameter biases that we observe for the data created by FORGE at the expense of enlarging error bars.

6.1. Adding uncorrelated theoretical errors to the likelihood

The implementation presented in this work is based on the work of Audren et al. (2013). We have adjusted the recipe to the full 3 × 2 pt likelihood, although we only show the result for the WL analysis in this paper. Here, we discuss the idea behind the formulation. We calculated the angular power spectrum from the power spectrum using Eq. (1). If there is some uncertainty in the modelling of the nonlinear power spectrum, this propagates to the angular power spectrum. As Eq. (1) is a linear functional of the nonlinear matter power spectrum, ℱ[P_δδ], we can propagate the error on the power spectrum, ΔP_δδ, and find

$$\begin{array}{c}\hfill \mathcal{F}[P_{\delta \delta }+\mathrm{\Delta }{P}_{\delta \delta }]=\mathcal{F}[P_{\delta \delta }]+\mathcal{F}[\mathrm{\Delta }{P}_{\delta \delta }]C_{\mathit{ij}}^{\mathit{XY}}(\ell )+E_{\mathit{ij}}^{\mathit{XY}}(\ell ),\end{array}$$(34)

where we have defined the angular power spectrum error, $E_{\mathit{ij}}^{\mathit{XY}}(\ell )$. Following Audren et al. (2013), we defined the relative error on the power spectrum:

$$\begin{array}{c}\hfill \mathrm{\Delta }{P}_{\delta \delta }(k,z)\alpha (k,z)P_{\delta \delta }(k,z).\end{array}$$(35)

We took the most conservative approach and assumed that the error is uncorrelated between different values of ℓ. To account for this uncorrelated theoretical error, we added a new nuisance parameter, ε_ℓ, for each multipole and define the shifted covariance as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\mathit{XY}}(\ell ){\widehat{C}}_{\mathit{ij}}^{\mathit{XY}}(\ell )+{\epsilon }_{\ell }{L}^{1/2}{E}_{\mathit{ij}}^{\mathit{XY}}(\ell ).\end{array}$$(36)

The free parameters, ε_ℓ, have the following meaning: they quantify how many standard deviations the shifted covariance matrix is from the theoretically predicted one. We treat them as a random Gaussian variable with zero mean and a standard deviation of one. The normalisation factor, L^1/2, where L = ℓ_max − ℓ_min + 1, will become clear soon. The log-likelihood function becomes

$$\begin{array}{c}\hfill {\stackrel{\sim }{\chi }}^2\left({\epsilon }_{\ell }\right)={\displaystyle \sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}}[(2\ell +1)f_{\mathrm{sky}}(\frac{{\stackrel{\sim }{d}}_{\ell }^{\mathrm{mix}}\left({\epsilon }_{\ell }\right)}{{\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\left({\epsilon }_{\ell }\right)}+ln\frac{{\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\left({\epsilon }_{\ell }\right)}{d_{\ell }^{\mathrm{obs}}}-N)+{\epsilon }_{\ell }^2].\end{array}$$(37)

The quantities ${\stackrel{\sim }{d}}^{\mathrm{th}}$ and ${\stackrel{\sim }{d}}^{\mathrm{mix}}$ are constructed in the same way as d^th and d^mix [see Eqs. (15) to (17)] using the shifted covariance matrix ${\stackrel{\sim }{C}}_{\mathit{ij}}^{\mathit{XY}}(\ell )$. To include the theoretical error, we vary ε_ℓ and marginalise over them. To a very good approximation, this is equivalent to minimising the ${\stackrel{\sim }{\chi }}^2$ with respect to ε_ℓ at the level of the likelihood. Thus, we define our new χ² as the minimum of the ${\stackrel{\sim }{\chi }}^2$ with respect to ε_ℓ:

$$\begin{array}{c}\hfill {\chi }^2{\displaystyle \underset{{\epsilon }_{\ell }\in {\mathbb{R}}^L}{min}}{\stackrel{\sim }{\chi }}^2.\end{array}$$(38)

The normalisation factor, L^1/2, can be explained as follows. If we were to measure ${\widehat{C}}^{\mathrm{obs}}(\ell )={\widehat{C}}^{\mathrm{th}}(\ell )+E(\ell )$, the minimisation would find ε_ℓ = L^−1/2. The resulting χ² = ∑_ℓε_ℓ² = ∑_ℓL⁻¹ = 1 would match our expectation that a one-sigma theoretical error for each ℓ results in an increase in χ² by one.

The main ingredient of this formulation is the relative error function α(k, z) defined in Eq. (35). The construction of α(k, z) using FORGE and ReACT will be discussed in Appendix B.

6.2. Numerical implementation

The numerical computation of the theory error covariance $E_{\mathit{ij}}^{\mathit{XY}}(\ell )$ follows the prescription presented by Casas et al. (2024). For the minimisation, we can use the fact that all the free ε_ℓ are independent of each other, and we can do the minimisation for each multipole separately. We use Newton’s method to find the minimum. For this, we have to compute the first and second derivatives of the likelihood with respect to ε_ℓ. For any single multipole, we find

$$\begin{array}{c}\hfill \frac{\mathrm{d}\stackrel{\sim }{\chi }}{\mathrm{d}{\epsilon }_{\ell }}=(2\ell +1)f_{\mathrm{sky}}(\frac{\left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{mix}}\right)\prime +\left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\right)\prime }{{\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}}-\frac{{\stackrel{\sim }{d}}_{\ell }^{\mathrm{mix}}\left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\right)\prime }{{\left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\right)}^2})+2{\epsilon }_{\ell }.\end{array}$$(39)

The derivatives of the determinants are computed using Jacobi’s formula. This gives, for example

$$\begin{array}{c}\hfill \left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{th}}\right)\prime =det\left({\widehat{C}}_{\ell }^{\mathrm{th}}\right)L^{1/2}\mathrm{Tr}\left[{\left({\widehat{C}}_{\ell }^{\mathrm{th}}\right)}^{-1}{E}_{\ell }\right],\end{array}$$(40)

and a similar expression for $\left({\stackrel{\sim }{d}}_{\ell }^{\mathrm{mix}}\right)\prime$. We compute the inverse of the covariance numerically and obtain the second derivatives numerically from the first derivative by doing a double-sided three-point stencil. The minimisation would then need to be done for each multipole. For the pessimistic settings, this would correspond to a minimisation in a 1500-dimensional parameter space. To save time, we only do the minimisation on a logarithmically-spaced grid with 100 discrete values. The other values are obtained from an interpolating function. We checked that the effect of the interpolation does not change the results by more than 1%. The obtained ε_ℓ tend to vary continuously with ℓ, as they try to mimic the effects of changing other theory parameters to the observed power spectrum.

6.3. Application to ReACT

We first applied this method to the case in which the data are created by FORGE and the parameter fitting is done by ReACT. Due to the significant bias in the recovered parameters for both |f_R0| and cosmological parameters, we find that the inclusion of theoretical errors is not enough to mitigate the bias. Therefore, we additionally correct the prediction of ReACT by FORGE for the fiducial cosmology as

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{\mathtt{ReACT}}-\mathtt{FORGE}={\mathrm{\Xi }}_{\mathtt{ReACT}}{\left(\frac{{\mathrm{\Xi }}_{\mathtt{FORGE}}}{{\mathrm{\Xi }}_{\mathtt{ReACT}}}\right)}_{\mathrm{fiducial}}.\end{array}$$(41)

The constraints on cosmological parameters are shown in the left panel of Fig. 12. In this case, as expected, we recover the input parameters in an unbiased way when the data are created by FORGE. We see that theoretical errors affect mainly log₁₀|f_R0| and h. In order to check a non-trivial case, next we create the data with e-Mantis. The result is shown in the right panel of Fig. 12. Without theoretical errors, the means of 1D marginalised constraints for log₁₀|f_R0|, Ω_ch², and ln(10¹⁰A_s) are slightly biased compared with the input values. The inclusion of the theoretical errors largely resolves these biases not only by enlarging the error bars but also by making the means closer to the input values. We note that the inclusion of theoretical errors in this case is important to justify the rescaling of Ξ in Eq. (41). We also note that the theoretical error included here overestimates the errors significantly now that we corrected the prediction of ReACT by FORGE for the fiducial cosmology.

Fig. 12. Effect of adding theoretical errors. Left panel: Data are created by FORGE in the WL pessimistic setting and fitted by ReACT corrected by FORGE for fiducial cosmology, Eq. (41), with and without theoretical errors. Right panel: Data are created by e-Mantis in this case.

Finally, we study the impact of adding baryonic effects in the presence of theoretical errors. We again use e-Mantis as the data. The result is shown in Fig. 13. Due to the enlarged errors, constraints on parameters are not affected significantly by baryonic effects. Also, the means of the 1D marginalised constraint remain consistent with the input parameters. However, the inclusion of the theoretical errors changes the degeneracies between cosmological parameters, baryonic parameters and log₁₀|f_R0|. This leads to a tighter lower bound on log₁₀|f_R0| compared with the case where the same data are fitted by e-Mantis itself. This reinforces our conclusion that the 1D marginalised constraint on log₁₀|f_R0| is sensitive to degeneracies among parameters, and the difference in the theoretical predictions strongly affects the constraint.

Fig. 13. Effect of adding theoretical errors in the presence of baryonic effects. Data are created by e-Mantis and fitted by ReACT corrected by FORGE for the fiducial cosmology with theoretical errors. For a comparison, we also show the result obtained by using e-Mantis to fit the same data.

The result shown here is just an illustration of the inclusion of theoretical errors and their effects on the parameter constraints. Our implementation is flexible and it can be applied to any theoretical error described by the relative error function α(k, z) defined in Eq. (35).

7. Conclusion

In this paper, we have studied the effect of using different nonlinear predictions for the dark matter power spectrum on the parameter constraints in the Hu–Sawikci f(R) gravity model obtained from Euclid primary photometric probes. We implemented four different models in the MontePython pipeline to predict angular power spectra for WL, photometric galaxy clustering (GC_ph) and their cross-correlation (XC_ph). Comparing with the N-body simulation data obtained in Paper 1, we found that e-Mantis agrees very well with ECOSMOG, which was used to run simulations to construct the emulator, while FORGE had larger errors compared with the MG-AREPO simulation. The agreement is better for one of the N-body simulations used to construct FORGE obtained in Paper 2 run by MG-AREPO. This indicates that the difference between e-Mantis and FORGE is larger than the one in the baseline N-body simulations (i.e. ECOSMOG and MG-AREPO) mainly due to the way FORGE was constructed. In the Euclid reference cosmology, FORGE gives a larger Ξ, the ratio between the power spectrum in f(R) and in ΛCDM, on all scales compared with e-Mantis and the fitting formula at the 2% level. ReACT underestimates Ξ compared with FORGE more than e-Mantis and the fitting formula at intermediate k.

We used the fiducial value of |f_R0| = 5 × 10⁻⁶ (log₁₀|f_R0| = −5.301) and ran MCMC in the Euclid fiducial cosmology defined in EC20. For the fitting formula, the Fisher Matrix forecast and MCMC results generally agree well although the MCMC result is non-Gaussian. This is partly caused by a lack of cosmological parameter dependence of Ξ, which affects the degeneracy between log₁₀|f_R0| and cosmological parameters for large log₁₀|f_R0|. e-Mantis gives more Gaussian constraints with smaller errors. When FORGE is used to create the data, the 1D mean of log₁₀|f_R0| is not strongly biased in the case of e-Mantis and the fitting formula and the 1D bias is at most 0.6σ. Even for the 3 × 2 pt analysis including all the probes and their cross-correlations, the 1D bias is 0.15σ for log₁₀|f_R0| in the case of e-Mantis.

The impact of baryonic physics on WL was studied by using a baryonification emulator BCemu. For the optimistic setting, the f(R) parameter and two main baryonic parameters are well constrained despite the degeneracies among these parameters. However, the difference in the nonlinear dark matter prediction can be compensated by the adjustment of baryonic parameters as well as cosmological parameters, and the 1D marginalised constraint on log₁₀|f_R0| is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the ΛCDM synthetic data for WL, we obtained the prior-independent bound of log₁₀|f_R0| < −5.6 using e-Mantis.

ReACT shows a large bias in log₁₀|f_R0| as well as cosmological parameters when FORGE was used to create the data. This is because the prediction of ReACT is furthest away from FORGE. We implemented a method to include uncorrelated theoretical errors proposed in Audren et al. (2013) to address this issue. The method is based on the relative error function for the nonlinear dark matter power spectrum. We estimated this using the difference between FORGE and ReACT. We found that the inclusion of theoretical errors alone was not enough to mitigate the bias. We then corrected the prediction of ReACT with FORGE for the fiducial model. We applied this model to the data created by e-Mantis. We found that theoretical errors, in this case, helped reduce the bias not only by enlarging errors but also by making the means of 1D marginalised constraint closer to the input values. When we added baryonic effects with BCemu, errors were not significantly affected and the input values were still recovered. However, the lower bound on log₁₀|f_R0| is tighter than the one obtained by applying e-Mantis itself. This reinforces our conclusion that the 1D marginalised constraint on log₁₀|f_R0| is sensitive to the degeneracies among parameters.

Based on the result of this paper, we draw the following conclusions:

• It is important to check the agreement of different N-body codes that are used to create theoretical predictions. This means not only the code itself, but also accuracy settings. In Paper 1, we found that the accuracy setting such as the refinement criteria in the adaptive mesh refinement method has a large effect on the power spectrum. With the controlled accuracy setting, it is possible to realise 1% agreements between different N-body codes in terms of Ξ(k, z) in the Hu–Sawicki f(R) gravity model.

• We then need to check the accuracy of emulators. We found that FORGE suffers from larger emulation errors and this leads to a larger difference between FORGE and e-Mantis compared with the difference in their baseline N-body code MG-AREPO and ECOSMOG in some cosmologies. Improving the emulation technique will make the agreement better.

• Including baryonic effects can worsen the bias in the 1D marginalised constraint on log₁₀|f_R0|. This is because the difference in the nonlinear dark matter power spectrum prediction can be compensated by the adjustment of baryonic parameters, and the best-fit values are biased. In addition, the degeneracy between baryonic parameters, cosmological parameters and log₁₀|f_R0| leads to a stronger volume effect. This bias can be avoided if there is a prior on baryonic parameters from external datasets for example.

• To account for the uncertainty in the theoretical prediction for the nonlinear power spectrum, it is safer to include theoretical errors. In this paper, we used a conservative error estimation using FORGE and ReACT. An improvement in emulators will make the theoretical error smaller, leading to better constraints. However, we still need to check whether or not the inclusion of baryons will change this conclusion.

The pipeline developed in this paper can be used to test the readiness of the nonlinear power spectrum prediction for the application to real data from Euclid in other extended cosmologies. For the Data Release 1 of Euclid and future data releases, we plan to improve the FORGE emulator and make the pipeline ready for the data analysis. We also plan to extend the analysis to the models studied in Frusciante (2024).

Finally, we note that our forecasts have ignored observational systematics such as shear and redshift measurement biases. For Stage-IV surveys like Euclid it is known that these need to be very well characterised in order to get accurate cosmological parameter estimates. This is a challenge that needs to be addressed for both ΛCDM and exotic cosmologies. A detailed study of this issue is beyond the scope of our paper, but we refer to Euclid Collaboration: Mellier et al. (2025) for details. In summary, the mean galaxy redshifts within the bins need to be known with an accuracy better than ∼0.002(1 + z) if errors in cosmological parameters are not to be degraded. On the other hand, Spurio Mancini & Bose (2023) included a parameter that shifts the mean of the redshift distribution in each redshift bin and found that these parameters are not strongly degenerate with the f(R) parameter.

Acknowledgments

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org). K. K. is supported by STFC grant ST/W001225/1. B. B. is supported by a UK Research and Innovation Stephen Hawking Fellowship (EP/W005654/2). A. P. is a UKRI Future Leaders Fellow [grant MR/X005399/1]. P. C. is supported by grant RF/ERE/221061. M. C. acknowledges the financial support provided by the Alexander von Humboldt Foundation through the Humboldt Research Fellowship program, as well as support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. F. P. acknowledges partial support from the INFN grant InDark and the Departments of Excellence grant L.232/2016 of the Italian Ministry of University and Research (MUR). FP also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021. GR’s research was supported by an appointment to the NASA Postdoctoral Program administered by Oak Ridge Associated Universities under contract with NASA. GR was supported by JPL, which is run under contract by the California Institute of Technology for NASA (80NM0018D0004). Numerical computations were done on the Sciama High Performance Compute (HPC) cluster which is supported by the ICG, SEPNet, and the University of Portsmouth. This work has made use of the Infinity Cluster hosted by Institut d’Astrophysique de Paris. We acknowledge open libraries support IPython (Perez & Granger 2007), Matplotlib (Hunter 2007), Numpy (Walt et al. 2011), and SciPy (Virtanen et al. 2020). For the purpose of open access, we have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. Supporting research data are available on reasonable request from the corresponding author.

References

Appendix A: e-Mantis fiducial data

In this Appendix, we show the result of the cases where we use e-Mantis as fiducial data for the pessimistic setting of WL with and without baryons. Figure A.1 and Table A.1 show the results without baryons while the Figure A.2 and Table A.2 show the results with baryons. We observe that the agreements of the three models are generally better as expected although without baryons, 1D bias is still at the 3σ level for ReACT due to smaller errors of e-Mantis. On the other hand, in the case with baryons, 1D bias is reduced to less than 1σ.

Fig. A.2. Bias due to different nonlinear modelling. The synthetic data are created by e-Mantis and fitted by two different models in the WL pessimistic case with baryons.

Fig. A.1. Bias due to different nonlinear modelling. The synthetic data are created by e-Mantis and fitted by two different models in the WL pessimistic case without baryons.

Appendix B: Construction of the relative error

Fig. B.1. Relative difference between the power spectra from FORGE and ReACT. We show the relative deviation as the lines and the 68.3% confidence level as the shaded area. The boosts were all computed at the fiducial cosmology.

To construct the relative error function, α(k, z), we can use the fact that Ξ(k, z) is constructed as the ratio to the same ΛCDM nonlinear power spectrum calculated with the halofit ‘Takahashi’ prescription in all models. This means that any deviations of the ΛCDM power spectrum from N-body simulations is also absorbed in Ξ. We can thus write the relative error as

$$\begin{array}{c}\hfill \alpha (k,z)=\frac{\mathrm{\Delta }{P}_{\delta \delta }(k,z)}{P_{\delta \delta }(k,z)}=\frac{\mathrm{\Delta }\mathrm{\Xi }(k,z)}{\mathrm{\Xi }(k,z)}.\end{array}$$(B.1)

To construct ΔP_δδ we follow the basic idea that all predictions for Ξ are equally accurate in calculating the modified gravity power spectrum. The true power spectrum is thus only known up to the spread of the predictions. To be conservative with our forecast, we choose to calculate the error from the difference between the FORGE predictions and the ReACT predictions, which is the largest among these predictions. This is shown in Fig. B.1. The ReACT prescription underpredicts the power spectrum at intermediate scales and overpredicts slightly at smaller scales. This is most likely due to ReACT handling the one- and two-halo power spectra separately. HMcode2020 predictions for the pseudo-cosmology nonlinear power spectrum also contribute to these inaccuracies as we can see by comparing the left and right panels of Figure 4 in Cataneo et al. (2019). To simplify the exact difference between these two codes, we argue that the theoretical error should plateau at the maximum deviation. This is because at nonlinear scales, the modes of different wave numbers are no longer independent of each other. By plateauing, we essentially say that we can no longer be more precise in the computation of the power spectrum on smaller scales.

Finally, we construct α using only the difference of the predictions at a fiducial of |f_R0| = 5 × 10⁻⁶ in order not to bias our results. Our final fit has the form

$$\begin{array}{cc}\hfill \alpha (k,z)& =A(z)\frac{x^2+x}{x^2+x+1},\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill A(z)& =\frac{A_1}{exp\left(\frac{z-A_2}{A_3}\right)+1}+A_4,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill k_{\mathrm{p}}(z)& =B_{1}exp\left[\tanh \left(\frac{z-B_2}{B_3}\right)\right],\hfill \end{array}$$(B.4)

where we use x ≔ k/k_p. We separately fit the amplitude function A(z) to the maximum deviation and the plateau wavenumber k_p(z) to the wavenumber of maximum deviation. The best-fit values can be found in Table B.1 and are shown as the shaded region in Fig. B.1. We see that the zero line is within the 68.3% confidence bounds most of the time. This is in accordance with our construction that the two emulators should differ by the theoretical error.

All Tables

All Figures

Fig. 1. Left panel: Comparisons between N-body simulations, e-Mantis and FORGE for |fR0| = 10−5. N-body data from simulations run by MG-AREPO and ECOSMOG codes are taken from the latest comparison project presented in Paper 1. We also include the prediction of FORGE where the cosmological parameters in the fiducial ΛCDM simulation, Ωm = 0.31315 and h = 0.6737, are used to make the prediction (FORGE Node 0). Right panel: Comparison of four prescriptions with an N-body simulation run by MG-AREPO. This simulation is one of the simulations used for training to construct the FORGE emulator (Node 13). The measurement of the power spectrum was done in Paper 2.

In the text

Fig. 2. Comparisons between all models considered for the HS6 fiducial cosmology for angular power spectra of WL (top left), galaxy clustering (top right) and their cross-correlation (bottom right) as well as for the matter power spectrum (bottom left). We show these for the WL redshift bin 6 and galaxy clustering bin 2, with the matter power spectrum plotted at z = 0.5, the redshift at which the kernels of both observables peak for the chosen bins. For all angular spectra, the error bands taken from the diagonal of the covariance are also shown.

In the text

	Fig. 3. Constraints on parameters where the same model is used to create the synthetic data and perform the fitting such that the input parameters are reproduced. Left panel: WL optimistic case. Right panel: WL pessimistic case.
In the text

	Fig. 4. Bias due to different nonlinear modelling. The synthetic data are created by FORGE and fitted by four different models. Left panel: WL optimistic case. Right panel: WL pessimistic case.
In the text

	Fig. 5. Constraints from 3 × 2 pt analysis in the pessimistic setting. In the left panel, the same nonlinear model is used for the data and the fitting while in the right panel, the data are generated by FORGE.
In the text

Fig. 6. Effect of log10|fR0|, two baryonic parameters log10Mc and θej, and h on the WL angular power spectrum in bin 10. The dotted lines show predictions of e-Mantis with varying these parameters individually. The plot shows the ratio to the ΛCDM power spectrum. The plot includes the prediction of FORGE and e-Mantis with the fiducial parameters and the e-Mantis prediction using the best-fit values (log10|fR0| = −5.765, 100h = 65.736, log10(Mc) = 12.999) with the data created by FORGE. We can see that the e-Mantis prediction can fit the FORGE by adjusting these parameters. We note that the cosmological parameters for the ΛCDM power spectrum, $C_{\ell }^{\mathrm{\Lambda }\text{CDM}}$ is fixed to be the fiducial ones when we vary 100h.

In the text

Fig. 7. Constraints on parameters from the WL optimistic setting with baryons. We show the case where e-Mantis is used to create both the data and fitting and the case where FORGE is used to create the data and it is fitted by e-Mantis. We also include the case where we impose a Gaussian prior on log10 Mc with the width of 0.2 in the case of FORGE fiducial data. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.

In the text

Fig. 8. Profile likelihood for log10|fR0| where FORGE is used to create data fitted by e-Mantis with baryons in the optimistic setting of WL. The square points show Δχ2 from the best-fit indicated by a diamond. The 95.5% confidence interval from the 1D marginalised posterior is shown as the error bar. The vertical dotted line indicates the input value of −log10|fR0| = 5.301. We see that the 1D marginalised posterior is shifted to smaller values of log10|fR0|.

In the text

	Fig. 9. Constraints on parameters from the WL pessimistic setting with baryons. As in Fig. 7, we show the case where e-Mantis is used to create both the data and fitting as well as the case where FORGE is used to create the data and it is fitted by e-Mantis, ReACT and the fitting formula. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.
In the text

	Fig. 10. Constraints on parameters from the WL pessimistic setting with baryons where the data are created in ΛCDM and it is fitted by e-Mantis. Left panel: Constraints on baryonic parameters. Right panel: Constraints on cosmological parameters.
In the text

	Fig. 11. Same as in Fig. 10, but we show the constraint on the derived quantity Ωm.
In the text