On the degeneracies between baryons, massive neutrinos and f(R) gravity in Stage IV cosmic shear analyses

Modelling nonlinear structure formation is essential for current and forthcoming cosmic shear experiments. We combine the halo model reaction formalism, implemented in the REACT code, with the COSMOPOWER machine learning emulation platform, to develop and publicly release REACTEMU-FR, a fast and accurate nonlinear matter power spectrum emulator for $f(R)$ gravity with massive neutrinos. Coupled with the state-of-the-art baryon feedback emulator BCEMU, we use REACTEMU-FR to produce Markov Chain Monte Carlo forecasts for a cosmic shear experiment with typical Stage IV specifications. We find that the inclusion of highly nonlinear scales (multipoles between $1500\leq \ell \leq 5000$) only mildly improves constraints on most standard cosmological parameters (less than a factor of 2). In particular, the necessary modelling of baryonic physics effectively damps most constraining power on the sum of the neutrino masses and modified gravity at $\ell \gtrsim 1500$. Using an approximate baryonic physics model produces mildly improved constraints on cosmological parameters which remain unbiased at the $1\sigma$-level, but significantly biases constraints on baryonic parameters at the $>2\sigma$-level.


INTRODUCTION
One of the main goals of forthcoming Stage IV cosmic shear surveys, such as Euclid (Laureĳs et al. 2011) 1 , the Vera Rubin Observatory's Legacy Survey of Space and Time (VRO-LSST, Ivezić et al. 2019) 2 and the Nancy Grace Roman Telescope (Spergel et al. 2015) 3 , is to perform precise and accurate tests of gravity on cosmological scales.Constraints on deviations from General Relativity, the gravity theory underlying our standard LCDM model of cosmology, will particularly benefit from measurements of the cosmic shear power spectrum from unprecedentedly large numbers of galaxies.Contributions to the cosmic shear power spectrum from small, nonlinear scales, will be of the utmost importance to distinguish between competing gravity models with different nonlinear screening mechanisms (see e.g.Koyama 2016;Harnois-Déraps et al. 2022).
On the theoretical end, the achievement of this goal hinges on various factors.We list some of the primary criteria below: (i) we require accurate prescriptions for the key observables of Stage IV surveys, at both linear and nonlinear scales; (ii) we require these prescriptions to be computationally efficient, to use them within rigorous Bayesian analysis pipelines; ★ E-mail: a.spuriomancini@ucl.ac.uk 1 https://www.euclid-ec.org/ 2 https://www.lsst.org/ 3 https://roman.gsfc.nasa.gov/(iii) we need to theoretically model all relevant and known physical phenomena in order to distinguish any potential new physics.So far, these criteria have only been met for a subset of theoretical descriptions of gravity and cosmology, and all suffer from various limitations.Emulators trained on -body simulations (e.g.Lawrence et al. 2017;Knabenhans et al. 2019Knabenhans et al. , 2021;;Angulo et al. 2021;Ramachandra et al. 2021;Arnold et al. 2021) are restricted to the parameter space probed by their training simulations, and cannot be safely extended without running new simulations.Analytic prescriptions (e.g.Takahashi et al. 2012;Zhao 2014;Mead et al. 2016;Winther et al. 2019;Mead et al. 2020) are also restricted to the simulations which they have originally been fit to.
A method potentially capable of overcoming some of these issues is the halo model reaction, a theoretically flexible framework introduced in Cataneo et al. (2019) (see also Mead 2017, for initial motivations) based on 1-loop perturbation theory and the halo model (see e.g.Asgari et al. 2023, for a recent review), shown to be consistent with -body simulations at the ∼ 2% level down to  ∼ 3ℎ/Mpc for a wide range of non-standard cosmological scenarios (Srinivasan et al. 2021;Bose et al. 2021;Carrilho et al. 2022;Parimbelli et al. 2022;Adamek et al. 2022).However, this method is not as accurate as -body based emulators (see e.g.Arnold et al. 2021), nor is the associated code REACT (Bose et al. 2020(Bose et al. , 2022, ) , ) as fast as current emulators.
This being said, the sacrifice in accuracy of the halo model reaction may be a non-issue in real data analyses, once we take into account relevant physical effects such as baryon feedback processes and massive neutrinos.Furthermore, the computational inefficiency of RE-ACT can be circumvented by employing the recent COSMOPOWER machine learning platform (Spurio Mancini et al. 2022, ), which trains fast neural network emulators of key cosmological quantities starting from a set of 'slow' theoretical predictions.Once trained, these emulators can be used in Bayesian inference pipelines to replace the call to the original, computationally intensive prescription, massively accelerating parameter estimation.Developing similar accelerated pipelines is critical to the success of forthcoming Stage IV cosmic shear analyses, as these will have to be carried out over unprecedentedly large parameter spaces for accurate modelling of baryons, massive neutrinos, systematics and deviations from LCDM.
In this paper we present REACTEMU-FR, an emulator of the nonlinear effects on the matter power spectrum caused by massive neutrinos in Hu-Sawicki  () gravity (Hu & Sawicki 2007), a popular non-standard model of gravity being considered by Stage IV surveys (e.g.Amendola et al. 2018;Ramachandra et al. 2021).This emulator is constructed by COSMOPOWER using REACT predictions, demonstrating the ease by which non-standard model emulators can be generated using the two codes, and seamlessly integrated into analysis pipelines.We use this emulator to forecast the constraining power on cosmology and gravity, as well as to investigate parameter degeneracies using cosmic shear, in the context of Stage IV surveys.We check the capability of such surveys to detect new physics given a full host of secondary physical phenomena including baryon effects and massive neutrinos.Conversely, we also investigate the importance of including all non-standard physical effects in order to remain unbiased in the final parameter constraints.
This paper is outlined as follows.In section 2 we present and validate REACTEMU-FR and discuss baryonic effects and cosmic shear.In section 3 we perform forecasts for a Stage IV survey using different sets of degrees of freedom.We summarise our findings and conclude in section 4.

Massive neutrinos and 𝑓 (𝑅) gravity: REACTEMU-FR
Massive neutrinos entered the regime of standard physics two decades ago with the measurements of flavour oscillations (Fukuda et al. 1998;Ahmed et al. 2004), and were shown to have a significant effect on cosmological structure formation (Lesgourgues & Pastor 2006).They are usually parameterised in terms of the combined neutrino mass,   .Neutrino oscillation experiments have placed a lower bound on this parameter of   ≥ 0.06 eV.There is also an upper bound from the KATRIN collaboration of   ≤ 0.8 eV (Aker et al. 2022).
Modified gravity theories are extensions of General Relativity developed to deviate from standard gravity on cosmological scales, where it remains largely untested.These models can provide interesting solutions to fundamental physical problems (see e.g.Bernardo et al. 2023, for a recent review).In this paper we select a specific theory of modified gravity -the popular Hu-Sawicki  () model (Hu & Sawicki 2007).In the  () class of models, the Einstein-Hilbert action is generalised to include an arbitrary function of the scalar curvature, .The Hu-Sawicki form is given by where we choose  = 1 and  2 = Ω m  2 0 . 1 and  2 are parameters of the theory, and as is commonly done, we reparametrise these in terms of the the value of  R =  ()/ today R0 being the background scalar curvature today.In the following we take the effective LCDM-limit to be |  R 0 | = 10 −9 .We have checked that this value induces effects in the observable of interest much smaller than all other sources of error, making it indistinguishable from the true LCDM limit of  R 0 = 0.This model predicts non-trivial scale dependencies in the growth of structure.It has been shown to have degeneracies with massive neutrinos (Mead et al. 2016;Wright et al. 2019).The degree of deviation from LCDM in this model is parameterised by the value of the additional scalar degree of freedom today,  R 0 , with the LCDM limit of the theory given by  R 0 → 0. Current cosmological bounds on  R 0 are of O (10 −5 ) (Cataneo et al. 2015;Lombriser 2014;Desmond & Ferreira 2020;Brax et al. 2021).
We train an emulator for the  () + -boost, i.e. the nonlinear modification to a LCDM spectrum caused by  () and massive neutrinos, as a function of Fourier mode  and scale factor : where  is the scale factor and   ()+ NL and  ΛCDM NL are the nonlinear matter power spectra in the  () with massive neutrino cosmology and LCDM model, respectively.The nonlinear power spectrum under the halo model reaction prescription (Cataneo et al. 2019) is given by R is the halo model reaction and  pseudo NL is the pseudo power spectrum.The latter is defined as originating from a LCDM cosmology whose initial conditions are tuned such that the linear clustering matches the non-standard cosmology at a target redshift.In practice, we can construct this by supplying prescriptions such as HALOFIT (Takahashi et al. 2012) or HMCODE (Mead et al. 2020) with a modified linear power spectrum (see Bose et al. 2020, for details).
R encodes nonlinear corrections to the pseudo power spectrum due to  () and massive neutrinos -in essence, it is a ratio between halo model predictions for the non-standard cosmology and a pseudo cosmology, with a correction factor coming from a calibration with 1-loop perturbation theory at quasi-nonlinear scales.We refer the reader to Bose et al. (2021) for more details, and we follow this reference in our modelling of R. Note that the background cosmology is assumed to be LCDM, with the effects of modified gravity and massive neutrinos only entering at the level of the perturbations.
In this paper we use HMCODE (Mead et al. 2020) to model the pseudo and LCDM power spectrum appearing in Equation 3and Equation 4. To compute the halo model reaction we use the REACT code (Bose et al. 2022).Our COSMOPOWER emulator, REACTEMU-FR, is then built on Equation 3 computed within the parameter ranges in Table 1, sampled using a Latin Hypercube.The ranges are chosen as the rough 5 bounds around the Planck 2018 (Aghanim et al. 2020) best fit values including a varying   .We take an upper bound on  R 0 to be the 3 limit found in the cosmic shear forecast of Bose et al. (2020) using only linear scales.These represent very conservative ranges, with any deviation outside these being in strong disagreement with current observations.Note that the emulator can always be trained on extended ranges as the halo  2020), and we model the boost, which is far more accurate as the ratio of HMCODE predictions reduces systematics significantly (Cataneo et al. 2019).Given that we wish to primarily investigate the bias incurred from incomplete baryonic physics modelling, we believe this prescription is sufficiently accurate for our purposes.REACTEMU-FR is deemed to have an accuracy of ∼ 3% at scales of  ≤ 1ℎ/Mpc within its prior range as shown by comparisons with emulators and -body simulations (Cataneo et al. 2019;Bose et al. 2021;Arnold et al. 2021), with worsening accuracy for larger deviations from LCDM.For reference, for |  R 0 | = 10 −6 and   ≤ 0.1, the boost was shown to be 2% consistent with -body measurements (Bose et al. 2021;Arnold et al. 2021).
We train a COSMOPOWER emulator on ∼ 10 5 REACT predictions and test its performance on ∼ 2 • 10 4 testing samples, finding ≤ 0.05% residuals across the entire emulated -range  ∈ [0.01, 3]ℎ/Mpc, i.e. the entire validity range of REACT. Figure 1 shows a quantile plot of the accuracy of our emulator on the testing set.This accuracy level is more than an order of magnitude bigger than that used in Spurio Mancini et al. (2022) to validate their matter power spectrum emulators, which were shown to be comfortably accurate to guarantee unbiased emulated Stage IV posterior contours.The speed up factor provided by this emulation is above O (10 4 ) as the standard boost computation requires multiple Boltzmann solver calls, multiple HMCODE calls and a call to REACT, with a typical boost taking up to ∼ 25 s to compute compared to REACTEMU-FR's ∼ 1 ms.

Baryons: BCEMU
Beyond massive neutrinos and modified gravity, we must also consider the effects of baryons as we move into the nonlinear regime.These have been shown to have a significant impact on the matter power spectrum and the cosmic shear spectrum (Semboloni et al. 2011;Mead et al. 2020;Martinelli et al. 2021;Schneider et al. 2020a,b;Aricò et al. 2021).
BCEMU (Giri & Schneider 2021, ) is a state-of-the-art baryonic boost,  baryons , emulator giving the nonlinear effects of baryons on the matter power spectrum.It requires a set of cosmology-independent parameters as input, augmented with the baryon fraction  b = Ω b /Ω m .Some BCEMU parameters have also been shown to have a time evolution to varying degrees, depending on which hydrodynamical simulation they are fit to.In Giri & Schneider (2021) the authors propose a power law dependence where   is a free parameter giving the time dependence of the BCEMU baryon input parameter   , with the value today,  ,0 =   ( = 1), also being a free parameter.Giri & Schneider (2021) investigate two parameterisations: • a 7-parameter model,   ∈ {log 10  c ,   ,  ej , , , , }; • a 3-parameter model,   ∈ {log 10  c ,   ,  ej }.
We refer to Giri & Schneider (2021) for an in-depth discussion of the meaning of each of these parameters.We note that the 3-parameter model is motivated by its ability to fit within 1% a large suite of hydrodynamical simulation measurements at a fixed redshift.One of our goals (see section 3) is to investigate to what degree this reduced parameterisation improves or biases cosmological inference.Note that an -parameter model actually has  × 2 degrees of freedom when including the time dependence.
In Table 1 we list the chosen BCEMU parameters along with their fiducial values today ( = 1) and prior ranges that we select for our analysis.The ranges are 20% tighter than the defaults implemented in the BCEMU software.The reason for this is to prevent any deviations outside the BCEMU allowed ranges when we consider higher redshifts, given the time dependence we impose.These ranges are also the edges of the uniform prior distributions used for these baryon parameters in our Markov Chain Monte Carlo (MCMC) analyses of section 3. Priors on the power law parameters   are given in Table 2, which are also chosen so as not to exceed the allowed BCEMU range at any redshift considered in this work (0 ≤  ≤ 5).As fiducial values for these parameters, we take  log 10   = 0.01 and    = 0.06, with all other   = 0.

Cosmic Shear
We model the nonlinear matter power spectrum  NL as where nonlinear LCDM effects are modelled using a COS-MOPOWER emulator of HMCODE (Mead et al. 2016), the baryonic effects are modelled using BCEMU and the massive neutrino and  ()-gravity effects by REACTEMU-FR.The cosmic shear power spectrum     (ℓ) can be obtained from the nonlinear matter power spectrum  NL (, ) as an integral along the line of sight (e.g.Bartelmann & Schneider 2001;Kilbinger 2015;Kilbinger et al. 2017;Kitching et al. 2017;Asgari et al. 2021): where  is the comoving distance as a function of the scale factor,  H is the Hubble radius  H = / 0 , and    are weighting functions for each redshift bin , where   is the redshift distribution for bin ,  0 is the Hubble constant and Ω m is the total matter density parameter.Note that throughout the paper we assume the extended Limber approximation where the superscript I refers to the IA weighting function ( ) in Equation 10represents the linear growth factor,  cr denotes the critical density,  1 is a constant, and  pivot is an arbitrary redshift that has been set to 0.3.The IA contamination is parameterised by an amplitude  IA and a power law redshift dependence parameter  IA (see e.g.Bridle & King 2007).We remark that intrinsic alignment modelling and priors are also an open question (see Samuroff et al. 2021, for a recent study).We defer this issue to future work and focus on the impact of baryonic, massive neutrino and beyond-LCDM modelling.
We vary these parameters assuming a Gaussian prior N (0, 10 −4 ) for each of them, following Spurio Mancini et al. (2022).
We perform MCMC analyses on a set of mock cosmic shear data vectors, adopting the fiducial cosmology and baryonic boost parameters given in Table 1.We then consider two scenarios: In all analyses we vary all 5 cosmological parameters as well as To investigate the importance of the baryon parameters when looking to constrain cosmology and gravity, we perform analyses varying both the 7 × 2-and 3 × 2-parameter sets (see subsection 2.2).Note that all fiducial spectra are generated using the 7 × 2-parameter set.Lastly, to examine the impact that nonlinear scales have on the degeneracies, biases and constraints of the various parameters, we adopt three scale cuts, ℓ max = 1500 (pessimistic), ℓ max = 3000 (realistic) and ℓ max = 5000 (optimistic).The pessimistic case is well within the current capabilities of REACT and available fitting formulae such as HMCODE, as shown in Bose et al. (2020), which is what we employ here.The realistic case is also possibly within the current capabilities of available codes, for example using REACT to construct the boost with the newest version of HMCODE (Mead et al. 2020), and combining with state-of-the-art LCDM emulators such as the Euclid Emulator (Knabenhans et al. 2021).Finally, the optimistic case represents a possible modelling potential given, for example, high quality emulators for the pseudo power spectrum (Giblin et al. 2019).To implement these choices we use 50 ℓ bins in the range [30, 5000], which for our ℓ max = 1500, 3000 runs we cut to multipoles up to the corresponding values.For each setup we assume a simple Gaussian covariance matrix (Joachimi et al. 2007).Note that we leave for future work the inclusion of non Gaussian (Sato & Nishimichi 2013;Krause & Eifler 2017) and super sample (Takada & Hu 2013;Li et al. 2014;Lacasa et al. 2018) contributions (see also a discussion on their possible impact in section 4).
In summary, we study 3 scale cuts and 2 baryonic models for each of the 2 scenarios, giving a total 12 analyses.We vary 33 (25) parameters in total: 5 cosmological +   & log 10 |  R 0 |, 14 (6) baryonic, 2 intrinsic alignment and 10 redshift distribution shift parameters.We sample our posterior distribution using P C (Handley et al. 2015a,b), using 1000 live points to ensure accurate contours despite the large parameter space.

Results
The 2 constraints on all cosmological and baryonic parameters are summarised in Table 3 for all cases considered in this work.We also show the constraints and shifts of the recovered cosmological parameters in Figure 2. Note that log 10 |  R 0 | takes values below the effective LCDM-limit of −9 as we emulate down to −10.In section A we give the full set of posterior distributions for all parameters except the redshift-dependency parameters of the baryonic effects, which we find to be largely unconstrained.This is provided for both scenarios, all scale cuts and both choices of baryonic parameter sets.

Scenario A
In Figure 3 we show the marginalised 2-dimensional 95% (2) and 68% (1) distributions for the cosmological parameters as well as the beyond-LCDM parameters   and |  R 0 |.We summarise the main findings below and forward the reader to Table 3 for the exact constraints.
• In the conservative case (ℓ max = 1500) we find that cosmic shear is able to constrain all LCDM cosmological parameters at least at the ∼ 10% level (95 % C.L.), with almost an order of magnitude better constraints on Ω m and  8 . s is also better constrained, but to a lesser extent (3.8%).
• The inclusion of scales 1500 ≤ ℓ ≤ 5000 does not significantly improve the base LCDM cosmological parameters, with the most significant improvement being in constraints on ℎ (up to a 6% measurement) and Ω b (roughly a factor of 2 with respect to the conservative case).
• The inclusion of scales 1500 ≤ ℓ ≤ 5000 does not significantly improve the constraints on   and |  R 0 |.
• Reducing the number of baryonic degrees of freedom from 7×2 to 3 × 2 does not notably improve constraints on any parameter nor does it lead to any significant (up to 1) biases in the recovered cosmological parameters.We find at most ∼ 1 shifts on the   and |  R 0 | marginalised posteriors in the ℓ max = 5000 case.
These results show that the highly nonlinear scales have only a marginal amount of information on massive neutrinos and modified gravity which is not degenerate with baryonic degrees of freedom.Similarly, the baryonic degrees of freedom dampen the cosmological information one is able to extract from small, nonlinear scales.This was already shown in Copeland et al. (2018).There the authors illustrate how cosmological and dark energy constraints degrade with marginalisation of baryonic effects and that there is little improvement by going to smaller scales.Further, highly consistent results were also found in Bose et al. (2020), where similar poor gains in constraining power on log 10 |  R 0 | was found when including scales 1500 ≤ ℓ ≤ 3000.They do however obtain worse overall constraints compared to ours, by almost half an order of magnitude on  R 0 , despite their analysis not including baryonic degrees of freedom nor massive neutrinos.This is likely due to the difference in setup, as for example we use 10 tomographic bins while they only considered 3.
To investigate further the origin of this lack of constraining power on cosmological parameters at highly nonlinear scales, we plot the posteriors for the cosmological, intrinsic alignment and baryonic degrees of freedom.This is shown in Figure 4.It is clear from these plots that much of the information at nonlinear scales is associated with the baryonic degrees of freedom, with all of these parameters seeing significant gains in constraints when moving from ℓ max = 1500 to ℓ max = 5000.This is also clearly seen in Table 3. Notably, in the 7 × 2 parameter case, there are no large degeneracies between baryonic degrees of freedom and cosmological ones, which supports the lack of improvement in constraining power on cosmological parameters when moving to nonlinear scales.Baryonic parameters are in general less constrained than cosmological parameters, with the exception of log 10  c which can be constrained down to the ∼ 5% level even in conservative configurations.However, baryonic parameters are in general more heavily biased, up to 2 for  ej , when using the reduced model with 3 × 2 parameters over the one with 7 × 2 parameters.
In Table 3 we also show the marginalised constraints for Scenario A analyses but with all baryonic degrees of freedom fixed to their fiducial values.As in Schneider et al. (2020a), this represents an ideal situation where baryonic effects are completely known 4 .In Figure 5 we plot the posterior contours for the cosmological parameters in this setup where we assume perfect knowledge of the baryonic parameters, for the realistic configuration ℓ max = 3000, and we compare the results against those obtained with the same ℓ-cut, but varying either all of the 7 × 2 or all of the 3 × 2 baryonic parameters.We see that knowing the baryonic effects improves constraining power, relative to the 7 × 2 parameter model, on Ω  , ℎ and   by 50%, and by almost an order of magnitude on  R 0 .Unlike Schneider et al. (2020a) we see little improvement on   , which is likely due to the mild degeneracy of   with  R 0 and   .As expected, the constraint on Ω b is significantly improved, by roughly factor of 2, which just indicates that the baryon fraction is degenerate with baryonic effects. 4One can constrain these degrees of freedom from other observations such as gas X-ray emissions (Merloni et al. 2012;Schneider et al. 2020b).

Scenario B
In Figure 6 we show the marginalised 2-dimensional distributions for the cosmological parameters as well as the beyond-LCDM parameters.The main findings are consistent with Scenario A, with marginally worse constraints on  8 and   .This is driven by the strong degeneracies between  8 , |  R 0 | and   .At all scales, increasing log 10 |  R 0 | will naturally increase  8 as the modification works to enhance power, giving the correlation we see between these parameters.Conversely, massive neutrinos are known to suppress power and so we expect an anti-correlation between   and  8 .
We note that weak lensing alone is able to provide a clear detection of  () gravity, given a modification of log 10 |  R 0 | = 10 −6 , even in the presence of massive neutrinos.This detection is not changed notably by including scales 1500 ≤ ℓ ≤ 5000, and is moderately improved by reducing the baryonic parameters.These results are highly consistent with the results of Bose et al. (2020), which performed a similar analysis without massive neutrinos and no baryonic effects.
It is clear that even without baryons and massive neutrinos, most constraining power comes from scales ℓ ≤ 1500.
We also show the marginalised posterior distributions between the cosmological and baryonic parameters in Figure 7.This uncovers clear degeneracies between  () effects and baryonic physics in the reduced baryonic parameter model, which accounts for the gain in constraining power when moving to very nonlinear scales, which breaks these degeneracies.On the other hand, the full baryonic model does not exhibit such strong degeneracies and we gain less in this case when moving to highly nonlinear scales.
In Table 3 we report the expected constraints obtained in the realistic case ℓ max = 3000 assuming perfect knowledge of the baryonic parameters.In Figure 8 we compare these results against the contours for the 7 × 2 and 3 × 2 parameter models.From this comparison we observe that opening up the parameter space for baryonic parameters leads to significantly worse constraints on all cosmological parameters, particularly on |   0 |.

CONCLUSIONS
In this work we have presented a lightening fast pipeline to perform cosmic shear analyses in beyond-LCDM scenarios.To accomplish this we have created an emulator, REACTEMU-FR for the modification to the nonlinear LCDM power spectrum coming from Hu-Sawicki  () gravity and massive neutrinos based on halo model reaction predictions.Together with the state-of-the-art baryonic feedback emulator, BCEMU, the pipeline has sufficient accuracy to be applied to Stage IV surveys' mildly nonlinear data (roughly multipoles ℓ max ≤ 1500).Combining with available tools such as the Euclid Emulator (Knabenhans et al. 2021) would upgrade our pipeline to be safely applied to nonlinear scales, which we roughly define as multipoles ℓ ≤ 3000.Given the development of high quality emulators for the pseudo spectrum (Giblin et al. 2019), a necessary component of the  ()-massive neutrino boost, one could push to even larger multipoles, say ℓ ≤ 5000.We consider all these scale cuts in our paper and denote them as pessimistic, realistic and optimistic cases respectively.
Using our pipeline, we have run a number of mock data analyses using cosmic shear, providing forecasts for forthcoming surveys on modified gravity and massive neutrinos.Our main findings are as follows: (i) Scales between 1500 ≤ ℓ ≤ 5000 offer only mild improvements to constraining power on cosmology, with most power going into constraining non-degenerate baryonic parameters.
(ii) Using a reduced baryonic parameter set (3 × 2 = 6 instead of 7 × 2 = 14) does not impact constraints on any cosmological parameters significantly, but does improve massive neutrino constraints marginally.Further, using this reduced model may lead to significant biases in the baryonic degrees of freedom.
A number of papers have provided similar forecasts in the literature.First, we compare with the LCDM forecasts of Blanchard et al. (2020).We obtain surprisingly consistent constraints on Ω  ,  s and  8 , but significantly better constraints on Ω  and ℎ, of factors of 4 and 2 respectively.The differences we observe may be caused by our higher redshift range (they stop at  = 2.5) and our slightly higher sky fraction.Further, our nonlinear prescriptions are also very different.We make use of HMCODE while they employ , and we include baryonic feedback via BCEMU while they do not consider feedback at all.These differences to the absolute power can have significant impacts on the constraining power (Martinelli et al. 2021).
In any case, it is expected that the inclusion of cross-correlations with galaxy clustering, and the full 3×2pt analysis, will give far better constraints on all cosmological parameters (Blanchard et al. 2020).Future surveys will also offer constraints on baryonic physics (Huang et al. 2021) which can be combined with other probes of baryonic physics.As we have shown, this will lead to improved constraints on massive neutrinos and  () gravity, besides being interesting in its own right.However, we also note that our analysis has optimistically assumed a Gaussian covariance, ignoring super-sample contributions which in all likelihood will degrade the overall constraining power.A detailed investigation of this aspect is left for future work, as is the inclusion of optimal weighting schemes for physical scales (e.g.Bernardeau et al. 2014;Taylor et al. 2018a;Deshpande et al. 2020).
Moving on to studies of  () gravity, in Bose et al. (2020) they do not include baryonic or massive neutrino physics and have slightly different survey specifications, including far fewer tomographic bins (3 bins vs our 10 bins).Despite these differences, we find similar contributing constraining power when including multipoles between 1500 ≤ ℓ ≤ 3000 for both the LCDM fiducial and the  () fiducial mock data analyses as well as similar overall constraints on  ().
Similarly, Harnois-Déraps et al. ( 2022) do not include massive neutrinos nor baryonic effects, and have significant differences in the survey specifications and setup.Notably they use 5 tomographic bins up to  = 3, a quarter of our sky fraction, only consider autotomographic bins and have a restricted prior on  R 0 .Given this, they naturally report weaker overall constraints in their full MCMC analyses.Interestingly, their Fisher analysis finds that the constraints on  R 0 will improve when including highly nonlinear scales, in contrast to what we find.We remind the reader though that they do not consider baryons, which also dampen the  () constraints as well as their analysis having fixed the other cosmological parameters, and so representing an absolutely ideal scenario.More similar to our analysis is that of Schneider et al. (2020a,b) which include a similar baryonic model as used in this work, but without redshift evolution, and including massive neutrinos and  () gravity, albeit separately and using different theoretical prescriptions.They also use a different nonlinear power spectrum prescription, as well as only include 3 tomographic bins in the range 0.1 ≤  ≤ 1.5, and take ℓ max = 4000.The authors report significantly weaker constraints on cosmological,  () and neutrino mass parameters.The number of tomographic bins may lead to our improved constraints, particularly on cosmological parameters.Schneider et al. (2020b) do find almost 2 orders of magnitude worse constraints on  R 0 than Bose et al. (2020) but include baryonic degrees of freedom which significantly worsen the constraints.On the other hand, we also include baryonic degrees of freedom and obtain ∼2.5 orders of magnitude better constraints.
Another possibility could be the difference in tomographic binning and maximum redshift, which would contribute highly to these improved constraints as  R 0 ,  8 and   all have different redshift evolutions.The power spectrum becomes highly insensitive to  R 0 at early times (it being a late-time modification), breaking degeneracies with  8 and   .This means more redshift information will improve constraints, particularly high redshift information which was also pointed out in Harnois-Déraps et al. (2022), where they showed that most Fisher information on modified gravity is contained in the highest redshift bin.Finally, it was shown in Taylor et al. (2018b) that a large number of bins are needed to capture the high-redshift information when using equipopulated binning.
Lastly, a recent analysis (Aricò et al. 2023) of the Dark Energy Survey year 3 data including baryons, massive neutrinos and the same intrinsic alignment model employed here has shown consistent constraints on the  8 parameter after considering they employ a sky fraction that is roughly 4 times less than that employed here.Their constraints on Ω  are still significantly worse than ours, and the forecast of Blanchard et al. (2020) for example.We also employ 6 times the number density of their analysis, which will vastly improve constraining power on small scales.Survey specifications and the full power of Stage IV make such comparisons difficult, but we deem these results to be well within our expectations.It is also worth noting that this analysis (Aricò et al. 2023) was possible due to improvements in the prior ranges of the BACCO emulator.Priors can be informative and it is important to have accurate nonlinear modelling across a wide parameter space (see Carrilho et al. 2023, for example), making REACT emulation an easy and viable means of performing data analyses for beyond-LCDM scenarios.
We conclude by noting that such a high dimensional parameter space analysis is both necessary and challenging when looking to optimise reliable information extraction from Stage IV surveys.It is also worthwhile as nonlinear scales can provide considerable information on interesting extensions to LCDM.Our analyses has focused on a minimal extension to LCDM, and serves as a 'proof of principle' and 'proof of worth' for the tools used in this paper, notably COSMOPOWER and REACT, which we combined to provide REACTEMU-FR.One key goal of forthcoming Stage IV surveys is to provide clues for solving the looming theoretical issues of Dark Energy and, even more broadly, the cosmological constant problem (see Bernardo et al. (2023) for a recent review on modified gravity solutions).Thus, pipelines able to accurately model comprehensive extensions to LCDM becomes an absolute necessity.Such extensions have already been proposed (see Bose et al. 2022, for example), but include a number of new degrees of freedom, making the need for emulation even more clear, as is the development of fast inference pipelines based on gradient-based sampling techniques (see e.g.Ruiz-Zapatero et al. 2023;Campagne et al. 2023).It is straightforward to include existing models of gravity and dark energy into this emulation framework.It is also the goal of ongoing work to integrate these and perform similar tests as part of necessary preparation for the data releases of the largest galaxy surveys to date.model, both when considering a fiducial LCDM data vector and  () data vector (Scenarios A & B respectively in the main text).All scale cuts are also shown.Similarly, in Figure A2 we show the results for the reduced 3 × 2 baryonic parameter model.
These plots visually highlight the vast parameter space that one needs to probe in order to comprehensively account for all known physics, and probe minimal extensions to LCDM.This space can forseeably increase when looking to move deeper into the nonlinear regime or when wishing to probe more general deviations to LCDM.

Figure 2 .
Figure 2. The marginalised means and 2 constraints on the cosmological and beyond-LCDM parameters, normalised to their fiducial values for Scenario A (left) and Scenario B (right), and for all scale cuts and baryonic parameter models.

Figure 6 .
Figure 6.Same as Figure 3 but for Scenario B.

Figure 7 .
Figure 7. Same as Figure 4 but for Scenario B.

Figure A1 .
Figure A1.Full set of marginalised 2-dimensional posterior distributions for the 7 × 2 baryonic parameter model.We show all scale cuts as well as the case where we fit to a LCDM fiducial data vector (bottom left) and  () fiducial data vector (top right).

Figure A2 .
Figure A2.Same as Figure A1 but for the 3 × 2 baryonic parameter model case.

Table 1 .
Varied cosmological, beyond-LCDM and baryonic parameters at  = 0, with prior ranges and fiducial values.The cosmological parameter fiducials are taken to be the Planck 2018 best-fits, while the baryonic parameter fiducials are the default values of BCEMU.

Table 2 .
Priors and fiducial values for time dependent baryonic parameters.

Table 3 .
Percent 2 constraints on cosmological and baryonic parameters (at  = 0) for both Scenarios, both baryonic parameter models and all scale cuts.For Scenario A, we quote the upper 2 bound for log 10 |  R 0 | as the fiducial is  R 0 = 0.