We study the constraint on fR gravity that can be obtained by photometric primary probes of the Euclid mission. Our focus is the dependence of the constraint on the theoretical modelling of the nonlinear matter power spectrum. In the HuSawicki fR gravity model we consider four different predictions for the ratio between the power spectrum in fR and that in LambdaCDM a fitting formula the halo model reaction approach ReACT and two emulators based on dark matter only Nbody simulations FORGE and eMantis. These predictions are added to the MontePython implementation to predict the angular power spectra for weak lensing WL photometric galaxy clustering and their crosscorrelation. By running Markov Chain Monte Carlo we compare constraints on parameters and investigate the bias of the recovered fR parameter if the data are created by a different model. For the pessimistic setting of WL one dimensional bias for the fR parameter log10fR0 is found to be 0.5 sigma when FORGE is used to create the synthetic data with log10fR0 5.301 and fitted by eMantis. The impact of baryonic physics on WL is studied by using a baryonification emulator BCemu. For the optimistic setting the fR parameter and two main baryon 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 baryon parameters and the onedimensional marginalised constraint on log10fR0 is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting using the LambdaCDM synthetic data for WL we obtain the priorindependent upper limit of log10fR0lt 5.6. Finally we implement a method to include theoretical errors to avoid the bias.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Euclid preparation. Simulations and nonlinearities beyond LambdaCDM. 4. Constraints on fR models from the photometric primary probes
{'Simulations and nonlinearities beyond Λ CDM. 4. Constraints on f ( R ) models from the photometric primary probes': "Euclid Collaboration: K. Koyama ⋆ 1 , S. Pamuk 2 , S. Casas 2 , B. Bose 3 , P. Carrilho 3 , I. Sáez-Casares 4 , L. Atayde 5 , 6 , M. Cataneo 7 , 8 , B. Fiorini 1 , C. Giocoli 9 , 10 , A. M. C. Le Brun 4 , F. Pace 11 , 12 , 13 , A. Pourtsidou 3 , 14 , Y. Rasera 4 , 15 , Z. Sakr 16 , 17 , 18 , H.-A. Winther 19 , E. Altamura 20 , J. Adamek 21 , M. Baldi 22 , 9 , 23 , M.-A. Breton 24 , 25 , 4 , G. Rácz 26 , F. Vernizzi 27 , A. Amara 28 , S. Andreon 29 , N. Auricchio 9 , C. Baccigalupi 30 , 31 , 32 , 33 , S. Bardelli 9 , F. Bernardeau 27 , 34 , C. Bodendorf 35 , D. Bonino 13 , E. Branchini 36 , 37 , 29 , M. Brescia 38 , 39 , 40 , J. Brinchmann 41 , 42 , A. Caillat 43 , S. Camera 11 , 12 , 13 , V. Capobianco 13 , C. Carbone 44 , J. Carretero 45 , 46 , M. Castellano 47 , G. Castignani 9 , S. Cavuoti 39 , 40 , A. Cimatti 48 , C. Colodro-Conde 49 , G. Congedo 3 , C. J. Conselice 20 , L. Conversi 50 , 51 , Y. Copin 52 , F. Courbin 53 , 54 , 55 , H. M. Courtois 56 , A. Da Silva 6 , 5 , H. Degaudenzi 57 , G. De Lucia 31 , M. Douspis 58 , F. Dubath 57 , C. A. J. Duncan 20 , X. Dupac 51 , S. Dusini 59 , M. Farina 60 , S. Farrens 61 , S. Ferriol 52 , P. Fosalba 62 , 24 , M. Frailis 31 , E. Franceschi 9 , S. Galeotta 31 , B. Gillis 3 , P. Gómez-Alvarez 63 , 51 , A. Grazian 64 , F. Grupp 35 , 65 , L. Guzzo 66 , 29 , M. Hailey 67 , S. V. H. Haugan 19 , W. Holmes 26 , F. Hormuth 68 , A. Hornstrup 69 , 70 , P. Hudelot 34 , S. Ili'c 71 , 17 , K. Jahnke 72 , M. Jhabvala 73 , B. Joachimi 74 , E. Keihänen 75 , S. Kermiche 76 , A. Kiessling 26 , M. Kilbinger 61 , B. Kubik 52 , M. Kunz 77 , H. Kurki-Suonio 78 , 79 , P. B. Lilje 19 , V. Lindholm 78 , 79 , I. Lloro 80 , G. Mainetti 81 , D. Maino 66 , 44 , 82 , E. Maiorano 9 , O. Mansutti 31 , O. Marggraf 8 , K. Markovic 26 , M. Martinelli 47 , 83 , N. Martinet 43 , F. Marulli 84 , 9 , 23 , R. Massey 85 , E. Medinaceli 9 , S. Mei 86 , M. Melchior 87 , Y. Mellier 88 , 34 , M. Meneghetti 9 , 23 , E. Merlin 47 , G. Meylan 53 , M. Moresco 84 , 9 , L. Moscardini 84 , 9 , 23 , E. Munari 31 , 30 , C. Neissner 89 , 46 , S.-M. Niemi 90 , C. Padilla 89 , S. Paltani 57 , F. Pasian 31 , K. Pedersen 91 , W. J. Percival 92 , 93 , 94 , V. Pettorino 90 , S. Pires 61 , G. Polenta 95 , M. Poncet 96 , L. A. Popa 97 , L. Pozzetti 9 , F. Raison 35 , A. Renzi 98 , 59 , J. Rhodes 26 , G. Riccio 39 , E. Romelli 31 , M. Roncarelli 9 , R. Saglia 65 , 35 , J.-C. Salvignol 90 , A. G. Sánchez 35 , D. Sapone 99 , B. Sartoris 65 , 31 , M. Schirmer 72 , T. Schrabback 100 , A. Secroun 76 , G. Seidel 72 , S. Serrano 62 , 101 , 24 , C. Sirignano 98 , 59 , G. Sirri 23 , L. Stanco 59 , J. Steinwagner 35 , P. Tallada-Crespí 45 , 46 , A. N. Taylor 3 , I. Tereno 6 , 102 , R. Toledo-Moreo 103 , F. Torradeflot 46 , 45 , I. Tutusaus 17 , L. Valenziano 9 , 104 , T. Vassallo 65 , 31 , G. Verdoes Kleijn 105 , A. Veropalumbo 29 , 37 , 106 , Y. Wang 107 , J. Weller 65 , 35 , G. Zamorani 9 , E. Zucca 9 , A. Biviano 31 , 30 , E. Bozzo 57 , C. Burigana 108 , 104 , M. Calabrese 109 , 44 , D. Di Ferdinando 23 , J. A. Escartin Vigo 35 , G. Fabbian 110 , 111 , R. Farinelli 9 , F. Finelli 9 , 104 , J. Gracia-Carpio 35 , S. Matthew 3 , N. Mauri 48 , 23 , A. Pezzotta 35 , M. Pöntinen 78 , V. Scottez 88 , 112 , M. Tenti 23 , M. Viel 30 , 31 , 33 , 32 , 113 , M. Wiesmann 19 , Y. Akrami 114 , 115 , S. Anselmi 59 , 98 , 4 , M. Archidiacono 66 , 82 , F. Atrio-Barandela 116 , M. Ballardini 117 , 9 , 118 , D. Bertacca 98 , 64 , 59 , A. Blanchard 17 , L. Blot 119 , 4 , H. Böhringer 35 , 120 , 121 , S. Bruton 122 , R. Cabanac 17 , A. Calabro 47 , B. Camacho Quevedo 62 , 24 , G. Cañas-Herrera 90 , 123 , A. Cappi 9 , 124 , F. Caro 47 , C. S. Carvalho 102 , T. Castro 31 , 32 , 30 , 113 , K. C. Chambers 125 , S. Contarini 35 , A. R. Cooray 126 , G. Desprez 127 , A. Díaz-Sánchez 128 , J. J. Diaz 129 , S. Di Domizio 36 , 37 , H. Dole 58 , S. Esco ffi er 76 , M. Ezziati 43 , A. G. Ferrari 48 , 23 , P. G. Ferreira 130 , I. Ferrero 19 , A. Finoguenov 78 , A. Fontana 47 , F. Fornari 104 , L. Gabarra 130 , K. Ganga 86 , J. García-Bellido 114 , T. Gasparetto 31 , V. Gautard 131 , E. Gaztanaga 24 , 62 , 1 , F. Giacomini 23 , F. Gianotti 9 , G. Gozaliasl 132 , 78 , C. M. Gutierrez 133 , A. Hall 3 , H. Hildebrandt 7 , J. Hjorth 91 , A. Jimenez Muñoz 134 , S. Joudaki 1 , J. J. E. Kajava 135 , 136 , V. Kansal 137 , 138 , D. Karagiannis 139 , 140 , C. C. Kirkpatrick 75 , J. Le Graet 76 , L. Legrand 141 , J. Lesgourgues 2 , T. I. Liaudat 142 , S. J. Liu 60 , A. Loureiro 143 , 144 , G. Maggio 31 , M. Magliocchetti 60 , F. Mannucci 145 , R. Maoli 146 , 47 , J. Martín-Fleitas 147 , C. J. A. P. Martins 148 , 41 , L. Maurin 58 , R. B. Metcalf 84 , 9 , M. Miluzio 51 , 149 , P. Monaco 150 , 31 , 32 , 30 , A. Montoro 24 , 62 , A. Mora 147 , C. Moretti 33 , 113 , 31 , 30 , 32 , G. Morgante 9 , C. Murray 86 , S. Nadathur 1 , Nicholas A. Walton 110 , L. Pagano 117 , 118 , L. Patrizii 23 , V. Popa 97 , D. Potter 21 , P. Reimberg 88 , I. Risso 106 , P.-F. Rocci 58 , M. Sahlén 151 , E. Sarpa 33 , 113 , 32 , A. Schneider 21 , M. Sereno 9 , 23 , A. Silvestri 123 , A. Spurio Mancini 152 , 67 , J. Stadel 21 , K. Tanidis 130 , C. Tao 76 , N. Tessore 74 , G. Testera 37 , R. Teyssier 153 , S. Toft 154 , 155 , S. Tosi 36 , 37 , A. Troja 98 , 59 , M. Tucci 57 , J. Valiviita 78 , 79 , D. Vergani 9 , G. Verza 156 , 157 , and P. Vielzeuf 76  \n(A ffi liations can be found after the references)  \nSeptember 6, 2024", 'ABSTRACT': 'We study the constraint on f ( R ) gravity that can be obtained by photometric primary probes of the Euclid mission. Our focus is the dependence of the constraint on the theoretical modelling of the nonlinear matter power spectrum. In the Hu-Sawicki f ( R ) gravity model, we consider four di ff erent predictions for the ratio between the power spectrum in f ( R ) and that in Λ CDM: a fitting formula, the halo model reaction approach, ReACT and two emulators based on dark matter only N -body simulations, FORGE and e-Mantis . These predictions are added to the MontePython implementation to predict the angular power spectra for weak lensing (WL), photometric galaxy clustering and their cross-correlation. By running Markov Chain Monte Carlo, we compare constraints on parameters and investigate the bias of the recovered f ( R ) parameter if the data are created by a di ff erent model. For the pessimistic setting of WL, one dimensional bias for the f ( R ) parameter, log10 | fR 0 | , is found to be 0 . 5 σ when FORGE is used to create the synthetic data with log10 | fR 0 | = -5 . 301 and fitted by e-Mantis . The impact of baryonic physics on WL is studied by using a baryonification emulator BCemu . For the optimistic setting, the f ( R ) parameter and two main baryon parameters are well constrained despite the degeneracies among these parameters. However, the di ff erence in the nonlinear dark matter prediction can be compensated by the adjustment of baryon parameters, and the one-dimensional marginalised constraint on log10 | fR 0 | is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the Λ CDMsynthetic data for WL, we obtain the prior-independent upper limit of log10 | fR 0 | < -5 . 6. Finally, we implement a method to include theoretical errors to avoid the bias due to inaccuracies in the nonlinear matter power spectrum prediction.  \nKey words. Cosmology: theory; large-scale structure of Universe; cosmological parameters; dark energy. Gravitational lensing: weak', '1. Introduction': "In 1998, astronomers made a 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. Is the acceleration 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 2019). Testing gravity is one of the main objectives of stage-IV dark energy surveys (Albrecht et al. 2006).  \nOf particular importance in these surveys is the Euclid mission (Euclid Collaboration: Mellier et al. 2024). The Euclid satellite undertakes a spectroscopic survey of galaxies and an imaging survey (targeting weak lensing 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).  \nModified 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-GabadadzePorratti (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. 2019; 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  \nenergy and modified gravity models has been developed (Cataneo et al. 2019; Bose et al. 2021; Bose et al. 2023; Carrilho et al. 2022). These nonlinear predictions were used to study their modifications of the Weak Lensing (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. 2023) for example in f ( R ) gravity.  \nIn 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 | fR 0 | = 5 × 10 -6 , Euclid alone will be able to constrain the additional parameter log10 | fR 0 | 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 linearly scale-independent growth (Frusciante et al. 2023).  \nFor real data analysis, it is imperative to check the e ff ect 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:  \n- 1. Numerical methods and validation (Euclid Collaboration: Adamek et al. 2024, paper 1 hereafter)\n- 2. Results from non-standard simulations (Euclid Collaboration: Racz et al. 2024, paper 2 hereafter)\n- 3. Cosmological constraints on non-standard cosmologies from simulated Euclid probes (D'Amico et al. in prep.)\n- 4. Constraints on f ( R ) models from the photometric primary probes (this work)  \nIn 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 will perform Markov Chain Monte Carlo (MCMC) simulations using synthetic data for Euclid photometric probes, and assess the impact of using di ff erent nonlinear models for the  \nmatter power spectrum at the level of parameter constraints. In addition, we will add baryonic e ff ects using a baryonification method, which 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 Bose et al. (2023) and D'Amico et al. (paper 3) will perform a similar analysis to this paper's for the spectroscopic probes.  \nThis 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 will be presented based on the synthetic data created by four di ff erent 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 di ff erent nonlinear model. Section 5 is devoted to the study of baryonic e ff ects using the BCemu baryonification method. We show how the bias is a ff ected by baryonic e ff ects and obtain the upper bound on | fR 0 | 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 (GCph) and their crosscorrelation (XCph). In EC20, these were calculated using the Limber and flat-sky approximations in a flat Λ CDM Universe, as  \nC XY i j ( ℓ ) = c Z z max z min d z W X i ( z ) W Y j ( z ) H ( z ) r 2 ( z ) P δδ ( k ℓ, z ) , (1)  \nwhere W X i ( z ) is the window function in each tomographic redshift bin i with X = { L , G } corresponding to WL and GCph, respectively, k ℓ = ( ℓ + 1 / 2) / r ( z ), r ( z ) is the comoving distance to redshift z , P δδ ( k ℓ, z ) is the nonlinear power spectrum of matter density fluctuations, δ , at wave number k ℓ and redshift z , in the redshift range of the integral from z min = 0 . 001 to z max = 4 and H ( z ) is the Hubble function.  \nHowever, 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)  \nP ϕ W ( k , z ) = " -3 Ω m GLYPH<18> H 0 c GLYPH<19> 2 (1 + z ) Σ ( k , z ) # 2 P δδ ( k , z ) ,  \n(2) N ℓ = " N LL i j ( ℓ ) N GL i j ( ℓ ) N LG i j ( ℓ ) N GG i j ( ℓ ) # , (9)  \nwhere Σ ( 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 0 is the Hubble parameter - both at the present time. Here we assumed a standard background evolution of the matter component, i.e. ρ m( z ) = ρ m , 0(1 + z ) 3 .  \nWe 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)  \nW G i ( k , z ) = 1 bi ( k , z ) ni ( z ) H ( z ) ,  \nc (3) W L i ( k , z ) = 3 2 Ω m GLYPH<18> H 0 c GLYPH<19> 2 (1 + z ) r ( z ) Σ ( k , z ) × Z z max z d z \' ni ( z ) r ( z \' ) -r ( z ) r ( z \' ) + W IA i ( k , z ) , (4)  \nwhere ni ( z ) is the normalised galaxy number-density distribution in a tomographic redshift bin i such that R z max z min ni ( z ) dz = 1, and W IA i ( k , z ) encodes the contribution of intrinsic alignments (IA) to the WL power spectrum. We follow EC20 in assuming an e ff ective scale-independent galaxy bias, constant within each redshift bin, and its values bi are introduced as nuisance parameters in our analysis, with their fiducial values determined by bi = √ 1 + ¯ zi , where ¯ zi is the mean redshift of each redshift bin.  \nThe IA contribution is computed following the nonlinear alignment model with a redshift dependent amplitude (EC20), in which  \nW IA i ( k , z ) = -A IA C IA Ω m F IA( z ) δ ( k , z ) / δ ( k , 0) ni ( z ) H ( z ) c , (5)  \nwhere  \nF IA( z ) = (1 + z ) η IA . (6)  \nThe parameters A IA and η IA are the nuisance parameters of the model, and C IA is a constant accounting for dimensional units. The galaxy distribution is binned into 10 equipopulated redshift bins with an overall distribution following  \nn ( z ) ∝ z z 0 ! 2 exp        -z z 0 ! 3 / 2        , (7)  \n√  \nChanges 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 GCph case are accounted for in the new P δδ ( k ℓ, z ).  \nwith z 0 = 0 . 9 / 2 and the normalisation set by the requirement that the surface density of galaxies is ¯ n g = 30 arcmin -2 (EC20).  \nFinally, we present the likelihood that we will use. We follow the approach presented in Casas et al. (2024). We first construct an ( N G + N L ) × ( N G + N L ) angular power spectrum matrix for each multipole, where the di ff erent N correspond to the number of redshift bins for each probe (WL and GCph):  \nC ℓ = " C LL i j ( ℓ ) C GL i j ( ℓ ) C LG i j ( ℓ ) C GG i j ( ℓ ) # , (8)  \nwhere lower-case Latin indexes i , . . . run over all tomographic bins. Similarly, the noise contributions can also be condensed into one noise matrix N :  \nwhere the noise terms N AB i j ( ℓ ) are given by  \nK  \nN LL i j ( ℓ ) = δ i j ¯ ni σ 2 ϵ , (10) N GG i j ( ℓ ) = δ K i j ¯ ni , (11) N GL ( ℓ ) = 0 , (12)  \ni j  \nwhere σ 2 ϵ = 0 . 3 2 is the variance of observed ellipticities. We can further define ˆ C ℓ B C ℓ + N ℓ . This is the covariance of the spherical multipole moments alm .  \nIn 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  \nˆ C obs i j ( ℓ ) = 1 2 ℓ + 1 ℓ X m = -ℓ ( a ℓ m ) i ( a ℓ m ) ∗ j , (13)  \nand the theoretically predicted one ˆ C th i j ( ℓ ) as  \nχ 2 = f sky ℓ max X ℓ = ℓ min (2 ℓ + 1) " d mix d th + ln d th d obs ! -N # , (14)  \nwhere the determinants d are defined as  \nd th ( ℓ ) = det h ˆ C th i j ( ℓ ) i , (15) d obs ( ℓ ) = det h ˆ C obs i j ( ℓ ) i , (16)  \nd mix ( ℓ ) = N X k det            ˆ C obs i j ( ℓ ) for k = j ˆ C th i j ( ℓ ) for k , j       . (17)  \nHere N is the number of bins and thus either ( N G + N L ) for multipoles where we include the cross correlation, or the respective N for multipoles where we treat the two probes separately. The additional factor f sky comes from an approximation to account for having less available independent ℓ modes due to partial sky coverage. In this paper, we set the observed covariance ˆ C 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. (2023).  \nIn 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 GCph are k max = [0 . 7 , 2 . 9] h Mpc -1 , respectively, while for the WL, maximum wavenumbers probed are k max = [1 . 4 , 4 . 8] h Mpc -1 at the peak redshift z = 0 . 25 (Casas et al. 2023). Here h denotes the dimensionless Hubble parameter h B H 0 / (100kms -1 Mpc -1 ). For smaller values of z , the values of k at a given ℓ increase, but the window functions in Eqs. (3) and (4)  \nWe consider two di ff erent scenarios: an optimistic and a pessimistic case. In the optimistic case, we consider ℓ max = 5000 for WL, and ℓ max = 3000 for GCph and XCph. Instead, in the pessimistic scenario, we consider ℓ max = 1500 for WL, and ℓ max = 750 for GCph and XCph.  \nTable 1. Euclid survey specifications for WL, GCph and GCsp used by EC20 and Casas et al. (2024).  \nsuppress the power spectrum and we set it to zero after a fixed k max = 30 h 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 di ff ers 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 ), i.e.  \nS = c 4 16 π G Z d 4 x √ -g GLYPH<2> R + f ( R ) GLYPH<3> + S m[ g µν, Ψ m] , (18)  \nwhere 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).  \nHow this modification changes gravity is more easily seen by formulating the theory in the so-called Einstein frame by performing a conformal transformation g µν = ˜ g µν A 2 ( ϕ ) where A ( ϕ ) = p 1 + fR = e ϕ/ √ 6 to obtain  \nS = c 4 16 π G Z d 4 x p -˜ g " ˜ R + 1 2 ( ∂ϕ ) 2 -V ( ϕ ) # (19) + S m[ A 2 ( ϕ ) ˜ g µν, Ψ m] ,  \nwith the potential V = [ fRR -f ( R )] / 2(1 + fR ) 2 and fR B d f ( R ) / d R . 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:  \nF fifth = 1 3 GMm r 2 (1 + r /λ C) e -r /λ C . (20)  \nThis e ff ectively changes G → 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 e ff ectively suppresses the fifth-force by a factor ∝ fR / Φ N where Φ N is the Newtonian gravitational potential.  \nNot all f ( R ) models one can write down possess such a screening mechanism as this 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 largecurvature limit is defined by  \nf ( R ) = -6 Ω DE H 2 0 c 2 + | fR 0 | ¯ R 2 0 R , (21)  \nwhere  \n¯ R 0 = 3 Ω m H 2 0 c 2 1 + 4 Ω DE Ω m ! , (22)  \nis 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 | fR 0 | . In the limit | fR 0 | → 0, we recover GR and the Λ CDM model. This parameter controls the range of the fifth-force and, in the cosmological background, we have λ C ≃ 32 p | fR 0 | / 10 -4 Mpc at the present time. Solar System constraints require | fR 0 | ≲ 10 -6 , cosmological constraints currently lie around 10 -6 -10 -4 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 -8 (Desmond & Ferreira 2020), but see Burrage et al. (2024) for a recent note of caution on such galactic-scale constraints.  \nThe 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 | fR 0 | Ω 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 a ff ected as follows  \nΣ ( z ) = 1 1 + fR ( z ) . (23)  \nSince the maximum value of | fR ( z ) | is given by | fR 0 | , for the values of | fR 0 | we consider in this paper, we can ignore this e ff ect. 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 e ff ect 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  \nΞ ( k , z ) B Pf ( R )( k , z ) P Λ CDM( k , z ) , (24)  \nto 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 e ff ect 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. 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 di ff erent prescriptions for Ξ ( k , z ). Given this, the Λ CDM nonlinear spectrum prescription does not matter and it is common to all nonlinear models.", '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 | fR 0 | . 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 ELEPHANT N-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 di ff erent σ 8, Ω 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.  \nThe fitting function has in total 54 parameters for the full scale, redshift and | fR 0 | dependence. This fitting formula is not defined outside the range 10 -7 < | fR 0 | < 10 -4 . The code is publicly available ( /github ).', '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 nonstandard physics. The nonlinear power spectrum is given by  \nP NL( k , z ) = R ( k , z ) P pseudo NL ( k , z ) , (25)  \nwhere P pseudo NL ( 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.  \nThe halo model reaction, R ( 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 et al. (2023) for the exact expressions for this term. The halo model reaction can be computed e ffi ciently  \nusing the publicly available ReACT (Bose et al. 2020, 2021; Bose et al. 2023, /github ) code.  \nDespite ReACT being highly e ffi cient, 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 create a neural network emulator using the Cosmopower package (Spurio Mancini et al. 2022, /github ) for the halo model reaction-based boost  \nΞ ( k , z ) = R ( k , z ) P pseudo HMCode2020 ( k , z ) P Λ CDM( k , z ) , (26)  \nwhere P Λ CDM( k , z ) and P pseudo HMCode2020 ( k , z ) are calculated using HMCode2020 (Mead et al. 2021, /github ) while R ( k , z ) is calculated using ReACT . We chose the HMCode2020 to model the pseudopower 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.  \nTo train the emulator, we follow the procedure of Spurio Mancini & Bose (2023) but widen the parameter priors. We produce ∼ 10 5 boost predictions in the range k ∈ [0 . 01 , 3] h Mpc -1 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 ( /github ). Finally, we note another small di ff erence 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 Ω m = Ω b + Ω c, Ω b and Ω c being the baryon and cold dark matter density parameters today respectively. 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 | fR 0 | , Ω m, σ Λ CDM 8 and h with all other parameters fixed. We note that σ Λ CDM 8 is the σ 8 we would obtain in a Λ CDM model with the same cosmological parameters and initial amplitude A s and not σ 8 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 -1 . f ( R ) simulations are performed by a modified version of the Arepo code, MG-AREPO (Springel 2010; Arnold et al. 2019; Weinberger et al. 2020) that solves the nonlinear scalar field equation using a relaxation method. See paper 1 for more details.  \nThe emulation was made for the ratio between the power spectrum in f ( R ) and the halofit prediction in Λ CDM. We note that this is di ff erent 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 e ff ect, 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 8 = 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.  \nThe latest version of the FORGE simulations has Λ CDMcounterparts 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 ( /f171 ) for this paper. We will 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 n | fR 0 | , Ω m , σ Λ CDM 8 o 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 -1 and cosmological parameters in the following range: GLYPH<12> GLYPH<12> GLYPH<12> fR 0 GLYPH<12> GLYPH<12> GLYPH<12> ∈ h 10 -7 , 10 -4 i , Ω m ∈ [0 . 2365 , 0 . 3941] and σ Λ CDM 8 ∈ [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 -1 and across the whole parameter space although in most cases, the accuracy is of order 1%. The code is publicly available ( /\\_621 ).', '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 di ff erent implementations have di ff erent 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 a ff ect the final results strongly. We checked that the contributions from high redshifts and high wavenumbers which used extrapolations were sub-dominant.  \nFor the modified gravity parameter | fR 0 | , we use flat priors in terms of log10 | fR 0 | 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 e ff ect of the fifth-force has to vanish on these scales. For small scales, we do a power-law extrapolation. To obtain the spectral index for the extrapolation, we proceed as follows: we calculate Ξ on a grid close to the edges of the emulators; we then fit, for fixed z , a linear function in the log-log space onto this grid. This is done to average out any numerical noise at the edge. We do the same for high redshifts  \nTable 2. A summary of four nonlinear models that we use in this paper. For the details of N-body simulations codes and simulation suites, see paper 1 and paper 2. The prediction for Ξ remains the same if we vary cosmological parameters which are not included in the prediction.Table 3. Ranges of wavenumbers, redshifts and log10 | fR 0 | in di ff erent emulators used in this work.  \nby fitting the power law for fixed wavenumbers. For regions of both high k and z , we do a constant extrapolation of the spectral index. 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 di ff erent emulators are found in Table 3. We do 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. The validity range of Λ CDM cosmological parameters is summarised in Table 4.  \nWe also add the e ff ect of baryonic feedback in the form of an additional correction Ξ BFM to our power spectrum prediction. The physics of the baryonic feedback e ff ects is discussed in Sect. 5. We refer the reader to Schneider et al. (2020) and Mead et al. (2021) for further details. We obtain 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 -1 . In this case, we use constant extrapolation for both k and z . We checked that this had little e ff ect on our results. This is because the Euclid main probes are mostly sensitive to redshifts around z ∼ 1. For these redshifts, the extrapolation only a ff ects very high multipoles ℓ ≳ 4500. Thus, the extrapolation has negligible e ff ect on the angular power spectrum. We estimate the overall e ff ect 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 e ff ect of modified gravity, we combine both boosts as  \nP BFM f ( R ) ( k , z ) = Ξ BFM ( k , z ) Ξ ( k , z ) P Λ CDM( k , z ) .  \nWe can make this approximation if both e ff ects 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. (2019) and Arnold & Li (2019).  \n8  \nOur 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 di ff erent 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 simulations to compare errors and investigate bias due to the di ff erence in the prediction of the nonlinear matter power spectrum using the settings defined in EC20.', '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 e ff ect 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 Λ CDMcosmology with Ω m = 0 . 269, h = 0 . 704, n s = 0 . 966 and σ 8 = 0 . 801. The simulations have N p = 512 3 particles of mass M p ≃ 8 . 756 × 10 9 h -1 M ⊙ in a box of size B = 250 h -1 Mpc and start at redshift z = 49. We picked a model called F5 with | fR 0 | = 10 -5 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 8 = 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.  \nTo further investigate this issue with FORGE , we use 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 | fR 0 | . This simulation has the following cosmological  \nTable 4. Ranges for cosmological parameters implemented in the MontePython code.Fig. 1. Left panel : Comparisons between N -body simulations, e-Mantis and FORGE for | fR 0 | = 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.  \n<!-- image -->  \nparameters: Ω m = 0 . 34671, h = 0 . 70056 and σ 8 = 0 . 78191 and we show a comparison at z = 0 . 652. Both Λ CDM and f ( R ) simulations with | fR 0 | = 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 -1 . In this case, both FORGE and e-Mantis agree with MG-AREPO within 1% up to k = 3 h Mpc -1 . The fitting formula and ReACT agree with MG-AREPO within 1% up to k = 1 h Mpc -1 . 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 Λ CDMpower 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  \nΘ = { Ω m , Ω b , h , n s , σ 8 , log10 | fR 0 |} , HS6 : Θ fid = { 0 . 32 , 0 . 05 , 0 . 67 , 0 . 96 , 0 . 853 , -5 . 301 } . (27)  \nThe cosmological parameters are the same ones adopted in EC20. Casas et al. (2023) showed that this value of | fR 0 | = 5 × 10 -6 can be well constrained by the Euclid photometric  \nprobes. Also the range of | fR 0 | 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 P m ν = 0 . 06 eV, but we keep P m ν fixed in the following analysis.  \nFigure 2 shows a comparison of the power spectrum and the angular power spectrum for WL, GCph and their crosscorrelation XCph. 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 -1 . 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 Λ CDM FORGE simulation ( Ω m = 0 . 31315). On the other hand, ReACT underestimates Ξ at k < 0 . 1 h Mpc -1 and even predicts Ξ < 1. As we discussed in the previous section, we enforce Ξ ≥ 1.  \nFig. 2. Comparisons between all models considered for the HS6 fiducial cosmology for angular power spectra of weak lensing ( 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 weak lensing 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.  \n<!-- image -->', '4.2.1. Forecasting errors for WL': 'We first compare 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 | fR 0 | parameter and cosmological parameters. We first consider the WL-only case. In this case, we impose a tight Gaussian prior on the spectral index, n s, taken from the Planck results (Planck Collaboration: Aghanim et al. 2020) and on the baryon density parameter, using Big Bang nucleosynthesis constraints (Pisanti et al. 2021):  \nn s = 0 . 96 ± 0 . 004 , Ω b h 2 = 0 . 022445 ± 0 . 00036 , (28)  \nas 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 use a GelmanRubin (Gelman & Rubin 1992) value of R -1 < 0 . 01 for each individual sampling parameter using MontePython . For postprocessing chains, we use GetDist (Lewis 2019, /github ).  \nFigure 3 shows the 2D contours of the constraints on parameters, and Table 6 and Table 7 summarise constraints on log10 | fR 0 | for the optimistic and pessimistic settings. The constraints on cosmological parameters are consistent among the four di ff erent nonlinear models, while we see some notable differences in the constraints on | fR 0 | . 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  \nTable 5. Fiducial values and flat prior ranges for cosmological parameters, | fR 0 | and IA parameters.  \nsmall log10 | fR 0 | but has a weaker constraint for larger log10 | fR 0 | . Constraints from FORGE agree with the fitting formula for large log10 | fR 0 | , but give a weaker constraint for small log10 | fR 0 | . The degeneracy between log10 | fR 0 | , Ω c h 2 and ln(10 10 A s) also presents some notable di ff erences. The degeneracy for larger log10 | fR 0 | is di ff erent 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 log10 | fR 0 | . ReACT gives a weaker constraint on log10 | fR 0 | , but the constraints on cosmological parameters are consistent among the four di ff erent 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.  \nWe compare 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 Ω b h 2 , thus no direct comparison is possible. Instead, we use MontePython as a Fisher Matrix forecast tool and compare the result with the MCMC analysis to validate the Fisher Matrix forecast. To be consistent with Casas et al. (2023), we use 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 log10 | fR 0 | . 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 create the synthetic data using FORGE and we fit it by di ff erent nonlinear models to assess the biases in the recovered parameters due to the di ff erence in the nonlinear modelling. We selected FORGE to create the data because it has a relatively narrow prior range for | fR 0 | and we encountered a problem with that range when using FORGE as the model to fit when including baryonic e ff ects. We should note that FORGE has a larger discrepancy from 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  \nFigure 4 shows the 2D contours of the constraints on parameters, and Table 8 and Table 9 summarise constraints on log10 | fR 0 | for the optimistic and pessimistic settings.  \nWefirst start from the optimistic setting. Since ReACT is valid only up to k = 3 h Mpc -1 , we will not include ReACT in this case.  \nThe 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 log10 | fR 0 | . 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  \nB 1D = µ -µ FORGE σ FORGE , (29)  \nwhere µ 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 log10 | fR 0 | is B 1D = 0 . 273 and 0 . 602 for the fitting formula and e-Mantis , respectively.  \nWe see a similar result in the pessimistic setting for e-Mantis and the fitting formula. In this case, the input log10 | fR 0 | is well within 1 σ although h is again slightly biased. The 1D bias for log10 | fR 0 | 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 log10 | fR 0 | 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 di ff erence from FORGE leads to biases also in the cosmological parameters leading to larger Ω c h 2 and smaller ln(10 10 A s) and h . This is due to the di ff erent k dependence of Ξ between ReACT and FORGE and this is compensated by adjusting cosmological parameters as well as log10 | fR 0 | leading to a stronger bias. The 1D bias for log10 | fR 0 | 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.  \nSince 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 log10 | fR 0 | 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 × 2pt analysis': 'We consider the constraints from 3 × 2pt statistics by adding GCph and its cross-correlation XCph 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. 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  \nFig. 4. Bias due to di ff erent nonlinear modelling. The synthetic data are created by FORGE and fitted by four di ff erent models. Left panel : WL optimistic case. Right panel : WL pessimistic case.  \n<!-- image -->  \n|  \n|  \n<!-- image -->  \n10  \n|  \n|  \nFig. 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.  \nTable 6. The mean, standard deviation and 68.3% upper and lower limit of log10 | fR 0 | for the WL optimistic setting where the same prediction is used to create the data and the model.  \n<!-- image -->  \n|  \n|  \nTable 7. The mean, standard deviation and 68.3% upper and lower limit of log10 | fR 0 | for the WL pessimistic setting where the same model is used to create the data and the model.  \n.  \n.  \n.  \n<!-- image -->  \n|  \n|  \nTable 8. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the optimistic setting where the data are generated by FORGE and it is fitted by di ff erent nonlinear models.  \n.  \n.  \n.  \ncomparison with the Fisher Matrix forecast and check if the bias when FORGE is used to create the data becomes worse with the  \nincreasing statistical power. We also do not include ReACT as we already observe a significant bias with WL only.  \nTable 9. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the pessimistic setting where the data are generated by FORGE and it is fitted by di ff erent nonlinear models.  \n.  \n.  \n.  \nWehave 10 scale-invariant bias parameters, one for each redshift bin. For the 3 × 2pt analysis, we vary n s, but we still impose a tight Gaussian prior on Ω b h 2 as it is unlikely to get a tighter constraint than the Big Bang nucleosynthesis constraints.  \nThe 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 e ff ect of nonlinear models and we vary the linear bias in the MCMC analysis resulting in di ff erent constraints depending on the nonlinear model used. The observed covariance is built from synthetic data. Thus the e ff ect of f ( R ) gravity is included in the covariance. We independently checked that the constraints on parameters do not change by using the Λ CDMcovariance. Thus we expect that the e ff ect of f ( R ) on galaxy bias has also negligible e ff ects on the covariance.  \nThe left-hand side of Fig. 5 shows the 2D contours of the constraints on parameters, and Table 11 summarises constraints on log10 | fR 0 | 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 log10 | fR 0 | as well as 1D bias for log10 | fR 0 | when FORGE is used to create the data.  \nIn 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 log10 | fR 0 | for the fitting formula as we observe for WL. FORGE gives slightly tighter constraints on log10 | fR 0 | . Constraints on cosmological parameters are very consistent among the three di ff erent nonlinear models. The Fisher Matrix forecast using the fitting formula is also very consistent with the MCMC results. We note that errors from 3 × 2pt analysis in the pessimistic setting are comparable to or better than WL alone in the optimistic setting. This demonstrates the strength of the 3 × 2pt analysis although we should bear in mind the limitation of the bias model used in this analysis.  \nWe 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 log10 | fR 0 | 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 log10 | fR 0 | also becomes slightly larger, with a value of B 1D = 0 . 667.', '5.1. Adding baryonic effects using BCemu': 'In this section, we study the impact of baryonic e ff ects on the constraints on the f ( R ) parameter. We use the seven-parameter emulator of baryonic e ff ects called BCemu (Giri & Schneider 2021) and we assume that the baryonic e ff ects and the modified gravity e ff ects 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 e ff ects. 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  \nfor these parameters at z = 0 to keep them within the prior range and we may miss important degeneracies between baryonic parameters and log10 | fR 0 | . We found that two baryonic parameters, M Giri etal . c and θ ej, have the strongest e ff ects 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 Giri et . al . c 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 Giri et . al . c . 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 baryon parameters. We will use the dimensionless parameter M c that is related to the one defined in Giri & Schneider (2021) via M c B M Giri et . al . c / M ⊙ where M ⊙ is the solar mass.', '5.2. Effects of adding baryons': 'Figure 6 shows the e ff ect of changing log10 | fR 0 | , h and two baryon parameters on the ratio of the WL angular power spectrum to the Λ CDM one. Baryonic e ff ects introduce scaledependent modifications to the angular power spectrum that are similar to the e ff ect 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.  \nFigure 7 shows the 2D contours of the WL constraints on parameters for the optimistic settings with baryons and Table 14 summarises constraints on log10 | fR 0 | .  \nIn 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 log10 | fR 0 | : the 1 σ error becomes 0.281 while it was 0.088 without baryons. However, it is important to stress that log10 | fR 0 | 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 highk tail of the matter power spectrum is not washed out by the baryonic feedback, and we can still distinguish between the e ff ect of f ( R ) and baryons.  \nNevertheless, if FORGE is used as the data, the di ff erence between FORGE and e-Mantis in the matter power spectrum can be absorbed by the shift in log10 M c and h , and the constraint on log10 | fR 0 | is shifted along the degeneracy direction between log10 | fR 0 | and log10 M c as well as h . This can be understood from Fig. 6. The di ff erence of the scale dependence between e-Mantis and FORGE at ℓ > 100 can be adjusted by decreasing log10 | fR 0 | and decreasing log10 M c. This leads to a lower amplitude that can be adjusted by decreasing h . Due to the combination of these e ff ects, the best-fit log10 | fR 0 | 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 8 . However, this dependence of h on Ξ ( k , z ) is very weak. To compute the power spectrum in f ( R ) we use Pf ( 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 ).  \nAs a result, we obtain a 95.5% upper limit of log10 | fR 0 | as log10 | fR 0 | < -5 . 477, which is incompatible with the input value of log10 | fR 0 | = -5 . 301. This is partly due to prior volume e ff ects shifting the contour to lower values of log10 | fR 0 | , caused by the strong degeneracy between log10 | fR 0 | and log10 M c as well as h .  \nTo confirm this, we obtain a profile likelihood for log10 | fR 0 | by fixing log10 | fR 0 | and finding a minimum χ 2 by varying other  \nTable 10. Additional parameters that are varied in the 3 × 2pt analysis.  \nFig. 5. Constraints from 3 × 2pt 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 .  \n<!-- image -->  \n|  \n|  \n<!-- image -->  \n|  \n|  \nTable 11. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the pessimistic setting where the same prediction is used to create the data and the model.  \nparameters. The profile likelihood was obtained by PROSPECT (Holm et al. 2023, /github ). This is shown in Fig. 8. The global best-fit value for log10 | fR 0 | is obtained as log10 | fR 0 | = -5 . 765, which is larger than the mean of the 1D marginalised posterior log10 | fR 0 | = -6 . 087. Figure 6 also shows the prediction of e-Mantis with the best-fit values, which agrees well with FORGE . The ∆ χ 2 for the input value of log10 | fR 0 | = -5 . 301 is found as ∆ χ 2 = 2 . 657, thus it is still within 2 σ . We also observe that the ∆ χ 2 curve becomes flat for smaller values of log10 | fR 0 | .  \nWe also tested the fitting formula as a model to fit to the data generated by FORGE . In this case, the posterior is highly nonGaussian and the chains did not converge well.  \nThis implies that we can break the degeneracy between f ( R ) gravity and baryonic e ff ects, 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 log10 M c from baryonic physics, we could break the degeneracy between log10 | fR 0 | and log10 M c. For example, in the case where FORGE was used as the data and they were modelled by e-Mantis , the mean of log10 M c is significantly shifted to a smaller value as we can see in Fig. 7. If we had a prior on log10 M c to prevent this, this bias could be avoided. To test this idea, we also ran an analysis imposing a prior on log10 M c and  \nTable 12. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the pessimistic setting where the data are generated by FORGE and it is fitted by di ff erent nonlinear models.  \n.  \n.  \n.  \nshowed the result in Fig. 7. We used a Gaussian prior with the width of 0 . 2, which was estimated from the weak lensing informed gas and stellar mass fraction measurements of massive haloes by Grandis et al. (2024). This prior is consistent with the error on log10 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 log10 M c and log10 | fR 0 | 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.  \nWe assumed that the e ff ect 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 e ff ect of baryonic and non-baryonic suppression mechanisms is greater than the sum of its parts for decaying dark matter (Elbers et al. 2024). The e ff ect of this coupling on the degeneracy needs to be studied in the future.  \nNext, 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 log10 | fR 0 | . In this case, due to the weaker constraining power, the 99 . 7% confidence level upper bound on log10 | fR 0 | reaches the prior boundary of log10 | fR 0 | = -7. Also, the 95 . 5% confidence level upper bound of θ ej is bounded by the prior. If the  \nFig. 6. The e ff ect of log10 | fR 0 | , two baryonic parameters log10 M c 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 | fR 0 | = -5 . 765 , 100 h = 65 . 736 , log10( M c) = 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 Λ CDMpower spectrum C Λ CDM ℓ is fixed to be the fiducial ones when we vary 100 h .  \n<!-- image -->  \nTable 13. Fiducial values and priors for BCemu parameters.Table 14. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the WL optimistic setting with baryons.  \n.  \n.  \n.  \ndata are generated by FORGE , strong degeneracies appear among log10 | fR 0 | and two baryonic parameters. The input log10 | fR 0 | is consistent within 1 σ for the fitting formula and e-Mantis although the prior bound log10 | fR 0 | = -7 is reached at 95.5% con-  \nTable 15. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the WL pessimistic setting with baryons.  \n.  \n.  \n.  \nfidence level lower bound. On the other hand, for ReACT , the input value is outside the 95.5% confidence level upper bound.  \nFigure 10 shows the 2D contours of the WL constraints on parameters for the pessimistic settings where Λ CDM is used  \n<!-- image -->  \n|  \n|  \n<!-- image -->  \n|  \n|  \nFig. 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 M c with the width of 0.2 in the case of FORGE fiducial data. Left panel : Constraints on baryon parameters. Right panel : Constraints on cosmological parameters.  \n<!-- image -->  \n10  \n|  \n|  \nFig. 8. The profile likelihood for log10 | fR 0 | 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 | fR 0 | = 5 . 301. We see that the 1D marginalised posterior is shifted to smaller values of log10 | fR 0 | .  \nAs we already observed without baryons, h is biased to a lower value.  \nWe also checked the case where fiducial data in 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 log10 | fR 0 | is now consistent with the input value within 1 σ .', '5.3. Constraints on | fR 0 | with the Λ CDM data': 'Based on these results, we consider the pessimistic setting and use e-Mantis to obtain an upper limit of log10 | fR 0 | in the presence of baryonic e ff ects to be conservative.  \nto 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 log10 | fR 0 | as  \nlog10 | fR 0 | < -5 . 21 . (30)  \nWe note that this bound depends on the prior log10 | fR 0 | = -7. To obtain the prior-independent bound, we follow the approach presented in Gordon & Trotta (2007), Piga et al. (2023) and Kou et al. (2023). The ratio of the marginalised posterior and prior is given by  \nb ( x | d , p ) = P ( x | d , p ) p ( x ) , (31)  \nwhere x is the parameter we are interested in (i.e. log10 | fR 0 | ), d is the data and p is the prior, and P is the Posterior. The Bayes factor B ( x 1 , x 2), which quantifies the support of the models with x = x 1 over the models with x = x 2 is given by  \nB ( x 1 , x 2) = b ( x 1 | d , p ) b ( x 2 | d , p ) = L ( d | x 1) L ( d | x 2) , (32)  \nwhere L ( d | x ) is the marginalised likelihood of the data for parameter x . For our purpose, we choose x = log10 | fR 0 | and fix x 1 to be the upper bound of the prior, x 1 = -7. Following Gordon &Trotta (2007), we then use B ( x 1 , x 2) = 2 . 5 to find x 2 so that the model with x = x 1 is favoured compared to the model with x = x 2 at 95 . 5% confidence level. This gives the 95 . 5% confidence interval of log10 | fR 0 | that does not depend on the prior. We note that this method applies only to a L ( d | x ) that is a monotonic function of x . Employing this technique, we obtain  \nlog10 | fR 0 | < -5 . 58 . (33)  \nFinally, in Fig. 11, we present the constraints on Ω m and log10 | fR 0 | 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  \n|  \n|  \nFig. 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 baryon parameters. Right panel : Constraints on cosmological parameters.  \n<!-- image -->  \n<!-- image -->  \n|  \n|  \ncomparison due to various di ff erences in the setting. Nonetheless, our constraint is comparable to the one presented in Schneider et al. (2020) [log10 | fR 0 | < -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 e ff ects are not included in their analysis, but again the constraint is comparable [log10 | fR 0 | < -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 log10 | fR 0 | when FORGE is used as data and this might be one of the reasons for this di ff erence. 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.', '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 will estimate the theoretical errors using ReACT and FORGE as these two give the most discrepant results. We will apply this method to ReACT for the pessimistic case of WL to check if we can 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 × 2pt likelihood, although we only show the result for the WL analysis in this paper. Here we will discuss the idea behind the formulation. We calculate 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  \nthe nonlinear matter power spectrum, F [ P δδ ], we can propagate the error on the power spectrum, ∆ P δδ , and find  \nF [ P δδ + ∆ P δδ ] = F [ P δδ ] + F [ ∆ P δδ ] C C XY i j ( ℓ ) + E XY i j ( ℓ ) , (34)  \nwhere we have defined the angular power spectrum error E XY i j ( ℓ ). Following Audren et al. (2013), we define the relative error on the power spectrum  \n∆ P δδ ( k , z ) B α ( k , z ) P δδ ( k , z ) . (35)  \nWe take the most conservative approach and assume that the error is uncorrelated between di ff erent values of ℓ . To account for this uncorrelated theoretical error, we add a new nuisance parameter εℓ for each multipole and define the shifted covariance as  \n˜ C XY i j ( ℓ ) B ˆ C XY i j ( ℓ ) + εℓ L 1 / 2 E XY i j ( ℓ ) . (36)  \nThe 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  \n˜ χ 2 ( εℓ ) = ℓ max X ℓ = ℓ min        (2 ℓ + 1) f sky        ˜ d mix ℓ ( εℓ ) ˜ d th ℓ ( εℓ ) + ln ˜ d th ℓ ( εℓ ) d obs ℓ -N        + ε 2 ℓ        . (37)  \nThe quantities ˜ d th and ˜ d mix are constructed in the same way as d th and d mix [see Eqs. (15) to (17)] using the shifted covariance matrix ˜ C XY i j ( ℓ ). To include the theoretical error, we vary εℓ and marginalise over them. To a very good approximation, this is equivalent to minimising the ˜ χ 2 with respect to εℓ at the level of the likelihood. Thus, we define our new χ 2 as the minimum of the ˜ χ 2 with respect to εℓ :  \nχ 2 B min εℓ ∈ R L ˜ χ 2 . (38)  \nThe normalisation factor L 1 / 2 can be explained as follows. If we were to measure ˆ C obs ( ℓ ) = ˆ C th ( ℓ ) + E ( ℓ ), the minimisation would  \nFig. 11. The same as in Fig. 10, but we show the constraint on the derived quantity Ω m.  \n<!-- image -->  \n|  \n|  \n<!-- image -->  \n|  \n|  \nFig. 10. Constraints on parameters from the WL pessimistic setting with baryons where the data are created in Λ CDMandit is fitted by e-Mantis . Left panel : Constraints on baryon parameters. Right panel : Constraints on cosmological parameters.  \n<!-- image -->  \n|  \n|  \nfind εℓ = L -1 / 2 . The resulting χ 2 = P ℓ ε 2 ℓ = P ℓ L -1 = 1 would match our expectation that a one-sigma theoretical error for each ℓ results in an increase of the χ 2 by one.  \nThe 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 XY i j ( ℓ ) 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  \nd˜ χ d εℓ = (2 ℓ + 1) f sky           GLYPH<16> ˜ d mix ℓ GLYPH<17> ' + GLYPH<16> ˜ d th ℓ GLYPH<17> ' ˜ d th ℓ -˜ d mix ℓ GLYPH<16> ˜ d th ℓ GLYPH<17> ' GLYPH<16> ˜ d th ℓ GLYPH<17> 2           + 2 εℓ . (39)  \nThe derivatives of the determinants are computed using Jacobi's formula. This gives, for example  \nGLYPH<16> ˜ d th ℓ GLYPH<17> ' = det GLYPH<16> ˆ C th ℓ GLYPH<17> L 1 / 2 Tr GLYPH<20> GLYPH<16> ˆ C th ℓ GLYPH<17> -1 E ℓ GLYPH<21> , (40)  \nand a similar expression for GLYPH<16> ˜ d mix ℓ GLYPH<17> ' . We compute the inverse of the covariance numerically and obtain the second derivatives numerically from the first derivative by doing a double-sided threepoint 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 logarithmicallyspaced grid with 100 discrete values. The other values are obtained from an interpolating function. We checked that the e ff ect 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 e ff ects of changing other theory parameters to the observed power spectrum.", '6.3. Application to ReACT': 'We first apply this method to the case where 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 | fR 0 | 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  \nΞ ReACT-FORGE = Ξ ReACT Ξ FORGE Ξ ReACT ! fiducial . (41)  \nThe 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 a ff ect mainly log10 | fR 0 | 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 log10 | fR 0 | , Ω c h 2 and ln(10 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.  \n|  \n|  \n<!-- image -->  \n|  \n|  \n<!-- image -->  \n|  \n|  \nFig. 12. The e ff ect 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.Fig. 13. The e ff ect of adding theoretical errors in the presence of baryonic e ff ects. 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.  \n<!-- image -->  \nWe 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.  \nFinally, we study the impact of adding baryonic e ff ects 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 a ff ected significantly by baryonic e ff ects. 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, baryon parameters and log10 | fR 0 | . This leads to a tighter lower bound on log10 | fR 0 | compared with the case where the same data are fitted by e-Mantis itself. This reinforces our conclusion that the 1D marginalised constraint on log10 | fR 0 | is sensitive to degeneracies among parameters, and the  \ndi ff erence in the theoretical predictions strongly a ff ects the constraint.  \nThe result shown here is just an illustration of the inclusion of theoretical errors and their e ff ects 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 studied the e ff ect of using di ff erent 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 di ff erent models in the MontePython pipeline to predict angular power spectra for weak lensing (WL), photometric galaxy clustering (GCph) and their cross-correlation (XCph). Comparing with the N -body simulation data obtained in paper 1, we found that e-Mantis agreed very well with ECOSMOG that 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 di ff erence 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 .  \nWe used the fiducial value of | fR 0 | = 5 × 10 -6 (log10 | fR 0 | = -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 a ff ects the degeneracy between log10 | fR 0 | and cosmological parameters for large log10 | fR 0 | . e-Mantis gives more Gaussian constraints with smaller errors. When FORGE is used to create the data, the 1D mean of log10 | fR 0 | is not strongly biased in the case of e-Mantis and the fitting formula and the 1D bias is at most 0 . 6 σ . Even  \nfor the 3 × 2pt analysis including all the probes and their crosscorrelations, the 1D bias is 0 . 15 σ for log10 | fR 0 | in the case of e-Mantis .  \nReACT shows a large bias in log10 | fR 0 | 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 di ff erence 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 e ff ects with BCemu , errors were not significantly a ff ected and the input values were still recovered. However, the lower bound on log10 | fR 0 | is tighter than the one obtained by applying e-Mantis itself. This reinforces our conclusion that the 1D marginalised constraint on log10 | fR 0 | is sensitive to the degeneracies among parameters.  \nThe 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 baryon parameters are well constrained despite the degeneracies among these parameters. However, the di ff erence in the nonlinear dark matter prediction can be compensated by the adjustment of baryon parameters as well as cosmological parameters, and the 1D marginalised constraint on log10 | fR 0 | is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the Λ CDMsynthetic data for WL, we obtained the prior-independent bound of log10 | fR 0 | < -5 . 6 using e-Mantis .  \nBased on the result of this paper, we draw the following conclusions:  \n- -It is important to check the agreement of di ff erent N -body codes that are used to create theoretical predictions. This is 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 e ff ect 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.\n- -Including baryonic e ff ects can worsen the bias in the 1D marginalised constraint on log10 | fR 0 | . This is because the difference in the nonlinear dark matter power spectrum prediction can be compensated by the adjustment of baryon parameters, and the best-fit values are biased. In addition, the degeneracy between baryon parameters, cosmological parameters and log10 | fR 0 | leads to a stronger volume e ff ect. This bias can be avoided if there is a prior on baryon parameters from external data sets for example.\n- -We then need to check the accuracy of emulators. We found that FORGE su ff ers from larger emulation errors and this leads to a larger di ff erence between FORGE and e-Mantis compared with the di ff erence in their baseline N -body code MG-AREPO and ECOSMOG in some cosmologies. Improving the emulation technique will be able to make the agreement better.\n- -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  \nestimation 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.  \nThe 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 et al. (2023).  \nFinally, 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 Λ CDMand exotic cosmologies. A detailed study of this issue is beyond the scope of our paper, but we refer to Euclid Collaboration: Mellier et al. (2024) 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) incldued 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.  \nAcknowledgements. 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 Luftund 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 O ffi ce (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': 'Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv:astro-ph / 0609591  \n- Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869\n- Arnold, C., Leo, M., & Li, B. 2019, Nature Astron., 3, 945\n- Arnold, C. & Li, B. 2019, MNRAS, 490, 2507\n- Arnold, C., Li, B., Giblin, B., et al. 2022, MNRAS, 515, 4161\n- Audren, B., Lesgourgues, J., Bird, S., Haehnelt, M. G., & Viel, M. 2013, JCAP, 01, 026  \nBaldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75  \n- Bose, B., Carrilho, P., Marinucci, M., et al. 2023, arXiv:2311.13529\n- Bose, B., Cataneo, M., Tröster, T., et al. 2020, MNRAS, 498, 4650\n- Bose, B., Tsedrik, M., Kennedy, J., et al. 2023, MNRAS, 519, 4780\n- Bose, B., Wright, B. S., Cataneo, M., et al. 2021, MNRAS, 508, 2479\n- Bose, S., Li, B., Barreira, A., et al. 2017, JCAP, 02, 050\n- Brax, P., Casas, S., Desmond, H., & Elder, B. 2021, Universe, 8, 11\n- Brax, P., van de Bruck, C., Davis, A.-C., & Shaw, D. J. 2008, Phys. Rev. D, 78, 104021  \nBuchdahl, H. A. 1970, MNRAS, 150, 1  \nBurrage, C., March, B., & Naik, A. P. 2024, JCAP, 04, 004  \n- Carrilho, P., Carrion, K., Bose, B., et al. 2022, MNRAS, 512, 3691\n- Carrion, K., Carrilho, P., Spurio Mancini, A., Pourtsidou, A., & Hidalgo, J. C. 2024, arXiv:2402.18562  \nCasas, S., Cardone, V. F., Sapone, D., et al. 2023, arXiv:2306.11053  \nCasas, S., Lesgourgues, J., Schöneberg, N., et al. 2024, A&A, 682, A90  \nCataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121  \nCautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195  \nClifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1  \nDesmond, H. & Ferreira, P. G. 2020, Phys. Rev. D, 102, 104060  \nElbers, W., Frenk, C. S., Jenskins, A., et al. 2024, arXiv:2403.12967  \nEuclid Collaboration: Adamek et al. 2024, in prep  \nEuclid Collaboration: Blanchard, A., Camera, S., Carbone, C., et al. 2020, A&A,', '642, A191': 'Euclid Collaboration: Knabenhans, M. et al. 2021, MNRAS, 505, 2840  \n- Euclid Collaboration: Mellier, Y., Acevedo Barroso, J., Achúcarro, A., et al. 2024, A&A, accepted, arXiv:2405.13491  \nEuclid Collaboration: Racz et al. 2024, in prep  \nEuclid Collaboration: Sciotti, D., Gouyou Beauchamps, S., Cardone, V. F., et al. 2023, arXiv:2310.15731  \nFerreira, P. G. 2019, Ann. Rev. Astron. Astrophys., 57, 335  \nFiorini, B., Koyama, K., & Baker, T. 2023, JCAP, 12, 045  \nFrusciante, N. Pace, F., Cardone, V. F., et al. 2023, arXiv:2306.12368  \n- Gelman, A. & Rubin, D. B. 1992, Statist. Sci., 7, 457\n- Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813  \nGiri, S. K. & Schneider, A. 2021, JCAP, 12, 046  \n- Gordon, C. & Trotta, R. 2007, MNRAS, 382, 1859  \nGrandis, S., Arico\', G., Schneider, A., & Linke, L. 2024, MNRAS, 528, 4379  \nGuillet, T. & Teyssier, R. 2011, Journal of Computational Physics, 230, 4756  \nHarnois-Déraps, J., Hernandez-Aguayo, C., Cuesta-Lazaro, C., et al. 2023, MN-  \nRAS, 525, 6336 Holm, E. B., Nygaard, A., Dakin, J., et al. 2023, arXiv:2312.02972 Hu, W. & Sawicki, I. 2007, Phys. Rev., D76, 064004 Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90 Ishak, M. 2019, Living Rev. Rel., 22, 1 Joyce, A., Jain, B., Khoury, J., & Trodden, M. 2015, Phys. Rept., 568, 1 Khoury, J. & Weltman, A. 2004, Physical Review Letters, 93 Kou, R., Murray, C., & Bartlett, J. G. 2023, arXiv:2311.09936 Koyama, K. 2016, Rept. Prog. Phys., 79, 046902 Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970 Li, B., Zhao, G.-B., Teyssier, R., & Koyama, K. 2012, JCAP, 01, 051 McKay, M. D., Beckman, R. J., & Conover, W. J. 1979, Technometrics, 21, 239 Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 Mead, A., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 Mead, A., Peacock, J., Heymans, C., et al. 2015, MNRAS, 454, 1958 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 Perez, F. & Granger, B. E. 2007, Computing in Science Engineering, 9, 21 Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 Piga, L., Marinucci, M., D\'Amico, G., et al. 2023, JCAP, 04, 038 Pisanti, O., Mangano, G., Miele, G., & Mazzella, P. 2021, JCAP, 04, 020 Planck Collaboration: Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020,  \nA&A, 641, A6, [Erratum: Astron.Astrophys. 652, C4 (2021)]  \nPuchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348  \n- Ramachandra, N., Valogiannis, G., Ishak, M., & Heitmann, K. 2021, Phys. Rev. D, 103, 123525  \nRasmussen, C. E. & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (The MIT Press)  \nRiess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009  \nSchneider, A., Refregier, A., Grandis, S., et al. 2020, JCAP, 04, 020  \nSotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451  \nSpringel, V. 2010, MNRAS, 401, 791  \nSpurio Mancini, A. & Bose, B. 2023, arXiv:2305.06350  \nSpurio Mancini, A., Köhlinger, F., Joachimi, B., et al. 2019, MNRAS, 490, 2155  \nSpurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771 Sáez-Casares, I., Rasera, Y., & Li, B. 2023, MNRAS, 527, 7242 Takahashi, R., Sato, M., Nishimichi, T., et al. 2012, ApJ, 761, 152 Teyssier, R. 2002, A&A, 385, 337 Tröster, T., Asgari, M., Blake, C., et al. 2021, A&A, 649, A88 Tsedrik, M., Bose, B., Carrilho, P., et al. 2024, arXiv:2404.11508 Valogiannis, G. & Bean, R. 2017, Phys. Rev. D, 95, 103515 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engg., 13, 22 Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32 Winther, H., Casas, S., Baldi, M., et al. 2019, Phys. Rev. D, 100, 123540 Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.-B. 2017, JCAP, 08, 006  \nWinther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208  \n- 1 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK\n- 3 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK\n- 2 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany\n- 4 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France\n- 6 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal\n- 5 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal  \n7  \nRuhr University Bochum, Faculty of Physics and Astronomy, As-  \ntronomical Institute (AIRUB), German Centre for Cosmological  \nLensing (GCCL), 44780 Bochum, Germany  \n- 9 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93 / 3, 40129 Bologna, Italy\n- 8 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany\n- 10 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 12 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 11 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 13 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy\n- 15 Institut universitaire de France (IUF), 1 rue Descartes, 75231 PARIS CEDEX 05, France\n- 14 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK\n- 16 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany\n- 18 Université St Joseph; Faculty of Sciences, Beirut, Lebanon\n- 17 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France\n- 19 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway\n- 21 Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland\n- 20 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK\n- 22 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93 / 2, 40129 Bologna, Italy\n- 24 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s / n, 08193 Barcelona, Spain\n- 23 INFN-Sezione di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 25 Institut de Ciencies de l\'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s / n Cerdanyola del Vallés, 08193 Barcelona, Spain  \n- 26 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA\n- 28 School of Mathematics and Physics, University of Surrey, Guildford, Surrey, GU2 7XH, UK\n- 27 Institut de Physique Théorique, CEA, CNRS, Université ParisSaclay 91191 Gif-sur-Yvette Cedex, France\n- 29 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy\n- 31 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy\n- 30 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy\n- 32 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy\n- 34 Institut d\'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France\n- 33 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy\n- 35 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany\n- 37 INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 36 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 38 Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126, Napoli, Italy\n- 40 INFN section of Naples, Via Cinthia 6, 80126, Napoli, Italy\n- 39 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy\n- 41 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal\n- 43 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France\n- 42 Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal\n- 44 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy\n- 46 Port d\'Informació Científica, Campus UAB, C. Albareda s / n, 08193 Bellaterra (Barcelona), Spain\n- 45 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain\n- 47 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy\n- 49 Instituto de Astrofísica de Canarias, Calle Vía Láctea s / n, 38204, San Cristóbal de La Laguna, Tenerife, Spain\n- 48 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 50 European Space Agency / ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy\n- 52 Université Claude Bernard Lyon 1, CNRS / IN2P3, IP2I Lyon, UMR 5822, Villeurbanne, F-69100, France\n- 51 ESAC / ESA, Camino Bajo del Castillo, s / n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain\n- 53 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland\n- 55 Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain\n- 54 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, 08028 Barcelona, Spain\n- 56 UCB Lyon 1, CNRS / IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France\n- 58 Université Paris-Saclay, CNRS, Institut d\'astrophysique spatiale, 91405, Orsay, France\n- 57 Department of Astronomy, University of Geneva, ch. d\'Ecogia 16, 1290 Versoix, Switzerland\n- 59 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy\n- 61 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France\n- 60 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy\n- 62 Institut d\'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain\n- 64 INAF-Osservatorio Astronomico di Padova, Via dell\'Osservatorio 5, 35122 Padova, Italy\n- 63 FRACTAL S.L.N.E., calle Tulipán 2, Portal 13 1A, 28231, Las Rozas de Madrid, Spain\n- 65 Universitäts-Sternwarte München, Fakultät für Physik, LudwigMaximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany\n- 67 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK\n- 66 Dipartimento di Fisica "Aldo Pontremoli", Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy\n- 68 Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany\n- 70 Cosmic Dawn Center (DAWN), Denmark\n- 69 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark\n- 71 Université Paris-Saclay, CNRS / IN2P3, IJCLab, 91405 Orsay, France\n- 73 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA\n- 72 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany\n- 74 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 76 Aix-Marseille Université, CNRS / IN2P3, CPPM, Marseille, France\n- 75 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland\n- 77 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, CH1211 Genève 4, Switzerland\n- 79 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland\n- 78 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland\n- 80 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands\n- 82 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy\n- 81 Centre de Calcul de l\'IN2P3 / CNRS, 21 avenue Pierre de Coubertin 69627 Villeurbanne Cedex, France\n- 83 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c / o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy\n- 85 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK\n- 84 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93 / 2, 40129 Bologna, Italy\n- 86 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France\n- 88 Institut d\'Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France\n- 87 University of Applied Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210 Windisch, Switzerland\n- 89 Institut de Física d\'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain\n- 91 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark\n- 90 European Space Agency / ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands\n- 92 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 94 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada\n- 93 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 95 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy\n- 96 Centre National d\'Etudes Spatiales - Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France  \n- 97 Institute of Space Science, Str. Atomistilor, nr. 409 M˘agurele, Ilfov, 077125, Romania\n- 98 Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy\n- 99 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile\n- 100 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25 / 8, 6020 Innsbruck, Austria\n- 101 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain\n- 102 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal\n- 103 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain\n- 104 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 105 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands\n- 106 Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova, via Dodecaneso 33, 16146, Genova, Italy\n- 107 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA\n- 108 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy\n- 109 Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020, Nus (Aosta Valley), Italy\n- 110 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK\n- 111 School of Physics and Astronomy, Cardi ff University, The Parade, Cardi ff , CF24 3AA, UK\n- 112 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France\n- 113 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy\n- 114 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain\n- 115 CERCA / ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA\n- 116 Departamento de Física Fundamental. Universidad de Salamanca. Plaza de la Merced s / n. 37008 Salamanca, Spain\n- 117 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 118 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 119 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 120 Ludwig-Maximilians-University, Schellingstrasse 4, 80799 Munich, Germany\n- 121 Max-Planck-Institut für Physik, Boltzmannstr. 8, 85748 Garching, Germany\n- 122 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA\n- 123 Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands\n- 124 Université Côte d\'Azur, Observatoire de la Côte d\'Azur, CNRS, Laboratoire Lagrange, Bd de l\'Observatoire, CS 34229, 06304 Nice cedex 4, France\n- 125 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA\n- 126 Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA\n- 127 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary\'s University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada\n- 128 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain  \n- 129 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica, Universidad de La Laguna (ULL), 38200, La Laguna, Tenerife, Spain\n- 131 CEA Saclay, DFR / IRFU, Service d\'Astrophysique, Bat. 709, 91191 Gif-sur-Yvette, France\n- 130 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK\n- 132 Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00 076, Finland\n- 134 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000, Grenoble, France\n- 133 Instituto de Astrofísica de Canarias, c / Via Lactea s / n, La Laguna E-38200, Spain. Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna, E-38200, Spain\n- 135 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland\n- 137 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia\n- 136 Serco for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 138 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, Victoria 3122, Australia\n- 140 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town, 7535, South Africa\n- 139 School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK\n- 141 ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil\n- 143 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm, SE-106 91, Sweden\n- 142 IRFU, CEA, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France\n- 144 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK\n- 146 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy\n- 145 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125, Firenze, Italy\n- 147 Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 149 HE Space for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 148 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal\n- 150 Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy\n- 152 Department of Physics, Royal Holloway, University of London, TW20 0EX, UK\n- 151 Theoretical astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden\n- 153 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA\n- 155 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark\n- 154 Cosmic Dawn Center (DAWN)\n- 156 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA\n- 157 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA  \n|  \n|  \nFig. A.1. Bias due to di ff erent nonlinear modelling. The synthetic data are created by e-Mantis and fitted by two di ff erent models in WL pessimistic case without baryons.  \n<!-- image -->  \nFig. A.2. Bias due to di ff erent nonlinear modelling. The synthetic data are created by e-Mantis and fitted by two di ff erent models in WL pessimistic case with baryons.  \n<!-- image -->  \n|  \n|  \nTable A.1. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the WL pessimistic setting.  \n.  \n.  \n.', '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 σ .  \nTable A.2. The mean, standard deviation and 68.3% lower and upper limit of log10 | fR 0 | for the WL pessimistic setting with baryons.  \n.  \n.  \n.', 'Appendix B: Construction of the relative error': '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  \nα ( k , z ) = ∆ P δδ ( k , z ) P δδ ( k , z ) = ∆Ξ ( k , z ) Ξ ( k , z ) . (B.1)  \nTo 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 di ff erence 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 twohalo power spectra separately. HMcode2020 predictions for the pseudo-cosmology non-linear 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 di ff erence 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 di ff erent 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.  \nFinally, we construct α using only the di ff erence of the predictions at a fiducial of | fR 0 | = 5 × 10 -6 in order not to bias our results. Our final fit has the form  \nα ( k , z ) = A ( z ) x 2 + x x 2 + x + 1 , (B.2)  \nA ( z ) = A 1 exp GLYPH<16> z -A 2 A 3 GLYPH<17> + 1 + A 4 , (B.3)  \nk p( z ) = B 1 exp " tanh z -B 2 B 3 !# , (B.4)  \nwhere we use x B 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 di ff er by the theoretical error.  \nFig. B.1. Relative di ff erence 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 are all computed at the fiducial cosmology.  \n<!-- image -->  \nTable B.1. Best fit values of the relative error function.  \n3'}