© The Authors 2026

1. Introduction

The European Space Agency mission Euclid (Euclid Collaboration: Mellier et al. 2025) is set to revolutionise our understanding of the large-scale structure (LSS) of the Universe by providing high-precision measurements of cosmic shear, galaxy clustering, and the distribution of matter on cosmological scales (e.g. Heymans et al. 2021; DES Collaboration 2022; Miyatake et al. 2023). These data will enable tight constraints on key cosmological parameters, such as the equation of state of dark energy, the sum of neutrino masses, and the amplitude of primordial density fluctuations (Euclid Collaboration: Mellier et al. 2025).

The classical cosmic shear analysis involves measuring the correlations between galaxy ellipticities as a function of their separation, quantified by the two-point correlation function of the shear (γ-2PCF, e.g. Asgari et al. 2020; Amon et al. 2022; Dalal et al. 2023). This statistic benefits from a well-established theoretical framework, primarily capturing the multiscale variance of the lensing field. However, gravitational collapse introduces non-Gaussian features into the shear field, which are not fully captured by γ-2PCF or other second-order statistics.

As a result, two-point statistics cannot extract all cosmological information, leading to parameter degeneracies, such as the one between the matter density parameter Ω_m and normalisation of the power spectrum σ₈, which limit constraints to combinations such as $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$. Various higher-order statistics (HOS) have been proposed in the literature to access additional information (see e.g. Euclid Collaboration: Ajani et al. 2023).

A compelling HOS for probing non-Gaussian features of the LSS is the matter bispectrum (e.g. Bernardeau et al. 2003; Takada & Jain 2003), which captures three-point correlations in the density field and provides complementary information to the widely studied two-point statistics, such as the power spectrum (Halder et al. 2021; Halder & Barreira 2022; Heydenreich et al. 2023; Burger et al. 2024; Gomes et al. 2025).

One of the biggest challenges of HOS is the precise modelling of intermediate to small scales, as HOS extract most additional information on nonlinear scales (Foreman et al. 2020). However, these small scales are heavily affected by baryonic feedback processes, such as gas cooling, star formation, and feedback from active galactic nuclei (AGN), which modify the distribution of matter (van Daalen et al. 2011; Chisari et al. 2018). These processes significantly affect the matter distribution, especially in the nonlinear regime, and must be accounted for in cosmological analyses to avoid biased constraints on fundamental physics (Semboloni et al. 2011; Schneider et al. 2019).

Semboloni et al. (2013) used a halo-model approach to predict the effect of baryons on second-order shear statistics and a perturbation theory approach to predict the effect of baryons on the three-point shear statistics. By applying their method to the OverWhelmingly Large Simulations (Schaye et al. 2010), they found that, similar to second-order statistics, the effect of baryonic processes on third-order statistics is not negligible. Moreover, the two statistics have different sensitivities and dependences on the baryonic parameters. Therefore, the combination of second- and third-order statistics has the potential to self-calibrate baryonic feedback parameters, thus improving constraining power on cosmological parameters and reducing projection effects that arise when marginalising high-dimensional posterior distributions over subsets of parameters, which lead to shifts in one- or two-dimensional marginal posteriors that do not reflect the true structure of the full joint distribution.

Hydrodynamical simulations (e.g. Le Brun et al. 2014; Nelson et al. 2015; Schaye et al. 2015, 2023; McCarthy et al. 2018; Pakmor et al. 2024) offer the most detailed insight into baryonic processes and their effects on the power spectrum (Schaller et al. 2025). However, their computational expense makes exploring the high-dimensional parameter space necessary for future cosmological surveys impractical. A promising method for incorporating baryonic effects involves modifying Gravity-Only (GrO) simulations to mimic the influence of baryons using physically motivated models. These approaches, often referred to as Baryon Correction Models (BCM) or baryonification techniques (Schneider & Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020), are computationally efficient and enable the generation of extensive training sets. Such datasets are instrumental in building fast emulators that predict baryonic suppression as a function of cosmological parameters and a few physically motivated baryonic parameters (Aricò et al. 2021a; Giri & Schneider 2021). Using these tools, the baryonic suppression effects observed in various hydrodynamical simulations can be replicated with percent-level accuracy down to scales of k = 5 − 10 h Mpc⁻¹. Although much progress has been made in emulating baryonic effects on the matter power spectrum, the bispectrum or higher-order statistics have remained comparatively underexplored (see e.g. Aricò et al. 2021b; Anbajagane et al. 2024; Zhou et al. 2025). Given the additional cosmological information contained in the bispectrum and its sensitivity to non-Gaussianities, we develop here, as an essential extension, an emulator that can accurately capture baryonic effects on the bispectrum and is fast enough to be used in Bayesian pipelines.

We do this by expanding the work of Aricò et al. (2021a) on the matter power spectrum and present the first emulator predicting the impact of baryonic processes on the matter bispectrum. We model the relevant baryonic processes employing high-resolution simulations from the BACCO suite (Angulo et al. 2021), post-processed using the baryonification technique developed in Aricò et al. (2020) and the cosmology scaling of Angulo & White (2010). We train deep neural networks to predict the power spectrum and bispectrum of our simulations. Although an emulator of baryonic effects on the power spectrum was already presented in Aricò et al. (2021a), we expand their work to have consistent modelling for the power and bispectrum with a broader parameter space and a scale that extends to larger k-values. We validate these emulators against the FLAMINGO set of high-resolution hydrodynamical simulations (Schaye et al. 2023; Kugel et al. 2023), each incorporating different physical prescriptions for feedback mechanisms, star formation, and gas dynamics. Training on various scales and redshifts ensures that our emulators cover the relevant parameter space for the Euclid mission.

Specifically, we show that our BCM emulator can accurately reproduce the FLAMINGO baryonic effects on the second- and third-order moments of the aperture mass lensing statistic (Jarvis et al. 2004; Schneider et al. 2005) across different configurations of box sizes, particle masses, and strengths of AGN feedback. Furthermore, we combine second- and third-order weak lensing statistics to forecast cosmological constraints under the impact of baryonic effects for a non-tomographic analysis of the Euclid first data release (DR1), quantifying the importance of accounting for baryonic effects in the power spectrum and bispectrum measurements. Our results demonstrate that ignoring baryonic effects could lead to significant biases in interpreting Euclid data or require scale cuts, leading to a severe reduction in constraining power, highlighting the necessity of incorporating these effects into cosmological analyses.

This paper is organised as follows. In Sect. 2, we describe the suites of simulations used for training, validating, and building the emulator. We also describe the GrO simulations used for estimating the Euclid-DR1 covariance matrix for the weak lensing statistics. Section 3 presents the estimation of the bispectrum and its validation. In Sect. 4, we present the training of the BCM emulator of the power spectrum and bispectrum. In Sect. 5, we present the methodology for converting the spectra to weak lensing statistics, which we then use in Sect. 6 to perform the non-tomographic Euclid-DR1 forecast. Finally, we discuss the implications of our findings and potential future extensions in Sect. 7.

2. Simulations

2.1. BACCO simulations

This section describes the BACCO N-body simulations used in this study to build a BCM and quantify the baryonic effects on the matter power and bispectrum. The BACCO simulation suite is specifically designed to provide highly accurate predictions of nonlinear cosmic structure formation as a function of cosmology. The gravitational evolution is computed using L-Gadget3 (Springel 2005; Angulo et al. 2012, 2021), a variant of the Gadget code, with a Plummer-equivalent softening length of ϵ = 5 h⁻¹ kpc. The numerical parameters of the suite, including force and mass resolution, were chosen to achieve convergence at the 1% level in the nonlinear power spectrum at z = 0 and k ∼ 10 h Mpc⁻¹. This setup has been validated against other N-body codes, performing exceptionally well in the Euclid Code Comparison project (Schneider et al. 2016), where it agreed within 2% with most codes up to k ∼ 10 h Mpc⁻¹ (Angulo et al. 2021). Each simulation explores a distinct cosmology, with parameters carefully selected to maximise the coverage of the parameter space when combined with rescaling algorithms (Contreras et al. 2020; Angulo et al. 2021).

In this work, we employ a single gravity-only simulation, a higher-resolution version of TheOne simulation presented in Aricò et al. (2021a). This simulation follows the evolution of 2288³ particles of mass m_p ≈ 9.55 × 10⁸ h⁻¹ M_⊙ within a periodic box of L = 512 h⁻¹ Mpc on a side. We report its cosmology in Table 1. The initial conditions were generated at z = 49 using second-order Lagrangian perturbation theory (2LPT). To minimise cosmic variance, the amplitudes of Fourier modes were fixed to the ensemble average of the linear power spectrum, and we average the summary statistics of two realisations having opposite phases, following the fix and pair approach described by Angulo & Pontzen (2016). The particle catalogue was downsampled by a factor of 4³ for computational and memory efficiency. We refer to Aricò et al. (2021b) for detailed convergence tests of resolution and cosmic variance, which show that baryonic effects on the power spectrum and bispectrum produced with the BCM on top of this simulation are converged at the percent level.

2.1.1. Cosmology rescaling

We explore the cosmological space with the cosmology rescaling algorithm originally proposed in Angulo & White (2010), and further refined by Zennaro et al. (2019) and Contreras et al. (2020). This method computes and applies optimal space and time scale factors, which minimise the linear mass variance between two different cosmologies. Furthermore, after applying these scale factors, one can also correct the large-scale flow motions and the small-scale internal structure of haloes. This algorithm can accurately reproduce the matter and halo clustering in a wide cosmology parameter space, including massive neutrinos and dynamical dark energy, with just a few simulations and in only a few seconds. In this work, we only vary the cosmological parameters Ω_m, Ω_b, σ₈, that have been shown in several works to be the dominant cosmological dependence of baryonic effects, at the level of the matter power spectrum (Schneider et al. 2019; Aricò et al. 2021a). We tested in Appendix C with the BACCO simulation suites that this statement is also valid for the matter bispectrum.

2.1.2. Baryonic correction model

To incorporate the effects of baryonic physics into the mass distribution in GrO simulations, we use the ‘baryonification’ framework (Schneider & Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020). This method applies physically motivated prescriptions for processes such as star formation, gas cooling, and AGN feedback to perturb the positions of particles in gravity-only N-body simulations, thereby modifying the mass distribution.

One key advantage of the baryonification approach is its flexibility in exploring possible modifications to the matter power and bispectrum, accounting for the uncertainties inherent in baryonic physics. Model parameters can be constrained using observational data and tested with hydrodynamical simulations. In this work, we adopt the implementation detailed in Aricò et al. (2021b), which has been successfully tested along with the cosmology rescaling algorithms in simultaneously reproducing the power spectra and bispectra of different hydrodynamical simulations. The model considers several baryonic components. The first component is the bound gas in haloes, whose density is described by a double power law with the transition radius and slopes being free parameters. In this work, we adopt only one free parameter for the bound gas shape, θ_inn, which represents the inner characteristic scale where the gas slope varies. The outer characteristic scale, θ_out, was fixed to one, which we verify in Appendix C. Both characteristic scales are reported in units of the halo radius, r₂₀₀, which in turn is defined as the radius where the halo’s density exceeds the Universe’s critical density by a factor of 200. The bound gas mass fraction, f_BG, is parametrised as

$$\begin{array}{c}\hfill f_{\mathrm{BG}}(M_{200})=\frac{f_{\mathrm{b}}-f_{\mathrm{CG}}-f_{\mathrm{SG}}}{1-{(M_{\mathrm{c}}/M_{200})}^{\beta }},\end{array}$$(1)

where f_b = Ω_b/Ω_m is the cosmic baryon fraction, M₂₀₀ is the mass inside radius r₂₀₀, M_c and β are free parameters, and f_CG, SG are the central and satellite galaxy stellar mass fractions, respectively. An abundance-matching parametrisation gives the central and satellite galaxy fractions (Behroozi et al. 2013), which have the same parametric form in halo mass and redshifts and their parameters are assumed to be linearly dependent.

In this work, we follow Aricò et al. (2021a) and free only one parameter for the stellar component, M_{1, z0, cen}. This parameter represents the characteristic halo mass for which the central galaxy to halo mass fraction is set to ϵ_0, cen = 0.023 at z = 0.

A fraction of the gas is assumed to be ejected by feedback, with a constant density and an exponential cutoff at a given scale, parametrised by a free parameter η. The sum of stellar and gas mass fractions, including the ejected gas, is by construction equivalent to the cosmic baryon fraction f_b. We report the priors on the varied BCM parameters in Table 2. Moreover, we fix to a fiducial value the BCM parameters $M_{\mathrm{inn}}=2.3\times {10}^{13}{h}^{-1}{\mathrm{M}}_{\odot }$, which are verified in Appendix C. Finally, we vary the scale factor a = 1/(1 + z) to account for the redshift evolution of baryonic effects. We note that in the following, we write M₁ instead of M_{1, z0, cen}.

2.2. FLAMINGO simulations

To test our baryon model, we use the FLAMINGO simulations¹, in particular, we use the convergence maps developed and described in Broxterman et al. (2024). A detailed description of the FLAMINGO simulations, their performance, and calibration strategy can be found in Schaye et al. (2023), Kugel et al. (2023), and McCarthy et al. (2023).

The simulations were conducted using the SWIFT hydrodynamic code (Schaller et al. 2024) with the SPHENIX smoothed particle hydrodynamics (SPH) implementation (Borrow et al. 2022). Neutrinos were modelled as massive particles using the method developed by Elbers et al. (2021), which minimises particle shot noise. The simulations incorporate radiative cooling and heating on an element-by-element basis (Ploeckinger & Schaye 2020), star formation following the prescription of Schaye & Dalla Vecchia (2008), and time-dependent stellar mass loss as outlined by Wiersma et al. (2009). Supernova and stellar feedback are implemented kinetically, with neighbouring particles receiving kicks that conserve energy, linear momentum, and angular momentum, as described in Chaikin et al. (2023). Gas accretion onto supermassive black holes and subsequent thermal AGN feedback are modelled based on Booth & Schaye (2009), while kinetic jet feedback follows the AGN jet implementation of Huško et al. (2022), where gas particles are accelerated to a fixed jet velocity aligned with the black hole’s spin.

In this work, we use the box with 2.8 comoving Gpc side length, 5040³ dark matter (DM) and baryon particles, and 2800³ massive neutrino particles, where the baryonic particle mass is 1.07 × 10⁹ M_⊙ and the DM particle mass is 5.65 × 10⁹ M_⊙. We use all eight available light cones with different observers. The cosmology of these runs was taken from the Dark Energy Survey year three “3 × 2pt+ All Ext”. ΛCDM cosmology (DES Collaboration 2022), which we denote as ‘D3A’. Complementary to this run, we use all nine baryon descriptions available for the D3A cosmology given in Table 1. These additional runs have the same resolution, but with a box size of 1 Gpc, 1800³ baryonic particles, and 1000³ massive neutrino particles. We also use the higher and lower mass resolution for the fiducial baryon feedback description (L1_m8 and L1_m10). In the following, we label the gravity-only simulations as ‘GrO’ and those including baryonic feedback as ‘Hydro’.

At the D3A cosmology, eight runs were calibrated to different galaxy stellar mass functions (M_★) and/or gas fractions in clusters (f_gas), or they were varied in their overall AGN subgrid feedback prescriptions. The values of the subgrid parameters were adjusted so that the stellar mass function and/or the gas fractions in clusters deviated by specific standard deviations (ΔM_★ and Δf_gas) from the fiducial model, as shown in the second and third columns of Table 3.

Four of these eight runs focused exclusively on varying gas fractions in clusters. These are identified as fgas± nσ, where n represents the number of standard deviations by which the gas fraction is modified. One additional run was calibrated to match a stellar mass function shifted to a lower mass by 1σ (M^*–σ). Another run, labelled M^*–σ_fgas–4σ, incorporated variations in both observables. The fiducial AGN feedback model uses a thermal implementation, but two alternative models employ a kinetic jet feedback mechanism. These are denoted as Jet and Jet_fgas–4σ, with the latter targeting reduced gas fractions in clusters. The jet models provide a way to evaluate the sensitivity of observables to variations in subgrid models calibrated to the same data. Lastly, we also make use of a simulation at a lower σ₈ than the D3A, which is reported in Table 1, denoted as LS8, and uses the baryon description of the fiducial case L1_m9.

2.2.1. Weak lensing convergence maps

The weak lensing convergence maps κ(ϑ, χ) were created using ray-tracing for each comoving slice, where ϑ refers to the sky coordinates and χ to the comoving distance (Broxterman et al. 2024). A key quantity of gravitational lensing is the magnification or Jacobian matrix defined as

$$\begin{array}{c}\hfill A_{\mathit{ij}}(\mathit{\vartheta },\chi )=\frac{\partial {\beta }_i}{\partial {\vartheta }_j}={\delta }_{\mathrm{K},ij}-\frac{{\partial }^2\varphi }{\partial {\theta }_i\partial {\theta }_j}\approx \left(\begin{array}{cc}1-{\varphi }_{,11}& -{\varphi }_{,12}\\ -{\varphi }_{,12}& 1-{\varphi }_{,22}\end{array}\right),\end{array}$$(2)

where δ_K is the Kronecker delta symbol, _,ij = ∂/∂θ_i ∂/∂θ_j, β denotes the original position of the ray-traced light ray, and can be computed as an integral along the comoving distance over the gradient of the deflection potential ϕ (see Bartelmann & Schneider 2001). Given the fact that ∇²ϕ = 2κ, we can also express the magnification as

$$\begin{array}{c}\hfill A(\mathit{\vartheta },\chi )=\left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2\\ -{\gamma }_2& 1-\kappa +{\gamma }_1\end{array}\right),\end{array}$$(3)

where

$$\begin{array}{c}\hfill {\gamma }_1(\mathit{\vartheta },\chi )=\frac{1}{2}[{\varphi }_{,11}-{\varphi }_{,22}](\mathit{\vartheta },\chi )\&{\gamma }_2(\mathit{\vartheta },\chi )={\varphi }_{,12}(\mathit{\vartheta },\chi ).\end{array}$$(4)

Based on the magnification matrix derived from density maps in HEALPix format (Górski et al. 2005) with an N_side = 8192, which corresponds to a pixel size of $0\stackrel{\prime }{.}43$, the shear signal κ(ϑ, χ) and γ(ϑ, χ) are computed at each pixel location. Next, these convergence maps are combined using the Limber integration (Limber 1954)

$$\begin{array}{c}\hfill \kappa (\mathit{\vartheta })={\displaystyle {\int }_0^{{\chi }_{\max }}}\mathrm{d}\chi n[z(\chi )]\frac{\mathrm{d}z}{\mathrm{d}\chi }\kappa (\mathit{\vartheta },\chi ),\end{array}$$(5)

where n(z) is the redshift probability distribution of source galaxies given in Euclid Collaboration: Blanchard et al. (2020) as

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

with $z_0=0.9/\sqrt{2}$. This redshift distribution is cut at z_max = 3, so that χ_max = χ(z_max) and it is normalised.

2.3. Takahashi simulations

As in Burger et al. (2023), we use the Takahashi et al. (2017) simulations to estimate a Euclid-like covariance matrix for our lensing statistics. The simulations by Takahashi et al. (2017) track the nonlinear evolution of 2048³ particles within a series of nested cosmological volumes. These volumes start with a side length of L = 450 h⁻¹ Mpc at low redshifts and grow progressively larger at higher redshifts, resulting in 108 full-sky realizations. These simulations were generated using the Gadget-3 N-body code (Springel 2005) and are publicly available². The cosmological parameters are set to values consistent with WMAP 9 (Hinshaw et al. 2013) and reported in Table 1. Lensing information from these simulations is provided using the Born approximation as 108 full-sky realisations of the convergence and shear, each subdivided into 38 redshift slices in ascending order. To construct a non-tomographic Euclid-DR1-like covariance matrix, we built the weighted average of the first 30 κ maps with ascending redshift. The weight for each redshift slice is measured by integrating the n(z) given in Eq. (6) over the corresponding width of the redshift slice. To have a footprint of roughly 2000 deg², we divided the 30 full-sky maps into 19 patches, each with a size of roughly 2146 deg². Shape noise is added to the convergence maps by drawing random numbers from a Gaussian distribution with a vanishing mean and a standard deviation, as follows:

$$\begin{array}{c}\hfill \sigma =\frac{{\sigma }_{ϵ}}{\sqrt{n_{\mathrm{eff}}{A}_{\mathrm{pix}}}},\end{array}$$(7)

with the pixel area as A_pix = 0.18 arcmin² (N_side = 8192), n_eff = 30 arcmin⁻² as the effective galaxy number density, and shape noise contribution σ_ϵ = 0.3 as stated in Euclid Collaboration: Blanchard et al. (2020).

3. Estimation of the bispectrum

We measure the bispectra and power spectra employing the estimator by Scoccimarro (2000), which defines the bispectrum at k-bins k = {k₁, k₂, k₃} for a simulation of volume V as

$$\begin{array}{c}\hfill B(k_1,k_2,k_3)=\frac{{\prod }_{i=1}^3{\displaystyle {\int }_{\mathrm{\Delta }{k}_i}}{\mathrm{d}}^3{q}_i\widehat{\delta }({\mathit{q}}_i){\delta }_{\mathrm{D}}({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)}{{\displaystyle {\prod }_{i=1}^3{\int }_{\mathrm{\Delta }{k}_i}}{\mathrm{d}}^3{q}_i{\delta }_{\mathrm{D}}({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)},\end{array}$$(8)

where $\widehat{\delta }$ is the Fourier transform of δ, δ_D is the Dirac delta distribution, and Δk_i are the bins of wavenumber k_i. The power spectrum is defined by

$$\begin{array}{c}\hfill P(k)=\frac{{\int }_{\mathrm{\Delta }k}{\mathrm{d}}^{3}q{|\widehat{\delta }(\mathit{q})|}^2}{{\int }_{\mathrm{\Delta }k}{\mathrm{d}}^{3}q}·\end{array}$$(9)

We further use

$$\begin{array}{c}\hfill {\delta }_{\mathrm{D}}({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)={\displaystyle {\int }_{L_{\mathrm{box}}}}\frac{{\mathrm{d}}^{3}x}{{(2\pi )}^3}{\mathrm{e}}^{-\mathrm{i}\mathit{x}·({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)},\end{array}$$(10)

so the q-integrals become separable, and we obtain

$$\begin{array}{c}\hfill B(k_1,k_2,k_3)=\frac{{\int }_{L_{\mathrm{box}}}\frac{{\mathrm{d}}^{3}x}{{(2\pi )}^3}{\displaystyle {\prod }_{i=1}^3}{I}_{k_i}^{\delta }(\mathit{x})}{{\displaystyle {\int }_{L_{\mathrm{box}}}}\frac{{\mathrm{d}}^{3}x}{{(2\pi )}^3}{\displaystyle {\prod }_{i=1}^3}{I}_{k_i}^{\mathbf{1}}(\mathit{x})},\end{array}$$(11)

where L_box is the length of the simulation box.

$$\begin{array}{c}\hfill I_k^{\delta }(\mathit{x})={\displaystyle {\int }_{\mathrm{\Delta }k}}{\mathrm{d}}^{3}q\widehat{\delta }(\mathit{q}){\mathrm{e}}^{-\mathrm{i}\mathit{q}·\mathit{x}}\mathrm{and}{I}_k^{\mathbf{1}}(\mathit{x})={\displaystyle {\int }_{\mathrm{\Delta }k}}{\mathrm{d}}^{3}q{\mathrm{e}}^{-\mathrm{i}\mathit{q}·\mathit{x}}.\end{array}$$(12)

The I_k^δ(x) can be quickly obtained using Fast Fourier Transforms (FFT). To make feasible the estimation of thousands of spectra using the large grids required for this project, we produced a new implementation of these FFTs that runs on GPUs using the JAX framework, which provides a speed-up of a factor up to 10 with respect to CPUs codes. We publicly released the code BiG (Linke 2024), and show in Appendix A how its estimates for the bispectrum agree with those of CPUs codes such as bskit (Foreman et al. 2020) to better than 0.01%.

This estimator averages over the number of triangle configurations inside a (k₁, k₂, k₃) bin. However, as we are interested in a bispectrum prediction for a specific configuration, we need to define which ‘effective’ triangle configuration any measured bispectrum value corresponds to. Following Oddo et al. (2020), we assign each bin an effective triangle configuration (k_1, eff, k_2, eff, k_3, eff) as

$$\begin{array}{c}\hfill k_{i,\mathrm{eff}}(k_1,k_2,k_3)=\frac{{\prod }_{j=1}^3{\displaystyle {\int }_{\mathrm{\Delta }{k}_j}}{\mathrm{d}}^3{q}_j{q}_i{\delta }_{\mathrm{D}}({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)}{{\displaystyle {\prod }_{i=1}^3{\int }_{\mathrm{\Delta }{k}_i}}{\mathrm{d}}^3{q}_i{\delta }_{\mathrm{D}}({\mathit{q}}_1+{\mathit{q}}_2+{\mathit{q}}_3)}·\end{array}$$(13)

We use these effective ks when creating the emulator in the following sections and, for simplicity, drop the ‘eff’ subscript.

For a given density field with N³ grid points and box length L_box, the theoretically accessible k modes are k ∈ [2k_ny, Nk_ny], where k_ny = π/L_box is the Nyquist frequency. In our case, we use N = 858, which leads to k_max = 5.3 h Mpc⁻¹.

To access smaller scales (larger k) without increasing the number of grid points, we use the folding technique presented in Jenkins et al. (1998) and Colombi et al. (2009) and validated for the bispectrum by Aricò et al. (2021b). The particle distribution is ‘folded’ by wrapping particles into a sub-box of side length L/f. Given the same grid and the k_ny defined in the full box, the new Nyquist frequency is k_ny′=f k_ny. In this work, we deploy three foldings f = [1, 2, 4], such that the lowest k_min is 0.0123 h Mpc⁻¹ and the largest k_max is 21.058 h Mpc⁻¹. While combining the three foldings for the power spectrum is straightforward, it is more complicated for the bispectrum. The reason is the missing squeezed triangles, where one of the k_i values can only be measured in one folding, while the other k_i values are only measurable in the other. In the following, we will allow the emulator to extrapolate to these missing squeezed triangles and validate this on the lensing analysis below, since we currently cannot measure bispectra where all k-configurations fit in a single grid. We also notice that we use the symmetries of the bispectrum B(k₁, k₂, k₃) = B(k₂, k₃, k₁) = B(k₃, k₂, k₁) = B(k₃, k₁, k₂) = B(k₁, k₃, k₂) = B(k₂, k₁, k₃) and that the three sides need to build a triangle, k₃ ≤ k₁ + k₂, which reduces the required measurements. We average the measurements from two opposite-phase realizations to suppress cosmic variance for both power spectra and bispectra. We also subtract the shot noise contribution, which is given as ${\overline{n}}^{-1}=L_{\mathrm{box}}^3/N_{\mathrm{par}}$ for the power spectrum and $B_{\mathrm{sn}}(k_1,k_2,k_3)={\overline{n}}^{-2}+{\overline{n}}^{-1}[P(k_1)+P(k_2)+P(k_3)]$ for the bispectrum (Peebles 1980). Lastly, we note that we do not deconvolve the window function introduced by the mass assignment scheme, but we checked that the effect is smaller than 0.02% when taking ratios of power spectra and bispectra as in our case.

4. Training and validation of the emulator

After measuring the power and bispectra for the GrO and baryonified cases for 1200 nodes distributed in a Latin Hypercube defined by the prior mentioned in Sect. 2, we are ready to build the BCM of the matter power spectrum, S(k) = P^Hydro/P^GrO, and bispectra, R(k) = B^Hydro/B^GrO. For our emulation, we use the CosmoPower neural network (NN) library described in Spurio Mancini et al. (2022). In principle, we could use only the six baryon parameters as input and then interpolate between different k-values. However, we found that it is more convenient, accurate and faster to use also the k-value as input parameters, as this allows us to predict the power and bispectra for any k-values even if they were not measured in the first place. Moreover, in this way, we can clean our training set more efficiently, removing the single B(k₁, k₂, k₃) configurations where B(k₁, k₂, k₃) < 2 B_sn(k₁, k₂, k₃) without having to exclude the whole sampling node. With this procedure, we discard only about 2.5% of the total measurements. We note that we emulate the three-dimensional power spectrum and bispectrum, as opposed to the projected two-dimensional lensing statistics presented in the next sections, to have the flexibility to employ the emulators with different survey setups and projections, without the need to perform the training of more emulators. To further speed up the cosmological analyses, specific emulators for a given survey configuration can be easily developed starting from the ones presented here.

We randomly split our samples of 1200 nodes into 200 points used for testing and 1000 for training. After some iterations, we found that an architecture with four hidden layers and 32 neurons, the Adam optimizer, a mean squared difference loss function, a learning rate that starts at 0.01 and decreases by a factor of 10 if the loss does not improve for 100 epochs of training, and a sigmoid activation, gives the best results. The reason for using as few layers and neurons as possible is to prevent overfitting and attempt to improve the extrapolation to k-values not present in the training. These are, for instance, the squeezed triangles that do not lay in any of the three foldings used and k-values smaller than k_min = 0.0123 h Mpc⁻¹. We note that, while for the power spectrum, the BCM to large scales, k < 0.0123 h Mpc⁻¹, can straightforwardly be set to one, this is not obvious for the bispectrum given its sensitivity to different scales (see Fig. C.2). Therefore, a more complex extrapolation is needed. We implicitly test the accuracy of our NN extrapolation by performing mock analyses using the expected characteristics of Euclid DR1 in Sect. 6, leaving a more thorough explicit test for future works.

We first show the accuracy of our BCM power spectrum emulator in Fig. 1. We find an accuracy on S(k), which we define as ΔS(k) = S^emu/S^true − 1, smaller than 0.5% for the 68% percentile, and 3% for the 95% percentile. In other words, this means that 68% of the predictions are better than 0.5% for the power spectrum. We show the accuracy of the bispectrum emulator, defined as ΔR(k) = R^emu/R^true − 1, in Fig. 2, as a scatter plot where the corresponding 68% and 95% percentiles of ΔR(k) are shown in the upper two and lower two rows, respectively. We used linear interpolation between the measured k-values for a better illustration. The linear extrapolation into regions without measurements, such as the upper left corner for k₃ = 5 h Mpc⁻¹, should be taken cautiously. Furthermore, we show in Figs. B.1 and B.2 two different illustrations of the accuracy of the bispectrum emulator, where it is seen that the accuracy is smaller than 2% for the 68% percentile and similar across the different foldings and not particularly problematic at the transition from one folding to the next. Since we took 200 nodes for testing, which are also distributed in a Latin hypercube, the accuracies shown here marginalize over all parameters. However, we found only a minimal dependency of the accuracy towards larger M_c. All the other parameters showed no dependency at all. Overall, we conclude that the emulator is less accurate for very squeezed triangles, especially where the measurement is at the edge of the k-range.

Fig. 1. Accuracy of the emulator for baryonic effects on the matter power spectrum. The y-axis is defined as ΔS(k) = Semu/Strue − 1, where S(k) = PBCM(k)/PGrO(k). The horizontal dashed lines show the 0.5% region.

Fig. 2. Accuracy of the emulator for baryonic effects on the matter bispectrum R(k). The upper two rows show the lower and upper 68% percentiles of ΔR(k), and the two lower rows show the lower and upper 95% percentiles of ΔR(k). Each column is for a different k3 value, which is by construction larger than k2 and k1 but smaller than k1 + k2. The black circles show the actual location where we have measured B(k1, k2, k3). The background is coloured using a linear interpolation/extrapolation of the percentiles between the measured k-values.

Remarkably, the overall accuracy of the emulator ΔR(k) is better than 2% for the 68% percentile. In the 95% percentile, the accuracy of isolated triangle configurations seems to be worse than 5%. We have checked that this error arises from residual noise in the measurements. Given that we made sure that our emulator does not overfit this noise, we expect its accuracy to be better than the quoted values in those regions. We stress that these values refer to the accuracy of the emulation process. The underlying physical model was found to reproduce hydrodynamical simulations at a level of 1–3% in Aricò et al. (2021b) down to scales of 5 h Mpc⁻¹. In the next sections, we explore whether this accuracy is enough to provide an unbiased cosmological inference in a Euclid DR1 lensing analysis. For completeness, we show in Figs. C.1 and C.2 how S(k) and R(k) depend on BCM and cosmological parameters.

5. From spectra to weak lensing statistics

This section reviews the second- and third-order weak lensing statistics used to perform a forecast of a non-tomographic analysis of Euclid DR1. For a detailed review of gravitational lensing statistics, we refer to Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010), and Kilbinger (2015). In this work, we follow the description in Burger et al. (2023) and repeat only the essential details.

5.1. Limber projections of the power spectrum and bispectrum

The projected power- and bispectrum at wavenumber, ℓ, can be computed using the extended Limber approximation (Limber 1954; Kaiser & Jaffe 1997; Bernardeau et al. 1997; Schneider et al. 1998; LoVerde & Afshordi 2008), up to the maximum distance defined by source redshift distribution (Eq. (6)),

$$\begin{array}{cc}\hfill P_{\kappa \kappa }(\ell ;\mathrm{\Theta })& ={\displaystyle {\int }_0^{{\chi }_{\max }}}\mathrm{d}\chi \frac{g^2(\chi )}{a^2(\chi )}P[\frac{\ell +1/2}{\chi },z(\chi );\mathrm{\Theta }],\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill B_{\kappa \kappa \kappa }({\ell }_1,{\ell }_2,{\ell }_3;\mathrm{\Theta })& ={\displaystyle {\int }_0^{{\chi }_{\max }}}\mathrm{d}\chi \frac{g^3(\chi )}{a^3(\chi )\chi }\hfill \\ \hfill & \times B[\frac{{\ell }_1+1/2}{\chi },\frac{{\ell }_2+1/2}{\chi },\frac{{\ell }_3+1/2}{\chi },z(\chi );\mathrm{\Theta }],\hfill \end{array}$$(15)

where Θ = (Θ_cos, Θ_BCM) is the set of cosmological and BCM parameters. The 3-dimensional P(k, z; Θ) is defined as

$$\begin{array}{c}\hfill P(k,z;\mathrm{\Theta })=P^{\mathrm{GrO}}(k,z;{\mathrm{\Theta }}_{\cos })S^{\mathrm{BCM}}(k,z;{\mathrm{\Theta }}_{\mathrm{BCM}})\end{array}$$(16)

and P^GrO(k, z; Θ_cos) is the matter power spectrum computed using Halofit (Takahashi et al. 2012) that depends only on cosmological parameters Θ_cos, and S^BCM = P^Hydro/P^GrO is the BCM described in Sects. 3 and 4 that depends on the BCM parameters Θ_BCM. In principle, we can use any implementation of P^Gro, but we wanted to make use of k > 10 h Mpc⁻¹ and a < 0.4, which the current baccoemu does not provide. In analogy, the 3-dimensional bispectrum B(k, z, Θ) is computed as

$$\begin{array}{c}\hfill B(\mathit{k},z;\mathrm{\Theta })=B^{\mathrm{GrO}}(\mathit{k},z;{\mathrm{\Theta }}_{\cos })R^{\mathrm{BCM}}(\mathit{k},z;{\mathrm{\Theta }}_{\mathrm{BCM}}),\end{array}$$(17)

where B^GrO(k, z; Θ_cos) is the nonlinear matter bispectrum computed using BiHalofit (Takahashi et al. 2020), and R^BCM = B^Hydro/B^GrO is the BCM for the bispectrum.

Lastly, we are only considering closed triangles, ℓ₃ = |ℓ₁ + ℓ₂|, and g(χ) denotes the lensing efficiency and is defined as

$$\begin{array}{c}\hfill g(\chi )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}{\displaystyle {\int }_{\chi }^{{\chi }_{\max }}}\mathrm{d}\chi \prime n(\chi \prime )\frac{\chi \prime -\chi }{\chi \prime },\end{array}$$(18)

with n(χ) = n(z) dz/dχ and n(z) being the redshift probability distribution given in Eq. (6).

5.2. Aperture mass statistics

The aperture mass ℳ_ap at position ϑ with filter radius θ_ap is defined through the convergence κ, as follows:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_{\mathrm{ap}})={\displaystyle {\int }_{{\mathbb{R}}^2}}{\mathrm{d}}^2\vartheta \prime U_{{\theta }_{\mathrm{ap}}}(|\mathit{\vartheta }\mathbf{\prime }|)\kappa (\mathit{\vartheta }+\mathit{\vartheta }\mathbf{\prime }),\end{array}$$(19)

where U_θap(ϑ′) is a compensated filter such that ∫_ℝdϑ′ ϑ′ U_θap(ϑ′) = 0.

We define U_θap(ϑ) = θ_ap⁻²u(ϑ/θ_ap), denote by $\widehat{u}(\alpha )$ the Fourier transform of u, and use the filter function introduced in Crittenden et al. (2002),

$$\begin{array}{c}\hfill u(x)=\frac{1}{2\pi }(1-\frac{x^2}{2})exp(-\frac{x^2}{2}),\widehat{u}(\alpha )=\frac{{\alpha }^2}{2}exp(-\frac{{\alpha }^2}{2}).\end{array}$$(20)

5.3. Modelling aperture mass moments

While the expectation value of the aperture mass ⟨ℳ_ap⟩(θ_ap) vanishes by construction, the second-order (variance) of the aperture mass is nonzero and can be calculated as

$$\begin{array}{c}\hfill ⟨{\mathcal{M}}_{\mathrm{ap}}^2⟩({\theta }_{\mathrm{ap}})={\displaystyle {\int }_{{\mathbb{R}}_{+}}}\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\kappa \kappa }(\ell ){\widehat{u}}^2({\theta }_{\mathrm{ap}}\ell ).\end{array}$$(21)

Equivalently, the third-order moment of the aperture statistics ⟨ℳ_ap³⟩, can be computed from the convergence bispectrum via (Jarvis et al. 2004; Schneider et al. 2005)

$$\begin{array}{c}\hfill ⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩({\theta }_{\mathrm{ap},1},{\theta }_{\mathrm{ap},2},{\theta }_{\mathrm{ap},3})={\displaystyle {\int }_{{\mathbb{R}}_{+}^2}}\frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}{\displaystyle {\int }_{{\mathbb{R}}_{+}^2}}\frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\\ \hfill \times B_{\kappa \kappa \kappa }({\ell }_1,{\ell }_2,{\ell }_3)\widehat{u}({\theta }_{\mathrm{ap},1}{\ell }_1)\widehat{u}({\theta }_{\mathrm{ap},2}{\ell }_2)\widehat{u}({\theta }_{\mathrm{ap},3}{\ell }_3),\end{array}$$(22)

where ℓ₃ = |ℓ₁ + ℓ₂|. In the two lower panels of Fig. C.7, we show how ⟨ℳ_ap²⟩ and ⟨ℳ_ap³⟩ react to changes in the baryonic parameters, while all other parameters are fixed to the fiducial parameters given in Tables 1 and 4. We divide the ⟨ℳ_apⁿ⟩ predictions by GrO predictions, ⟨ℳ_apⁿ⟩^GrO, to highlight the baryonic effects. As expected, M_c has the largest impact on both ⟨ℳ_ap²⟩ and ⟨ℳ_ap³⟩. This baryonic parameter sets the characteristic halo mass for which half of the gas mass is ejected by feedback, thus regulating the quantity of gas retained in haloes. As expected, the galaxy parameter M₁ and the inner gas slope scale θ_inn have a larger effect on small apertures, whereas the feedback scale η has a larger impact on larger apertures. We note that, comparatively, β has a significantly larger effect on ⟨ℳ_ap³⟩ than on ⟨ℳ_ap²⟩. Looking at the cosmological parameters, we find that the dependence on σ₈ is negligible for ⟨ℳ_ap³⟩, but we observe some effect in the second-order statistics. We note that the dependence on Ω_m and Ω_b can be better captured by their ratio, Ω_b/Ω_m, which comes from the fact that the impact of the baryon feedback, at fixed strength, depends to zeroth order on the relative amounts of baryons to the dark matter.

5.4. Modelling convergence correlation functions

Since the ⟨ℳ_ap²⟩ statistics smooth much of the small scales, we decided to add the information from the two-point correlation function ξ_κ itself, which is defined by

$$\begin{array}{c}\hfill {\xi }_{\kappa }(\theta )={\displaystyle {\int }_{{\mathbb{R}}_{+}}}\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\kappa \kappa }(\ell )J_0(\theta \ell ),\end{array}$$(23)

where J₀ is the zeroth Bessel function of the first kind. In the lower panel of Fig. C.7, we show how ξ_κ reacts to changes in the baryonic parameters. As for the ⟨ℳ_apⁿ⟩ moments, we divide the ξ_κ predictions by GrO predictions, ξ_κ^GrO, to highlight the role of baryons. Similarly as for ⟨ℳ_ap²⟩, we observe that the most important baryonic parameters are M_c and η.

5.5. Flamingofied data vector

In this section, we describe the process of ‘flamingofying’ a GrO reference model. A flamingofied reference model of ⟨ℳ_apⁿ⟩ is built by scaling the prediction of a GrO model at the cosmology of the corresponding FLAMINGO simulation with the ratio of a hydrodynamic-GrO pair of FLAMINGO simulations

$$\begin{array}{c}\hfill {⟨{\mathcal{M}}_{\mathrm{ap}}^n⟩}^{\mathrm{flamingofied}}={⟨{\mathcal{M}}_{\mathrm{ap}}^n⟩}^{\mathrm{GrO}}\frac{{⟨{\mathcal{M}}_{\mathrm{ap}}^n⟩}^{\mathrm{FLAMINGO},\mathrm{Hydro}}}{{⟨{\mathcal{M}}_{\mathrm{ap}}^n⟩}^{\mathrm{FLAMINGO},\mathrm{GrO}}}·\end{array}$$(24)

Analogously, we flamingofy the ξ_κ in the following way

$$\begin{array}{c}\hfill {\xi }_{\kappa }^{\mathrm{flamingofied}}={\xi }_{\kappa }^{\mathrm{GrO}}\frac{{\xi }_{\kappa }^{\mathrm{FLAMINGO},\mathrm{Hydro}}}{{\xi }_{\kappa }^{\mathrm{FLAMINGO},\mathrm{GrO}}}·\end{array}$$(25)

In Figs. 3 and 4, we show in blue the ratio Hydro/GrO of the FLAMINGO measurements for ⟨ℳ_apⁿ⟩ and ξ_κ for the L2p8_m9 scenario. Comparing these ratios to the estimated uncertainty in grey indicates which scales are affected by baryonic effects. Furthermore, we show the red dashed line as the best-fitting model to fit the flamingofied reference model if we fix the cosmology to D3A, where we performed a χ² minimisation for each statistic separately to find the best-fitting parameters. Our BCM model can perfectly fit the flamingofied reference model. This is an encouraging result for the cosmological analysis in the next section.

Fig. 3. Baryonic effects on ⟨ℳap2⟩ (upper panel) and on ξκ (lower panel), illustrated by ratios of measurements with baryonic physics to GrO measurements. We show with blue triangles the fiducial FLAMINGO measurements, and the grey shaded area represents the Euclid DR1 uncertainty. We overplot with a red dashed line the best-fitting BCM to the flamingofied data vector, found by employing our emulators with fixed cosmological parameters.

Fig. 4. Same as Fig. 3 but for the ⟨ℳap3⟩ statistics. The x-axis shows different combinations of θap, 1-θap, 2-θap, 3.

We also show the effects of baryons from the FLAMINGO description L2p8_m9 in each panel. The ⟨ℳ_apⁿ⟩ in the FLAMINGO and Takahashi simulations are measured by applying Eq. (19), using the Healpy smoothing option. We have decided to use aperture filter radii θ_ap = {2′, 4′, 8′, 16′, 32′}, where the lower value θ_ap = 2′ is constrained by the chosen pixel resolution $0\stackrel{\prime }{.}43$ of the FLAMINGO/T17 convergence maps. The ξ_κ in the FLAMINGO and Takahashi simulations are measured with TreeCorr (Jarvis et al. 2004), from $0\stackrel{\prime }{.}5<\theta <300\stackrel{\prime }{.}$ in 15 logarithmic bins, where the lower edge is determined again by the resolution of the chosen convergence maps.

6. Euclid forecasts

In this section, we perform a non-tomographic Euclid DR1 forecast, where we vary the matter density Ω_m, the baryon density Ω_b, the normalisation of the power spectrum σ₈, and the five baryonic parameters described in Sect 2.1.2. The dependence of these parameters on ⟨ℳ_apⁿ⟩ and the ξ_κ are already discussed in Sect. 5 and are shown in Fig. C.7. We have not yet mentioned that the computations of the integrals of these quantities are not efficient enough to perform several Markov Chain Monte Carlo (MCMC) analyses. Therefore, we decided to train another set of emulators on top of the three-dimensional power and bispectrum emulators, which predict the ⟨ℳ_apⁿ⟩ and ξ_κ. We computed the ⟨ℳ_apⁿ⟩ and the ξ_κ for 2200 points distributed in a Latin hypercube. We estimate the error of the emulation, by building a set of emulators iteratively leaving out from the training 200 of the 2200 points at a time. Since the power and bispectrum are noisy quantities themselves, their noise also propagates into the summary statistics. Furthermore, the integration to ⟨ℳ_ap³⟩ sometimes fails for very squeezed configurations, because we had to limit the number of ℓ values where we compute B_κκκ for the sake of computational speed. Therefore, the emulator error accounts for emulator noise in the summary statistics but also for noise in the power and bispectrum estimations and their subsequent integration. We note that this procedure, however, does not account for possible systematic errors in the emulation or integration. The error of this emulator is less than a tenth of the uncertainty estimated with the Takahashi et al. (2017) simulations. Although this might be neglected, we decided to build an emulator covariance matrix for correctly estimating the posterior, using the 2200 Δm = m^emu − m^true predictions and add it to the Takahashi et al. (2017) covariance matrix as

$$\begin{array}{c}\hfill \mathsf{C}={\mathsf{C}}^{\mathrm{T}17}+{\mathsf{C}}^{\mathrm{emu}}.\end{array}$$(26)

Since the estimated covariance matrix C is a random variable itself, we follow Percival et al. (2022) to compute the Likelihood. Given a data vector d and covariance matrix C of rank n_d measured from n_r realizations, the posterior distribution of a model vector m(Θ) that depends on n_Θ = 12 parameters is

$$\begin{array}{c}\hfill \mathit{P}(\mathit{m}(\mathrm{\Theta })|\mathit{d},\mathsf{C})\propto {|\mathsf{C}|}^{-\frac{1}{2}}{(1+\frac{{\chi }^2}{n_{\mathrm{r}}-1})}^{-m/2},\end{array}$$(27)

where

$$\begin{array}{c}\hfill {\chi }^2={[\mathit{m}(\mathrm{\Theta })-\mathit{d}]}^{\mathrm{T}}{\mathsf{C}}^{-1}[\mathit{m}(\mathrm{\Theta })-\mathit{d}].\end{array}$$(28)

The power-law index m is

$$\begin{array}{c}\hfill m=n_{\mathrm{\Theta }}+2+\frac{n_{\mathrm{r}}-1+B(n_{\mathrm{d}}-n_{\mathrm{\Theta }})}{1+B(n_{\mathrm{d}}-n_{\mathrm{\Theta }})},\end{array}$$(29)

and B is defined as

$$\begin{array}{c}\hfill B=\frac{n_{\mathrm{r}}-n_{\mathrm{d}}-2}{(n_{\mathrm{r}}-n_{\mathrm{d}}-1)(n_{\mathrm{r}}-n_{\mathrm{d}}-4)}·\end{array}$$(30)

and n_r = 570 is the number of realizations, n_Θ = 8 the number of free parameters, and n_d the rank of the covariance matrix. Although n_r for C^emu is not clearly defined, we used the same correction for C as we used for C^T17. Lastly, we used the nested sampler nautilus (Lange 2023) to perform the inference, employing 2000 live points. We validated the convergence of the posteriors by comparing them with a Hamiltonian Monte Carlo sampler and a standard Metropolis–Hastings algorithm.

To start our investigation of the impact of baryonic feedback, we show in Fig. 5 posteriors on the cosmological parameters S₈, Ω_m, and Ω_b, which we obtain by analysing a synthetic Euclid DR1 data vector built from a combination ξ_κ, ⟨ℳ_ap²⟩, and ⟨ℳ_ap³⟩. The investigation here is built solely using Halofit, BiHalofit, and the BCM, which means we always know the true input of cosmological and BCM parameters. In the first case (black), we created a data vector using the BCM with the parameters given in Table 4 and analysed it with the GrO model, meaning we turned the BCM off in the modelling. This results in a 2.5σ shift in S₈. Substantial scale cuts on small scales are required to avoid such a bias, which we show as the orange contours. In terms of scales these are θ > 8′ for ξ_κ, θ_ap = 32′ for ⟨ℳ_ap²⟩ and θ_ap, 1 ∈ {16′, 32′} with θ_ap, 2, 3 = 32′ for ⟨ℳ_ap³⟩. However, this drastically reduces the precision on S₈ by 8% and by almost 80% on Ω_m compared to modelling all five BCM parameters shown as the red contours. We obtain even tighter constraints by fixing these five BCM parameters to their true values, as shown in the blue contours. Having perfect knowledge of baryonic effects results in cosmological constraints almost identical to the case where we used a GrO data vector and analysed it with a GrO model, displayed with purple contours. We will show later that we probably do not need to vary all five BCM parameters for the expected Euclid DR1 uncertainty.

Fig. 5. Posteriors of the cosmological parameters S8, Ωm, and Ωb, analysing a synthetic Euclid DR1 joint data vector of second- and third-order weak lensing statistics. We consider several cases. First, a gravity-only (GrO) data vector was analysed with a GrO model (purple contours). We obtain the black contours when baryonic effects are included in the data vector but not in the model. In orange, we illustrate the constraints we get when all elements of the data vector are excluded where the baryonic data vector deviates by more than 0.4σ from the GrO data vector. Lastly, we show the results if the model includes baryonic feedback in red and blue, where we fixed the baryon parameters to the truth in the blue case and marginalised over all five parameters for the red. The grey dashed line indicates the underlying D3A cosmology.

6.1. Validating the emulator on the FLAMINGO simulations

First, we use a flamingofied data vector estimated from eight light cones of the 2.8 Gpc box that uses the fiducial baryon description (see Sect. 5.5). This test aims to see if the combination of the three different statistics improves the constraining power and check if the BCM is flexible enough to fit the fiducial FLAMINGO simulation. The constraints on both cosmological and baryonic parameters are shown in Fig. 6. This is one of the main results of this paper and shows that our baryon model can describe the modifications due to baryons while resulting in unbiased cosmological results for all three statistics and their joint analysis. We note that ⟨ℳ_ap²⟩ yields similar constraints as ξ_κ, since both result from integrating the power spectrum, so we decided to show only the combination. Next, we notice that ⟨ℳ_ap³⟩ constrains the cosmological parameters slightly better than the second-order statistics, and furthermore helps breaking parameter degeneracies. Compared to Burger et al. (2024), where second-order statistics were relatively more constraining, we expect that this comes from the fact that we include smaller scales with θ_ap = 2′, where higher-order statistics become more important. We do not show the ‘true’ FLAMINGO values for the BCM parameters. These can be measured in hydrodynamic simulations in principle; however, their measured values are redshift-dependent. Thus, interpreting our constraints is not straightforward since they are weighted in redshift according to our lensing kernel. We leave for future work a thorough test of the BCM parameter constraints, as a correct inference of baryonic physics is essential, for example, in cross-correlations and multi-probes analyses.

Fig. 6. Euclid DR1 forecast, where we used a flamingofied reference data vector (see Sect. 5.5) estimated from the fiducial baryon description of the FLAMINGO simulations with eight light cones of the 2.8 Gpc box. The stated constraints correspond to the joint analysis of second- and third-order statistics. The grey dashed line indicates the underlying D3A cosmology.

In Fig. 7, we compare the constraining power of second-order and third-order statistics in the Ω_m-$S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ plane. We highlight the different degeneracy axes of the second- and third-order statistics, which are broken when both are combined, in agreement with Semboloni et al. (2013). Therefore, combining second- and third-order statistics increases the constraining power by 22% for S₈ and even 50% for Ω_m compared to using second-order statistics alone. Lastly, we also show the case where we discard all the non-equal scale filter radii θ_ap. Interestingly, these non-equal scale filter radii θ_ap play an essential role in the precision and accuracy of cosmological parameters, since discarding them reduces the constraining power by around 16% for S₈ and Ω_m compared to using all filter combinations. Burger et al. (2023) found that these non-equal scale filter radii are less significant. We interpret this as non-equal scale filter radii being important to constrain the BCM parameters, which in turn helps to constrain cosmological parameters. Another explanation for the discrepancies with Burger et al. (2023) could also lie in the fact that we include the θ_ap = 2′ filter and that we used a different tomographic setup.

Fig. 7. Euclid DR1 forecast, where we used a flamingofied reference data vector (see Sect. 5.5) estimated from the fiducial baryon description of the FLAMINGO simulations with eight light cones of the 2.8 Gpc box. The stated constraints correspond to the joint analysis of second- and third-order statistics. The posteriors are the same as those used in Fig. 6, meaning we marginalised over all BCM parameters. The stated constraints correspond to the joint analysis of second- and third-order statistics (blue contours). The grey dashed line indicates the underlying D3A cosmology.

Testing only one astrophysical scenario is insufficient to ensure a robust cosmological inference. Therefore, we use all baryon feedback descriptions described in Sect. 2.2. We build flamingofied data vectors as described in Sect. 5.5. We run an MCMC for each baryon feedback description once using only second-order statistics ξ_κ + ⟨ℳ_ap²⟩ and once combining second- and third-order ξ_κ + ⟨ℳ_ap²⟩+⟨ℳ_ap³⟩, where we vary all baryon parameters. We present the result as the black symbols in Fig. 8, where the error bars represent the 1σ uncertainties estimated from the MCMC chains. The corresponding best-fitting parameters are shown as the stars, and we notice that the best-fitting parameters are much closer to the truth than the mean, which might indicate non-Gaussian parameter posteriors or significant projection effects. Remarkably, we get unbiased cosmological results for all FLAMINGO flavours and configurations for both the second-order alone and the second- and third-orders combined. This is the most important result of this paper, as it shows that our BCM model is flexible enough to get unbiased results for all FLAMINGO feedback variations. This also shows that the interpolation and extrapolation of our emulators for the BCM power spectrum and bispectrum are sufficient for the current setup. We expect the emulator accuracy requirements to be similar in a tomographic setup, such that our emulators are ready for a Euclid-like DR1 analysis. However, a tomographic analysis might require an explicit redshift dependence on the baryonic parameters. We highlight that the three-dimensional power spectrum and bispectrum emulators developed in this work would remain valid and could be used with parametrised redshift dependencies without creating new training sets or emulators.

Fig. 8. Estimated weighted mean (triangles) and 68% credible intervals from the MCMC chains for all FLAMINGO models described in Sect. 2.2. Since the LS8 and Planck cosmology is different we plot ΔΩm = Ωmbest − Ωmtrue and in analogy ΔS8 and ΔΩb. The different colours show cases where parameters are fixed. The stars indicate the best-fitting parameters resulting from a minimisation process that returns the lowest χ2. The figure is for ξκ + ⟨ℳap2⟩+⟨ℳap3⟩. The corresponding figure for ξκ + ⟨ℳap2⟩ is shown in Fig. C.6.

6.2. Assessing the importance of the BCM parameters

In this subsection, we exploit a wide range of feedback models of the FLAMINGO suite to assess the importance of each BCM parameter in a cosmological analysis. We aim to find the minimal BCM model required to deliver unbiased cosmological constraints in a Euclid DR1 lensing analysis. Fewer free parameters would help mitigate projection effects and allow for faster convergence of the posteriors. This analysis cannot be straightforwardly generalised to different hydrodynamical simulations and summary statistics. Nevertheless, this section provides us with valuable information about the most important BCM parameters that can capture the large variety of baryonic variations in FLAMINGO, even if care must be taken when fixing one or more baryonic parameters in data analyses.

To perform our test, we fix one or more BCM parameters to the weighted average of the best-fitting parameters over all FLAMINGO flavours and configurations, where the weights are estimated from the 1σ credible intervals of MCMC chains. We report those values in Table 5. Next, we run an MCMC for each scenario, and show the best-fitting values, mean and standard deviation of the marginalised posteriors in Fig. 8. Focusing first on ξ_κ + ⟨ℳ_ap²⟩+⟨ℳ_ap³⟩, we notice that the best-fitting values, especially for S₈, are closer to the truth than the weighted mean of the MCMC, indicating strong projection effects. We find that fixing M_c is particularly bad in some cases, where an incorrect estimation of M_c results in a biased S₈ inference. Next, we observe that fixing η leads to small biases for the strongest feedback descriptions (jet fgas–4σ and fgas–8σ). Fixing either M₁ or β does result in unbiased results for all scenarios, which is not surprising given that M₁ and β are unconstrained, as we can see in Fig. 6. Fixing η or θ_inn also results in small biases for the LS8 cosmology. Next, we test fixing all BCM parameters except M_c and η. We get unbiased results for all FLAMINGO flavours, except the Jet flavour shows 1σ shift for the Ω_m parameter. When fixing all the BCM except M_c, the effect on the Jet feedback description is enhanced. We note that varying only one or two BCM parameters results in much better-constrained cosmological parameters. If the baryonic feedback inferred from future observations will be consistent with the FLAMINGO descriptions that were tested here, it is probably sufficient for current Stage-III and first data releases of Stage-IV surveys to vary only M_c and η, or even only M_c. Lastly, we notice that varying parameters such as M₁ might not influence the estimation of cosmological parameters, but can play a vital role in constraining other BCM, such as M_c.

Focusing now on the difference between ξ_κ + ⟨ℳ_ap²⟩+⟨ℳ_ap³⟩ and ξ_κ + ⟨ℳ_ap²⟩, which we show in Fig. C.6, we see a huge difference in constraining power, particularly for Ω_m and some of the BCM parameters such as M_c. We note here that for ξ_κ + ⟨ℳ_ap²⟩ we use the best-fitting parameters obtained from ξ_κ + ⟨ℳ_ap²⟩ itself, which are also listed in Table 5. Furthermore, we observe an enhancement in the scatter of the best-fitting parameters for ξ_κ + ⟨ℳ_ap²⟩, which indicates stronger projection effects. The reason is that the ⟨ℳ_ap³⟩ part of the data vector helps to constrain the BCM parameters, which in turn reduces the projection effects. Lastly, we notice the gain in fixing all BCM parameters except M_c is greater if we include ⟨ℳ_ap³⟩ to the data vector. The reason is that ⟨ℳ_ap³⟩ significantly helps to constrain the M_c parameter, which in turn helps to constrain cosmological parameters.

7. Conclusions

This paper presents novel emulators for modelling the baryonic effects on the matter bispectrum and power spectrum, developed using high-resolution cosmological simulations. These emulators are designed to accurately predict the impact of baryonic processes such as gas cooling, star formation, and feedback from star formation and AGN on the matter power spectrum and bispectrum across a wide range of scales, redshifts, and cosmological scenarios. The ability to model these effects is crucial for large-scale structure surveys, such as Euclid, which aims to extract cosmological information from the nonlinear regime of structure formation.

To construct the emulators, we exploit the high-resolution BACCON-body simulations, which were post-processed to explore different cosmological scenarios using a cosmology rescaling algorithm. Since direct hydrodynamical simulations are computationally expensive, we incorporated baryonic effects such as gas cooling, star formation, and AGN feedback using the baryonification framework. This approach modifies the matter distribution in gravity-only simulations, effectively mimicking the impact of baryonic processes while maintaining computational efficiency.

Using a GPU-accelerated estimator, we then measured the matter power spectrum P(k) and bispectrum B(k₁, k₂, k₃) across a wide range of scales down to k_max ≲ 20 h Mpc⁻¹, with shot noise contribution carefully subtracted. The resulting measurements formed the training set of our emulators.

We trained deep neural networks to emulate the baryon correction model of the power spectrum, S(k) = P_Hydro(k)/P_GrO(k), and bispectrum R(k₁, k₂, k₃) = B_Hydro(k₁, k₂, k₃)/B_GrO(k₁, k₂, k₃), where ‘Hydro’ refers to hydrodynamical simulations and ‘GrO’ to Gravity-only simulations. The training dataset consisted of 1200 simulations, sampled using a Latin Hypercube across five baryonic and three cosmological parameters.

Our emulator is flexible enough to fit the state-of-the-art hydrodynamic FLAMINGO simulations and get unbiased cosmological results. It successfully captures the baryonic suppression of second- and third-order weak lensing statistics across various configurations and scales. This validation underscores the reliability of the emulator as a computationally efficient alternative to running full hydrodynamical simulations, which are prohibitively expensive for exploring the extensive parameter space required by surveys such as Euclid. We note that the expected accuracy of cosmology scaling+baryonification for the power spectrum is approximately 1% and 3% for the bispectrum, such that the emulator serves as an alternative to hydrodynamical simulations as long as the statistical uncertainty of the data dominates over the modelling error.

Using the emulator, we provided forecasts for the impact of baryonic effects on the matter power and bispectrum measured by Euclid. Our results highlight the importance of incorporating baryonic corrections when analysing future LSS data. Neglecting these effects could either lead to biased constraints on cosmological parameters of the order of 3σ in S₈ or to a significant reduction of the constraining power by up to 80% for Ω_m, if the scales affected by baryons are excluded. Both cases would be fatal for the high-precision measurements Euclid is aiming to deliver. By modelling these effects accurately, our emulator allows for more robust cosmological analyses and ensures that Euclid can achieve its scientific objectives.

We plan to assess the emulators’ accuracy against currently extrapolated scales in future work. This requires larger simulations and more resources to measure spectra from large grids. We also plan to employ our baryonification emulators to extend the forecast presented in this work to a mock Euclid tomographic analysis. Finally, it would be interesting to validate with hydrodynamical simulations not only the constraints on cosmological parameters, but also the inferred baryonic parameters.

Data availability

The emulators developed in this work are publicly available at https://baccoemu.readthedocs.io/en/latest/.

Acknowledgments

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMIMI, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org/consortium/community/). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This research was enabled in part by support provided by Calcul Quebec (https://docs.alliancecan.ca/wiki/Narval/en) and Compute Ontario (https://docs.alliancecan.ca/wiki/Graham) and the Digital Research Alliance of Canada (https://www.alliancecan.ca/en). This research used resources from the National Energy Research Scientific Computing Center, supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. PB acknowledge financial support from the Canadian Space Agency (Grant 23EXPROSS1), the Waterloo Centre for Astrophysics and the National Science and Engineering Research Council Discovery Grants program. LL is supported by the Austrian Science Fund (FWF) [ESP 357-N]. REA acknowledges sup- port from project PID2021-128338NB-I00 from the Spanish Ministry of Science and support from the European Research Executive Agency HORIZON-MSCA- 2021-SE-01 Research and Innovation programme under the Marie Skłodowska- Curie grant agreement number 101086388 (LACEGAL) MZ is supported by STFC. LP acknowledges support from the DLR grant 50QE2302. We acknowledge the use of the following software: NumPy (Harris et al. 2020), Numba (Lam et al. 2015), SciPy (Virtanen et al. 2020), JAX (Bradbury et al. 2021), TreeCorr (Jarvis et al. 2004), CosmoPower (Spurio Mancini et al. 2022), HEALPix (Górski et al. 2005).

References

Appendix A: Accuracy of bispectrum measurement code

To validate our bispectrum measurement code, we compare it to the established code bskit (Scoccimarro 2000; Foreman et al. 2020), which uses nbodykit (Hand et al. 2018). It implements the same algorithm as BiG (Linke 2024) but uses no GPU acceleration. We show in Fig. A.1 the relative difference of the BiG and bskit estimates for the bispectrum of an example simulation with boxsize L = 512h⁻¹Mpc and grid size N = 512. The difference is at the numerical precision level, with no deviation being larger than 2 × 10⁻⁴.

Fig. A.1. Relative difference of bispectrum estimate by our measurement code BiG and bskit for an example simulation. Left: For equilateral triangles, Right: for all triangle configurations

Appendix B: Additional accuracy plots of the bispectrum correction model

This section shows two additional plots that illustrate the accuracy of BCM bispectrum emulator. While in Fig. B.1 we show the flattened and concatenated measurements from all three folding, we show in Fig. B.2 the same as a three-dimensional scatter plot.

Fig. B.1. Accuracy of the emulator for baryonic effects on the matter bispectrum. The y-axis is defined as ΔR(k) = Remu/Rtrue − 1, where R(k) = BBCM(k)/BGrO(k). Here, we flattened and concatenated the three folding. The horizontal dashed lines show the 2% region.

Fig. B.2. Similar to Fig. B.1, but now illustrated as a three dimensional scatter plot.

Appendix C: Parameter dependence

In this appendix, we explore further the dependences of two- and three-point statistics on both cosmological and baryonic parameters. We begin in the upper panel of Fig. C.1, with the dependence of the power spectrum ratio, S(k) = P^Hydro/P^GrO, on baryonic feedback effects. In the lower panels of Figs. C.1 and C.2, we then show the dependence of bispectra ratio, R(k) = B^Hydro/B^GrO for three different k configurations. Here we especially see that setting R(k₁, k₂, k₃) = 1 if one k_i < k_min might result in biased results.

In Figs. C.3 and C.4, we display the dependence of power spectrum and bispectrum on the BCM parameter θ_out and M_inn. We apply the same baryonification, compatible with the FLAMINGO suite, to pairs of N-body simulations which vary one cosmological parameter at a time, in a wide range displayed in the legend. We note that, for simplicity, we use a simple linear interpolation and do not remove noisy measurements. However, this does not affect our conclusions. The dependence of θ_out and M_inn is significantly smaller compared to other BCM parameters, such as M₁ or β. Given their subdominant contribution, we decided to fix θ_out = 1 and M_inn = 2.3 × 10¹³ h⁻¹ M_⊙.

Furthermore, we show in Fig. C.5 the dependence of S(k) and R(k) on cosmological parameters. To do so, we apply the baryonification with a fixed set of baryonic parameters to a set of N-body simulations, where one cosmological parameter at a time is varied. We find that, analogously to the results with the power spectrum obtained by Aricò et al. (2021a), Ω_m and Ω_b are the most significant, followed by σ₈. The other parameters, including h, n_s, Σm_ν, w₀, w_a, have a smaller impact than the nominal accuracy of the cosmological scaling, and thus we do not model them. Lastly, we show in Fig. C.6 the equivalent to Fig. 8 but for second-order statistics alone, and in Fig. C.7 the dependence of the three measured weak lensing statistics on the effects of baryonic feedback.

Fig. C.1. In this figure, we show the power spectrum’s dependence on baryonic feedback effects in the three upper rows. In the bottom three rows, we show the dependence of the bispectrum for k1 = k2 = k3.

Fig. C.2. The same as the bottom three rows of Fig. C.1 but for two different settings of ki as indicated on the y-axis.

Fig. C.3. Dependence of the S(k) and R(k) on the BCM parameter θout. Yellow lines are smaller and blue lines are larger θout.

Fig. C.4. Dependence of the S(k) and R(k) on the BCM parameter Minn. Yellow lines are smaller and blue lines are larger Minn.

Fig. C.5. Dependence of the S(k) and R(k) on cosmological parameters. Here, we plot the ratio between each parameter’s lowest and largest values, as indicated in the legend.

Fig. C.6. Same as Fig. 8 but second-order statistics alone. The triangles with error bars are estimated from the MCMC chains for all FLAMINGO models described in Sect. 2.2. Since the LS8 and Planck cosmology is different we plot ΔΩm = Ωmbest − Ωmtrue and in analogy ΔS8 and and ΔΩb. The different colours show cases where parameters are fixed. The stars indicate the best-fitting parameters resulting from a minimisation process that returns the lowest χ2. The corresponding figure for ξκ + ⟨ℳap2⟩+⟨ℳap3⟩ is shown in Fig. 8.

Fig. C.7. Dependence of the ⟨ℳapn⟩ and ξκ to changes of baryonic and cosmological parameters. We scaled the model vectors using the same model just for a GrO case, disentangling the effect of the baryons. The θap = {θap, 1, θap, 2, θap, 3} values for the lowest panel are increasing from left to right as in Fig. 4, increasing first θap, 3, then θap, 2 and lastly θap, 1.

All Tables

All Figures

	Fig. 1. Accuracy of the emulator for baryonic effects on the matter power spectrum. The y-axis is defined as ΔS(k) = Semu/Strue − 1, where S(k) = PBCM(k)/PGrO(k). The horizontal dashed lines show the 0.5% region.
In the text

Fig. 2. Accuracy of the emulator for baryonic effects on the matter bispectrum R(k). The upper two rows show the lower and upper 68% percentiles of ΔR(k), and the two lower rows show the lower and upper 95% percentiles of ΔR(k). Each column is for a different k3 value, which is by construction larger than k2 and k1 but smaller than k1 + k2. The black circles show the actual location where we have measured B(k1, k2, k3). The background is coloured using a linear interpolation/extrapolation of the percentiles between the measured k-values.

In the text

Fig. 3. Baryonic effects on ⟨ℳap2⟩ (upper panel) and on ξκ (lower panel), illustrated by ratios of measurements with baryonic physics to GrO measurements. We show with blue triangles the fiducial FLAMINGO measurements, and the grey shaded area represents the Euclid DR1 uncertainty. We overplot with a red dashed line the best-fitting BCM to the flamingofied data vector, found by employing our emulators with fixed cosmological parameters.

In the text

	Fig. 4. Same as Fig. 3 but for the ⟨ℳap3⟩ statistics. The x-axis shows different combinations of θap, 1-θap, 2-θap, 3.
In the text

Fig. 5. Posteriors of the cosmological parameters S8, Ωm, and Ωb, analysing a synthetic Euclid DR1 joint data vector of second- and third-order weak lensing statistics. We consider several cases. First, a gravity-only (GrO) data vector was analysed with a GrO model (purple contours). We obtain the black contours when baryonic effects are included in the data vector but not in the model. In orange, we illustrate the constraints we get when all elements of the data vector are excluded where the baryonic data vector deviates by more than 0.4σ from the GrO data vector. Lastly, we show the results if the model includes baryonic feedback in red and blue, where we fixed the baryon parameters to the truth in the blue case and marginalised over all five parameters for the red. The grey dashed line indicates the underlying D3A cosmology.

In the text

	Fig. 6. Euclid DR1 forecast, where we used a flamingofied reference data vector (see Sect. 5.5) estimated from the fiducial baryon description of the FLAMINGO simulations with eight light cones of the 2.8 Gpc box. The stated constraints correspond to the joint analysis of second- and third-order statistics. The grey dashed line indicates the underlying D3A cosmology.
In the text

Fig. 7. Euclid DR1 forecast, where we used a flamingofied reference data vector (see Sect. 5.5) estimated from the fiducial baryon description of the FLAMINGO simulations with eight light cones of the 2.8 Gpc box. The stated constraints correspond to the joint analysis of second- and third-order statistics. The posteriors are the same as those used in Fig. 6, meaning we marginalised over all BCM parameters. The stated constraints correspond to the joint analysis of second- and third-order statistics (blue contours). The grey dashed line indicates the underlying D3A cosmology.

In the text

Fig. 8. Estimated weighted mean (triangles) and 68% credible intervals from the MCMC chains for all FLAMINGO models described in Sect. 2.2. Since the LS8 and Planck cosmology is different we plot ΔΩm = Ωmbest − Ωmtrue and in analogy ΔS8 and ΔΩb. The different colours show cases where parameters are fixed. The stars indicate the best-fitting parameters resulting from a minimisation process that returns the lowest χ2. The figure is for ξκ + ⟨ℳap2⟩+⟨ℳap3⟩. The corresponding figure for ξκ + ⟨ℳap2⟩ is shown in Fig. C.6.

In the text

	Fig. A.1. Relative difference of bispectrum estimate by our measurement code BiG and bskit for an example simulation. Left: For equilateral triangles, Right: for all triangle configurations
In the text

	Fig. B.1. Accuracy of the emulator for baryonic effects on the matter bispectrum. The y-axis is defined as ΔR(k) = Remu/Rtrue − 1, where R(k) = BBCM(k)/BGrO(k). Here, we flattened and concatenated the three folding. The horizontal dashed lines show the 2% region.
In the text

	Fig. B.2. Similar to Fig. B.1, but now illustrated as a three dimensional scatter plot.
In the text

	Fig. C.1. In this figure, we show the power spectrum’s dependence on baryonic feedback effects in the three upper rows. In the bottom three rows, we show the dependence of the bispectrum for k1 = k2 = k3.
In the text

	Fig. C.2. The same as the bottom three rows of Fig. C.1 but for two different settings of ki as indicated on the y-axis.
In the text

	Fig. C.3. Dependence of the S(k) and R(k) on the BCM parameter θout. Yellow lines are smaller and blue lines are larger θout.
In the text

	Fig. C.4. Dependence of the S(k) and R(k) on the BCM parameter Minn. Yellow lines are smaller and blue lines are larger Minn.
In the text

	Fig. C.5. Dependence of the S(k) and R(k) on cosmological parameters. Here, we plot the ratio between each parameter’s lowest and largest values, as indicated in the legend.
In the text

Fig. C.6. Same as Fig. 8 but second-order statistics alone. The triangles with error bars are estimated from the MCMC chains for all FLAMINGO models described in Sect. 2.2. Since the LS8 and Planck cosmology is different we plot ΔΩm = Ωmbest − Ωmtrue and in analogy ΔS8 and and ΔΩb. The different colours show cases where parameters are fixed. The stars indicate the best-fitting parameters resulting from a minimisation process that returns the lowest χ2. The corresponding figure for ξκ + ⟨ℳap2⟩+⟨ℳap3⟩ is shown in Fig. 8.

In the text

	Fig. C.7. Dependence of the ⟨ℳapn⟩ and ξκ to changes of baryonic and cosmological parameters. We scaled the model vectors using the same model just for a GrO case, disentangling the effect of the baryons. The θap = {θap, 1, θap, 2, θap, 3} values for the lowest panel are increasing from left to right as in Fig. 4, increasing first θap, 3, then θap, 2 and lastly θap, 1.
In the text