Simultaneous modelling of matter power spectrum and bispectrum in the presence of baryons
Giovanni Aricò, Raul E. Angulo, Carlos Hernández-Monteagudo, Sergio Contreras, Matteo Zennaro
MMNRAS , 1–14 (2020) Preprint 1 October 2020 Compiled using MNRAS L A TEX style file v3.0
Simultaneous modelling of matter power spectrum and bispectrumin the presence of baryons
Giovanni Aricò (cid:63) , Raul E. Angulo , , Carlos Hernández-Monteagudo , , Sergio Contreras , & Matteo Zennaro . Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal, 4, 20018, Donostia-San Sebastián, Guipuzkoa, Spain. IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain. Centro de Estudios de Física del Cosmos de Aragón, Unidad Asociada CSIC, Plaza San Juan 1, 44001 Teruel, Spain. Instituto de Astrofísica de Canarias, Spain. University of La Laguna, Spain.
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We demonstrate that baryonification algorithms, which displace particles in gravity-onlysimulations according to physically-motivated prescriptions, can simultaneously capture theimpact of baryonic physics on the 2 and 3-point statistics of matter. Specifically, we showthat our implementation of a baryonification algorithm jointly fits the changes induced bybaryons on the power spectrum and equilateral bispectrum on scales up to k = 5 h Mpc − andredshifts ≤ z ≤ , as measured in six different cosmological hydrodynamical simulations.The accuracy of our fits are typically ∼ for the power spectrum, and for the equilateraland squeezed bispectra, which somewhat degrades to ∼ for simulations with extremefeedback prescriptions. Our results support the physical assumptions underlying baryonifica-tion approaches, and encourage their use in interpreting weak gravitational lensing and othercosmological observables. Key words: large-scale structure of Universe – cosmological parameters – cosmology: theory
Despite large efforts of the scientific community, the nature of darkenergy and dark matter remains elusive. Even if the standard Λ CDMmodel has successfully passed many independent tests in the lastdecades, recent tensions in the estimated values of the Hubble con-stant and in the amplitude of the linear fluctuation have been pointedout as a possible window to physics beyond Λ CDM (e.g. Verde et al.2019; Wong et al. 2020). To successfully solve these tensions, it isparamount that current and upcoming cosmological surveys ex-tract the maximum amount of cosmological information at the lowredshifts, where dark energy and dark matter are more accessible(Planck Collaboration et al. 2018; Troxel et al. 2018; Benitez et al.2014; Laureijs et al. 2011; DESI Collaboration et al. 2016; Aiharaet al. 2018). For many observables and statistics, the limiting fac-tor will be the predictability and accuracy of theoretical modelsemployed to analyse the data.For the case of next-generation weak lensing surveys, thelargest theoretical uncertainty is given by baryonic physics – gascooling, star formation, and feedback, for instance, modify signif-icantly the total mass distribution in the universe in a way that is (cid:63)
E-mail:[email protected] (GA) not possible to accurately predict from first principles. On the otherhand, if these baryonic processes are modelled appropriately, thenwe could extract more cosmological information, and possibly alsoconstrain astrophysical processes.A promising approach to incorporate the baryonic effects inmodels for the cosmic density field is baryonification (Schneider &Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020). Briefly,these algorithms, a.k.a. baryon correction models (BCM), dis-place particles in gravity-only simulations according to physically-motivated recipes designed to mimic the effects produced bybaryons in the Large Scale Structure (LSS) of the universe. Thismethod has been extensively tested against many hydrodynamicalsimulations, and it is shown to be very accurate in capturing thechanges induced by baryons on the power spectrum (Schneideret al. 2020; Aricò et al. 2020).In general, baryons are expected to modify the full densityfield and thus the whole hierarchy of N -point functions, not onlythe power spectrum. Indeed, Foreman et al. (2019) recently showedthat baryonic effects on the bispectrum of hydrodynamical simula-tions are present, and that they carry extra information with respectto the power spectrum. Given the simplifications and assumptionsof baryonification methods – e.g. spherically symmetric displace-ments, dependences on halo mass, and neglected physical processes © 2020 The Authors a r X i v : . [ a s t r o - ph . C O ] S e p G. Aricò et al. – it is unclear whether they would be able to consistently modelbaryonic effects on the power spectrum and bispectrum.Exploring the predictions of baryonification for the bispec-trum is an important topic since it could highlight the pitfalls ofthe method or, instead, could support the correctness of the wholeapproach. Additionally, this comparison would represent an inde-pendent test of the method since higher-order statistics were neveremployed in the formulation of the baryonification algorithm.Motivated by these findings, in this paper we extend the analy-sis of the baryonification algorithm presented in Aricò et al. (2020)(hearafter A20), and here revisited, to the bispectrum. First, we ex-tend the model to account for gas that has been reaccreted by halos,and then show how different model parameters change the powerspectrum and bispectrum. Then, we show that our baryonificationimplementation can simultaneously reproduce, to better than ,the power spectrum and bispectrum measured in six state-of-the-arthydrodynamical simulations at k ≤ h Mpc − and at z ≤ . Wefurthermore explore the impact that the different components of themodel, e.g. central galaxy, ejected gas and back-reaction onto darkmatter, have on the matter bispectrum.This paper is structured as follow: in §2 we describe our nu-merical simulations, while in §3 we present our methodologies forbaryonic and cosmology modelling of the density field. In §4 wediscuss the impact of baryons in the bispectrum, whereas in §5we show our fits to the hydrodinamical simulations. We give ourconclusions in §6. In this work we use the same suite of N -body simulations used inA20. We refer the reader to it for further details, and here we onlyprovide a brief description.Our gravity-only simulations were carried out with l-gadget-3 (Angulo et al. 2012), a modified version of gadget (Springel2005). We employ simulations of box sizes: L=64, 128, and256 h − Mpc containing , , particles, respectively.We adopt the Nenya cosmology, as defined by Contreras et al.(2020): Ω cdm = 0 . , Ω b = 0 . , Ω Λ = 0 . , H = 60km s − Mpc − , n s = 1 . , σ = 0 . , τ = 0 . , (cid:80) m ν = 0 , w = − , and w a = 0 . These parameters are optimal to rescalethem to a large range of cosmologies (Contreras et al. 2020; Anguloet al. 2020). To test the accuracy of this rescaling, we will considertwo additional simulations of h − Mpc and particles: oneadopting the Nenya cosmology, and the other a massless neutrino
Planck cosmology (Planck Collaboration et al. 2018) .The initial conditions of all our simulations were computedwith the “fixed and paired” technique described in Angulo &Pontzen (2016), thus their cosmic variance is heavily suppressed.To compute the statistics of the density field, we use catalogue ofsimulation particles, selected homogeneously, diluted by a factor of . If not specified otherwise, our results will be computed with ourL=256 h − Mpc simulation, with which we expect our results to beconverged to about 2% for both bispectrum and power spectrum, asAppendix A shows. Ω cdm = 0 . , Ω b = 0 . , Ω Λ = 0 . , H = 67 . s − Mpc − , n s = 0 . , σ = 0 . , τ = 0 . , (cid:80) m ν = 0 , w = − , w a = 0 . Considering an overdensity field in Fourier space δ ( k ) , we definethe power spectrum as (cid:104) δ ( k ) δ ( k ) (cid:105) ≡ (2 π ) δ D ( k + k ) P ( k ) (1)and the bispectrum as (cid:104) δ ( k ) δ ( k ) δ ( k ) (cid:105) ≡ (2 π ) δ D ( k + k + k ) B ( k , k , k ) , (2)where (cid:104) ... (cid:105) denotes the ensemble average and δ D is the Dirac’s delta.To reduce the dependence of the bispectrum on the power spectrumand cosmology, we will mostly consider the reduced bispectrum (Scoccimarro 2000; Sefusatti & Komatsu 2007), defined as Q ( k , k , k ) ≡ B ( k , k , k ) P ( k ) P ( k ) + P ( k ) P ( k ) + P ( k ) P ( k ) . (3)We will mostly focus on the equilateral configuration, k = k = k , since it is expected to contain the most independentinformation from the power spectrum. In this case, Eq. 3 is reducedto: Q ( k ) = B ( k )3 P ( k ) . (4)We measure the bispectrum using bskit (Foreman et al.2019) , an extension of nbodykit (Hand et al. 2018) which usesa Fast Fourier Transform (FFT)-based bispectrum estimator (Scoc-cimarro 2000). Both the bispectrum and the power spectrum aremeasured in two interlaced grids (Sefusatti et al. 2016) employ-ing a triangular shaped cloud mass assignment scheme. The shotnoise contribution is estimated as / ¯ n for the power spectrum, andas / ¯ n + 1 / ¯ n [ P ( k ) + P ( k ) + P ( k )] for the bispectrum, andsubtracted. Finally, we have rebinned all the measurements in 25logarithmic bins over the interval [0 . , . h Mpc − . Additionally,when measuring the clustering on small scales, we use the “folding”technique (Jenkins et al. 1998; Colombi et al. 2009), described inAppendix B, which reduces CPU and memory usage.We will compare our results against the power spectra and bis-pectra from a number of cosmological hydrodynamical simulations,as measured by Foreman et al. (2019) . Specifically, we use fourstate-of-the-art hydrodynamical simulations: BAHAMAS (Mc-Carthy et al. 2017, 2018), EAGLE (Schaye et al. 2015; Crain et al.2015; McAlpine et al. 2016; Hellwing et al. 2016; The EAGLEteam 2017), Illustris (Vogelsberger et al. 2013, 2014a,b; Sijackiet al. 2015), and Illustris TNG-300 (Springel et al. 2018; Pillepichet al. 2018; Nelson et al. 2018; Naiman et al. 2018; Marinacci et al.2018; Nelson et al. 2019). In the case of BAHAMAS, we considertwo additional AGN feedback calibrations, dubbed as “low-AGN”and “high-AGN”: in the first one the temperature at which the AGNis activated is lower and thus the AGN feedback is weaker; whereasin the latter AGN feedback is stronger with respect to the standardrun. https://github.com/sjforeman/bskit https://github.com/sjforeman/hydro_bispectrum http://icc.dur.ac.uk/Eagle/ MNRAS000
Planck cosmology (Planck Collaboration et al. 2018) .The initial conditions of all our simulations were computedwith the “fixed and paired” technique described in Angulo &Pontzen (2016), thus their cosmic variance is heavily suppressed.To compute the statistics of the density field, we use catalogue ofsimulation particles, selected homogeneously, diluted by a factor of . If not specified otherwise, our results will be computed with ourL=256 h − Mpc simulation, with which we expect our results to beconverged to about 2% for both bispectrum and power spectrum, asAppendix A shows. Ω cdm = 0 . , Ω b = 0 . , Ω Λ = 0 . , H = 67 . s − Mpc − , n s = 0 . , σ = 0 . , τ = 0 . , (cid:80) m ν = 0 , w = − , w a = 0 . Considering an overdensity field in Fourier space δ ( k ) , we definethe power spectrum as (cid:104) δ ( k ) δ ( k ) (cid:105) ≡ (2 π ) δ D ( k + k ) P ( k ) (1)and the bispectrum as (cid:104) δ ( k ) δ ( k ) δ ( k ) (cid:105) ≡ (2 π ) δ D ( k + k + k ) B ( k , k , k ) , (2)where (cid:104) ... (cid:105) denotes the ensemble average and δ D is the Dirac’s delta.To reduce the dependence of the bispectrum on the power spectrumand cosmology, we will mostly consider the reduced bispectrum (Scoccimarro 2000; Sefusatti & Komatsu 2007), defined as Q ( k , k , k ) ≡ B ( k , k , k ) P ( k ) P ( k ) + P ( k ) P ( k ) + P ( k ) P ( k ) . (3)We will mostly focus on the equilateral configuration, k = k = k , since it is expected to contain the most independentinformation from the power spectrum. In this case, Eq. 3 is reducedto: Q ( k ) = B ( k )3 P ( k ) . (4)We measure the bispectrum using bskit (Foreman et al.2019) , an extension of nbodykit (Hand et al. 2018) which usesa Fast Fourier Transform (FFT)-based bispectrum estimator (Scoc-cimarro 2000). Both the bispectrum and the power spectrum aremeasured in two interlaced grids (Sefusatti et al. 2016) employ-ing a triangular shaped cloud mass assignment scheme. The shotnoise contribution is estimated as / ¯ n for the power spectrum, andas / ¯ n + 1 / ¯ n [ P ( k ) + P ( k ) + P ( k )] for the bispectrum, andsubtracted. Finally, we have rebinned all the measurements in 25logarithmic bins over the interval [0 . , . h Mpc − . Additionally,when measuring the clustering on small scales, we use the “folding”technique (Jenkins et al. 1998; Colombi et al. 2009), described inAppendix B, which reduces CPU and memory usage.We will compare our results against the power spectra and bis-pectra from a number of cosmological hydrodynamical simulations,as measured by Foreman et al. (2019) . Specifically, we use fourstate-of-the-art hydrodynamical simulations: BAHAMAS (Mc-Carthy et al. 2017, 2018), EAGLE (Schaye et al. 2015; Crain et al.2015; McAlpine et al. 2016; Hellwing et al. 2016; The EAGLEteam 2017), Illustris (Vogelsberger et al. 2013, 2014a,b; Sijackiet al. 2015), and Illustris TNG-300 (Springel et al. 2018; Pillepichet al. 2018; Nelson et al. 2018; Naiman et al. 2018; Marinacci et al.2018; Nelson et al. 2019). In the case of BAHAMAS, we considertwo additional AGN feedback calibrations, dubbed as “low-AGN”and “high-AGN”: in the first one the temperature at which the AGNis activated is lower and thus the AGN feedback is weaker; whereasin the latter AGN feedback is stronger with respect to the standardrun. https://github.com/sjforeman/bskit https://github.com/sjforeman/hydro_bispectrum http://icc.dur.ac.uk/Eagle/ MNRAS000 , 1–14 (2020) atter bispectrum with baryons Figure 1.
Baryonic effects at z = 0 on the power spectrum (left panel), equilateral bispectrum (central panel) and reduced bispectrum (right panel), measuredin 6 hydrodynamical simulations: BAHAMAS (standard, low and high AGN), EAGLE, Illustris and Ilustris TNG-300. We display the ratio of S = { P, B, Q } estimated in the full hydrodynamical simulation to that in their respective gravity-only counterpart. Figure 2.
Density profiles of gas (blue diamonds) and stars (orange circles) as measured in the Illustris TNG-300 simulation, at z = 0 . Each panel shows adifferent mass bin: . − h − M (cid:12) (left panel), − . h − M (cid:12) (central panel), and . − h − M (cid:12) (right panel). The baryonificationmodel that best fits simultaneously the three density profiles is shown as blue and orange shaded bands for gas and stars, respectively. The different gas andstellar subcomponents are displayed with different line styles, according to the legend. Note that the reaccreted gas density is consistent with zero, thus notappearing in the plot. In Fig. 1 we show the baryonic effects on the power spectrum,bispectrum, and reduced bispectrum. We display the ratio of theclustering measured in the full hydrodynamical simulations to thatin their gravity-only counterparts. Different colours show the resultsfor different simulations, as indicated by the legend.We see that the amplitude of baryonic effects considerablyvaries among simulations, in both power spectra and bispectra. Inparticular, Illustris and Bahamas high-AGN show the strongest sup-pression in both these statistics, likely due to their strong Supernovaeand AGN feedback. On the contrary, EAGLE and Illustris TNG-300show the smallest baryonic effects also likely related to their com- paratively weak feedback in massive halos. We highlight that bothof these simulations display an enhancement of the bispectrum at k ≈ − h Mpc − , which has been linked to the presence oflate-time reaccreted gas by Foreman et al. (2019). We will test thishypothesis with our baryonification framework later on.Interestingly, whereas baryons can either suppress or enhancethe gravity-only bispectrum, they appear simpler in the reducedbispectrum: baryons always enhance Q ( k ) on small scales. Quali-tatively, there seems to be a clear correlation between the baryoniceffects in the power spectrum and bispectrum. However, this cor-relation is not perfect: Illustris and BAHAMAS high-AGN show MNRAS , 1–14 (2020)
G. Aricò et al.
Figure 3.
Accuracy of the cosmology-rescaling algorithm when used witha baryonification procedure. [
Upper panel: ] Ratio of the baryonified andgravity-only mass power spectra and reduced bispectra at z = 0 . Symbolsshow the results using a simulation adopting the Planck cosmology, whereaslines indicate the results using a simulation rescaled to the same
Planck cos-mology. We provide results for 4 baryonification models roughly consistentwith the effects expected in BAHAMAS (red), EAGLE (brown), Illustris(green) and Illustris TNG-300 (blue) models. [
Lower panel: ] Differencebetween the baryonic effects measured in the target and scaled simulationshown in the upper panels. The grey shaded band marks a discrepancy of3%. similar effects on the reduced bispectrum, but the effects on thepower spectrum are clearly different . In the next sections we willexplore whether baryonification methods can successfully describeall these features at high precision. Given a particle field in a gravity-only (GrO) N-body simulation,we can obtain the mass field in arbitrary cosmologies and baryonicscenarios by manipulating the positions and masses of the particles.To do so, we use the framework described in A20, which we recapnext.We first apply a “cosmology rescaling” to obtain a simulationat a desired cosmological parameter set (Angulo & White 2010). Forthis, we scale the lengths, masses and time (by selecting differentsnapshots) of our simulation in order to match the amplitude ofthe linear density fluctuation of another cosmology. This techniquehas been extensively tested (Ruiz et al. 2011; Renneby et al. 2018;Mead & Peacock 2014a,b; Mead et al. 2015; Angulo & Hilbert However, note that Illustris simulates a box less than h − Mpc , thustheir results could be affected by cosmic variance and lack of long wave-modes. As showed in Appendix A, massive haloes contribute to baryoniceffects more in the bispectrum than in the power spectrum. As a consequence,the reduced bispectrum measured in relatively small boxes is suppressed atsmall scales with respect to larger boxes.
Figure 4.
Modifications to the matter power spectrum (upper panels), andreduced bispectrum (lower panels) at z = 0 caused by baryons, as predictedby our baryonification algorithm with parameters mimicking the effectsexpected in the Illustris TNG-300 simulations. The total baryonic effect isdecomposed into the contribution of each component, namely ejected gas,galaxies, hot bound gas, reaccreted gas, and dark matter, according to thelegend. < accuracy in the matter power spectrum and in the matterbispectrum up to k ∼ h − Mpc (Contreras et al. 2020, Zennaro etal. in prep), over a broad range of cosmologies, even beyond- Λ CDM.Note we expect a higher accuracy for the ratio of baryonified overgravity-only outputs, as we will show later.We then apply a “baryonification” algorithm to further dis-place the particles of the simulation, and mimic the effect of differ-ent baryonic components. In A20, each halo was assumed to havefour components: dark matter, a central galaxy, bound gas, and ex-pelled gas. In this work, we additionally model satellite galaxiesand late-time reaccreted gas, which we describe in detail in thenext subsection. The density profiles of all these components areparametrised with physically motivated functional forms, whereasthe GrO halo density profile is modelled with as a NFW profile(Navarro et al. 1997). Once we have the initial and the “baryonic”density profiles, we compute a displacement field which, applied tothe halo particles, distorts their distributions accordingly.
One of the main advantages of the baryon correction model is itsextreme flexibility, which allows us to make modifications or includenew physics according to various possible scenarios. In this work,we have implemented in the model of A20 the following four mainupdates: • The adoption of a more flexible functional form for the boundgas; • The inner slope of the power-law in the central galaxy is a newfree parameter; • The modelling of a satellite galaxies component;
MNRAS000
MNRAS000 , 1–14 (2020) atter bispectrum with baryons • The inclusion of a late-time reaccreted gas component;We find that the parametrisation of the bound gas density shapeused in Schneider & Teyssier (2015); Aricò et al. (2020) is notflexible enough to match the profiles measured in a wide range ofhalo masses of hydrodynamical simulations. We therefore use herea more flexible shape, with an explicit dependence on the halo mass.The shape of the bound gas now reads: ρ BG ( r ) = y (1 + r/r inn ) β i r/r out ) ) (5)where y is a normalisation factor, obtained by imposing (cid:82) r d r πr ρ BG ( r ) = f BG M . The profile is a double power-law with two characteristic scales, r inn and r out , defining where theslope changes at small and large radii, respectively. We define theinner radius r inn = θ inn × r and r out = θ out × r , with θ inn and θ out being free parameters of the model. The gas inner slopeexplicitly depends on halo mass as β i = 3 − ( M inn /M ) µ i ,with the characteristic mass M inn and µ i as free parameters. Afterchecking the small impact that µ i has on both power spectrum andbispectrum, we have fixed its value to µ i = 0 . , in agreement tothe Model A-avrg in Schneider et al. (2019).This profile is similar to that in Schneider et al. (2019), withthe main difference being that in our model the bound gas perfectlytraces the dark matter on scales beyond r out and the ejected gasdecay exponentially, whereas Schneider et al. (2019) models a singlegas component, with the slope at large radii as a free parameter.The central galaxy density profile is given by ρ CG ( r ) = y R h r α g exp (cid:34) − (cid:18) r R h (cid:19) (cid:35) , (6)where y is found imposing (cid:82) r d ρ CG (r) = f CG M , the half-mass radius is R h = 0 . × r and α g , the inner slope of thecentral galaxy, is a free parameter of the model, with a fiducial value α g = 2 .In addition to the central galaxy, we add the stellar componentof satellite galaxies. Stellar mass is, to a good approximation, col-lisionless and thus a good tracer of dark matter. For this reason wemodel the contribution of satellite galaxies as the dark matter. Thedark matter back-reacts to the baryonic potential well, and thereforethe satellite galaxies, being a linearly biased tracer of the dark mat-ter, are quasi-adiabatically relaxed. We refer the reader to A20 forthe details of the implementation of the back-reaction mechanism.Motivated by the hypothesis of Foreman et al. (2019), who sug-gested the presence of an overdensity of gas reaccreted at late timesto explain the maximum in the bispectrum around k ≈ . h Mpc − in the Illustris TNG-300 simulation, we added to our model a newgas component, which mimics such gas overdensity. We assume thisnew component to be Gaussian shaped: ρ RG ( r ) = y √ πσ r exp (cid:20) − ( r − µ r ) (2 σ r ) (cid:21) , (7)where y = f RG M / (cid:82) r πr ρ RG ( r ) dr .For simplicity, we assume the gas overdensity to have a fixedspatial distribution in terms of the halo virial radius, µ r = 0 . × r and σ r = 0 . × r , and after checking that our main results arenot affected by this choice. We let free instead the mass fraction, f RG , as explained in what follows.All the density profiles of the baryon correction model are Figure 5.
Upper panel:
Baryon suppression of the matter power spectrumat z = 0 , when applying the baryon correction model to haloes smaller than h − M (cid:12) (green dotted line), between − h − M (cid:12) (brownsolid line), − h − M (cid:12) (golden dashed line) and to all the haloes(blue dots) of our h − Mpc simulation.
Lower panel:
Same as theupper panel, but for the ratios between baryonic and gravity-only results inthe reduced bispectrum. normalised to M , with the abundance of each component deter-mined by its respective mass fraction. The dark matter fraction isfixed by cosmology, f DM = 1 − Ω b / Ω m .The central galaxy fraction is given by an abundance-matchingparametrisation (Behroozi et al. 2013): f CG ( M ) = (cid:15) (cid:18) M M (cid:19) g (log ( M /M )) − g (0) , (8) g ( x ) = − log (10 αx + 1) + δ (log (1 + exp( x ))) γ − x ) . (9)We use the best-fitting parameters at z = 0 given by Kravtsovet al. (2018), along with the redshift dependence given by Behrooziet al. (2013), both reported in Appendix A of A20 and not includedhere for the sake of brevity.Satellite and central mass fractions have the same parametricform, and their parameters are assumed to be linearly dependente.g. M , sat ( z = 0) = α sat M , cen ( z = 0) , with α sat as a freeparameter of the model, similar to the approach of Watson & Conroy(2013).The halo gas mass fraction, defined as the sum of the boundgas and the reaccreted gas, is f HG ( M ) = f BG + f RG = Ω b / Ω m − f CG − f SG M c /M ) β , (10) MNRAS , 1–14 (2020)
G. Aricò et al. with M c and β free parameters, and f CG , f SG , the central andsatellite galaxy mass fractions, respectively. The reaccreted gas massfraction is f RG ( M ) = Ω b / Ω m − f CG − f SG − f HG M r /M ) β r == f HG ( M c /M ) β M r /M ) β r , (11)with M r as a free parameter and β r fixed for simplicity to β r = 2 .Finally, the bound and ejected gas mass fractions are set bymass conservation: f BG = f HG − f RG ; (12) f EG = Ω b / Ω m − f CG − f SG − f HG . (13)As an example, we show in Fig. 2 how this updated modelis able to reproduce at the same time the gas and stellar den-sity profiles measured in three different halo mass bins of Illus-tris TNG-300, [10 . − ] h − M (cid:12) , [10 , . ] h − M (cid:12) and [10 . , ] h − M (cid:12) . Note that here our reaccreted mass fractionsare consistent with zero, thus not appearing in the plot. In A20 we showed that applying a baryonification algorithm togetherwith a cosmology-rescaled simulation leaded to percent-accurateresults in the power spectrum. We now perform an analogous testto validate the performance of the updated model and extend theanalysis to the bispectrum.In Fig. 3 we compare the baryonic effects on the power spec-trum and reduced bispectrum as measured in a simulation carriedout with a
Planck cosmology and a simulation carried our with a
Nenya cosmology and then rescaled to a
Planck cosmology (c.f.§2). These two cases are denoted as target and scaled, respectively,and displayed by symbols and lines as indicated by the legend.We display 4 different baryonification parameter sets, chosento roughly reproduce the clustering of EAGLE, Illustris, IllustrisTNG-300 and BAHAMAS. We can see that the difference betweenapplying the BCM on top a rescaled or target simulation is less than in the power spectrum and less than in the reduced bispec-trum. We show these results only for z = 0 , but we have explicitlychecked that at higher redshifts we obtain similar outcomes.We note that the initial conditions of the target simulation werenot set to match that of the simulation we scale, nor its volume havebeen chosen to match the volume of the rescaled simulation (whichcould have increase the agreement further). Nevertheless, the errorswe obtain are comparable to our target accuracies for reproducingthe baryonic effects on the power spectrum and bispectrum. In this section we systematically explore the effects that the variousbaryonic components, and the free parameters associated to them,produce on the clustering.We first isolate the effect of each baryonic component by se-lecting them one-by-one and considering all the others collisionless(thus behaving like dark matter). As shown in Fig. 4, we find that, in agreement with Schneider & Teyssier (2015); Aricò et al. (2020),the ejected gas largely dominates the suppression in the power spec-trum, despite its low mass fraction. Interestingly, the ejected gasshows the largest effect also in the reduced bispectrum, but in thiscase it causes an enhancement of the power at all scales.As an qualitative explanation, let us consider two overdensityfields, δ BCM and δ GrO . Assuming that one is suppressed with re-spect to the other, δ BCM = (1 − α ) δ GrO , it is easy to show thatthe ratios between their power spectra and equilateral bispectra are P BCM /P GrO = (1 − α ) , and B BCM /B GrO = (1 − α ) , respec-tively. Therefore, the reduced bispectrum ratio is Q BCM /Q GrO =(1 − α ) − . In other words, we observe an enhancement of thereduced bispectrum because the suppression in the bispectrumis smaller than the squared suppression of the power spectrum.The other components are, in this particular setting of the BCMwhich roughly mimics the BAHAMAS simulation, subdominant,contributing to about in the power spectrum and reduced bis-pectrum. The reaccreted gas, in particular, causes an enhancementat small scales in both the power spectrum and reduced bispectrum.It is interesting to explore which halo masses contribute themost to the baryonic effects on clustering. In order to answer thisquestion, we have split the halo catalogue of our simulation in differ-ent mass bins, and then we have applied our BCM separately to eachof them. In Fig. 5 we show how haloes between − h − M (cid:12) contribute more than of the effect on the power spectrum atsmall scales. Haloes of − h − M (cid:12) are dominant at largescales in the power spectrum, whereas at the small scales, slightlysmaller haloes contribute more. Haloes with M < h − M (cid:12) and M > h − M (cid:12) contribute for less than percent, andonly at small scales.The relative contribution of halos of different mass slightlychanges in the case of the reduced bispectrum. The dominantcontribution is still from haloes of − h − M (cid:12) , but therelative impact of the most massive haloes in the simulation( M > h − M (cid:12) ) is not as small as for the power spectrum.The fact that the bispectrum is more sensitive to the largest haloes isnot surprising (see e.g. Foreman et al. 2019), and has as a practicaloutcome the slower convergence of the bispectrum with simulatedvolume compared to that of power spectrum, which we investigatein Appendix A.We quantify now the impact of the free parameters in the powerspectrum and reduced bispectrum. To do so, we vary each parameterone by one, while keeping the others fixed to the value that best fitsthe BAHAMAS (standard AGN) simulation, described in §5. Theintervals in which we vary parameters, in log , are the following: M c ∈ [9 , , η ∈ [ − . , . , β ∈ [ − , . , M r ∈ [12 , , M , z0 , cen ∈ [9 , , θ inn ∈ [ − , − . , θ out ∈ [ − . , , M inn ∈ [12 , .Note we do not show any free parameters for satellite galaxies,as they have a negligible impact on the matter clustering. In fact,they are a biased tracer of the dark matter, and additionally theirmass fraction is very small. Given that the relaxation of the darkmatter contributes only for a few percent in the matter clustering, itis easy to see why the baryonic effect caused by satellite galaxies isnegligible. Thus, we fix their values to the best-fitting of the stellarprofile of the Illustris TNG-300 simulation found in §3.1.In Fig. 6 we display the mass power spectra, bispectra, and re-duced bispectra obtained after applying the BCM to a GrO N -bodysimulation. Each panel varies a single parameter of the model whilekeeping the others fixed to their fiducial value. Bluer (redder) col-ors represent low (high) parameter values. We can see that almostall the parameter combinations predict a suppression in the power MNRAS000
Planck cosmology (c.f.§2). These two cases are denoted as target and scaled, respectively,and displayed by symbols and lines as indicated by the legend.We display 4 different baryonification parameter sets, chosento roughly reproduce the clustering of EAGLE, Illustris, IllustrisTNG-300 and BAHAMAS. We can see that the difference betweenapplying the BCM on top a rescaled or target simulation is less than in the power spectrum and less than in the reduced bispec-trum. We show these results only for z = 0 , but we have explicitlychecked that at higher redshifts we obtain similar outcomes.We note that the initial conditions of the target simulation werenot set to match that of the simulation we scale, nor its volume havebeen chosen to match the volume of the rescaled simulation (whichcould have increase the agreement further). Nevertheless, the errorswe obtain are comparable to our target accuracies for reproducingthe baryonic effects on the power spectrum and bispectrum. In this section we systematically explore the effects that the variousbaryonic components, and the free parameters associated to them,produce on the clustering.We first isolate the effect of each baryonic component by se-lecting them one-by-one and considering all the others collisionless(thus behaving like dark matter). As shown in Fig. 4, we find that, in agreement with Schneider & Teyssier (2015); Aricò et al. (2020),the ejected gas largely dominates the suppression in the power spec-trum, despite its low mass fraction. Interestingly, the ejected gasshows the largest effect also in the reduced bispectrum, but in thiscase it causes an enhancement of the power at all scales.As an qualitative explanation, let us consider two overdensityfields, δ BCM and δ GrO . Assuming that one is suppressed with re-spect to the other, δ BCM = (1 − α ) δ GrO , it is easy to show thatthe ratios between their power spectra and equilateral bispectra are P BCM /P GrO = (1 − α ) , and B BCM /B GrO = (1 − α ) , respec-tively. Therefore, the reduced bispectrum ratio is Q BCM /Q GrO =(1 − α ) − . In other words, we observe an enhancement of thereduced bispectrum because the suppression in the bispectrumis smaller than the squared suppression of the power spectrum.The other components are, in this particular setting of the BCMwhich roughly mimics the BAHAMAS simulation, subdominant,contributing to about in the power spectrum and reduced bis-pectrum. The reaccreted gas, in particular, causes an enhancementat small scales in both the power spectrum and reduced bispectrum.It is interesting to explore which halo masses contribute themost to the baryonic effects on clustering. In order to answer thisquestion, we have split the halo catalogue of our simulation in differ-ent mass bins, and then we have applied our BCM separately to eachof them. In Fig. 5 we show how haloes between − h − M (cid:12) contribute more than of the effect on the power spectrum atsmall scales. Haloes of − h − M (cid:12) are dominant at largescales in the power spectrum, whereas at the small scales, slightlysmaller haloes contribute more. Haloes with M < h − M (cid:12) and M > h − M (cid:12) contribute for less than percent, andonly at small scales.The relative contribution of halos of different mass slightlychanges in the case of the reduced bispectrum. The dominantcontribution is still from haloes of − h − M (cid:12) , but therelative impact of the most massive haloes in the simulation( M > h − M (cid:12) ) is not as small as for the power spectrum.The fact that the bispectrum is more sensitive to the largest haloes isnot surprising (see e.g. Foreman et al. 2019), and has as a practicaloutcome the slower convergence of the bispectrum with simulatedvolume compared to that of power spectrum, which we investigatein Appendix A.We quantify now the impact of the free parameters in the powerspectrum and reduced bispectrum. To do so, we vary each parameterone by one, while keeping the others fixed to the value that best fitsthe BAHAMAS (standard AGN) simulation, described in §5. Theintervals in which we vary parameters, in log , are the following: M c ∈ [9 , , η ∈ [ − . , . , β ∈ [ − , . , M r ∈ [12 , , M , z0 , cen ∈ [9 , , θ inn ∈ [ − , − . , θ out ∈ [ − . , , M inn ∈ [12 , .Note we do not show any free parameters for satellite galaxies,as they have a negligible impact on the matter clustering. In fact,they are a biased tracer of the dark matter, and additionally theirmass fraction is very small. Given that the relaxation of the darkmatter contributes only for a few percent in the matter clustering, itis easy to see why the baryonic effect caused by satellite galaxies isnegligible. Thus, we fix their values to the best-fitting of the stellarprofile of the Illustris TNG-300 simulation found in §3.1.In Fig. 6 we display the mass power spectra, bispectra, and re-duced bispectra obtained after applying the BCM to a GrO N -bodysimulation. Each panel varies a single parameter of the model whilekeeping the others fixed to their fiducial value. Bluer (redder) col-ors represent low (high) parameter values. We can see that almostall the parameter combinations predict a suppression in the power MNRAS000 , 1–14 (2020) atter bispectrum with baryons Figure 6.
Modifications to the matter power spectrum (upper panels), bispectrum (central panels) and reduced bispectrum (lower panels) at z = 0 caused bybaryons, according to our baryonification algorithm. Each column varies one of the free parameters of the model while keeping the others fixed at their fiducialvalue. Parameter ranges are log M c ∈ [9 , , log η ∈ [ − . , . , log β ∈ [ − , . , log M r ∈ [12 , , log M , z0 , cen ∈ [9 , , log θ inn ∈ [ − , − . , log θ out ∈ [ − . , , log M inn ∈ [12 , . Blue to red colors denote low to high parameter values. spectrum and an enhancement on the reduced bispectrum, at all thescales. Specifically, by increasing η (the parameter which set themaximum range of the AGN feedback), the suppression (enhance-ment) of the power spectrum (reduced bispectrum) is pushed, asexpected, towards larger scales. The parameters M c and β set thefraction of gas which is retained in haloes of a given mass, and thusalso the mass of gas that is expelled. Therefore, is not surprisingthat these parameters have a big impact on both power spectrum andbispectrum, given that the ejected gas component is the dominantone. Varying M c we span a range in the clustering; in par-ticular, higher values mean that increasingly larger haloes are freeof gas, thus more ejected gas. In these cases we see, accordingly, alarger suppression in the power spectrum and enhancement in thereduced bispectrum.Varying the shape of the bound gas through the parameters θ inn , θ out and M inn has an impact only on small scales. Specifically,the model seems very sensitive to θ inn , for which we see a substantialenhancement of both power spectrum and reduced bispectrum whenchanging the inner gas slope at increasingly smaller radii. On theother hand, the dependence on M inn looks negligible. As expected,increasing M , z0 , cen , and thus having the peak of the star formationat higher halo masses, results in more power at small scales.Finally, we see that the impact of the late-time reaccreted gas isvery modest, despite we vary its mass fraction from practically zeroto a limit value of ≈ for some halo masses. Arguably, the effect of this parameter in both power spectrum and reduced bispectrum,can be absorbed by a combination of the other parameters of themodel, but might become more important on smaller scales.We note that the models shown in Fig. 6 are just an illustrativeexample, and do not encompass all the possible modifications givenby the BCM: even if the described trends would be likely similar,changing the underlying fiducial model would result in differentamplitude and shapes of baryonic effects. In this section, we explore whether the BCM is able to reproduce theimpact of baryons in six different state-of-the-art hydrodynamicalsimulations, namely EAGLE, Illustris, Illustris TNG-300, and threedifferent AGN implementations of BAHAMAS. We remind thereader that these simulations differ in cosmology, N -body code,sub-grid physics, box size, and observables with which they havebeen calibrated. They show a difference of at z = 0 in thepower spectrum and in the reduced bispectrum, thus being agood benchmark for the flexibility and realism of our model.We fit the power spectrum and the reduced bispectrum of eachhydrodynamical simulation over the range . < k/ ( h Mpc − ) < , both separately and jointly, varying seven free parameters: M c , η , β , M , z0 , cen , θ inn , θ out , M inn within the priors shown in §4. Weassume no correlation among power spectrum and bispectrum nor MNRAS , 1–14 (2020)
G. Aricò et al.
Figure 7.
Measurements of the baryonic impact to the matter power spectrum, S ( k ) ≡ P/P
GrO (upper panels), equilateral bispectrum, S ( k ) ≡ B/B
GrO (central panels), and reduced equilateral bispectrum, S ( k ) ≡ Q/Q
GrO (bottom panels), in different hydrodynamical simulations according to the legend(symbols), at z = 0 (left), z = 1 (centre) and z = 2 (right). The best-fitting baryonification model constrained using only the power spectrum is displayed asdashed lines, whereas the best-fitting model constrained on both the power spectrum and reduced bispectrum shown with solid lines. Grey vertical dotted linesmark the scales where the estimated shotnoise contributes to > / of the clustering amplitude. Figure 8.
Best-fitting model to Illustris TNG-300 power spectrum and bispectrum, at z = 0 . The contribution of the varyous baryon component are isolatedin the power spectrum (left), equilaetral bispectrum (centre) and reduced bispectrum (right). Note that the bound gas, and the back-reaction to the dark matter,are the principale causes of the bump visible in the bispectrum at small scales. MNRAS000
Best-fitting model to Illustris TNG-300 power spectrum and bispectrum, at z = 0 . The contribution of the varyous baryon component are isolatedin the power spectrum (left), equilaetral bispectrum (centre) and reduced bispectrum (right). Note that the bound gas, and the back-reaction to the dark matter,are the principale causes of the bump visible in the bispectrum at small scales. MNRAS000 , 1–14 (2020) atter bispectrum with baryons Figure 9.
Impact of the redshifts evolution in the baryonic correction param-eters. We have fitted the power spectra (upper panels) and reduced bispectra(lower panels) of the hydrodynamical simulations reported in the legendat z = 0 , and then applied the same model at higher redshifts z = 1 , ,assuming our best-fitting parameters to be redshift independent. among the measurements at different wavenumbers. Specifically,we use an empirical approach similarly to A20, where the covari-ance matrix is directly estimated by the intra-data variance, givingthe same weights to power spectrum and bispectrum. We expectthe errors associated to the bispectrum ratios to be larger than theone associated to the power spectrum (see for instance the errorsmeasured by Foreman et al. (2019) by dividing the hydrodynamicalsimulation volume in subboxes). Nevertheless, being the purposeof this test to asses the accuracy of the joint fit of power spectrumand bispectrum, we avoid to give more weight to the former to notdegrade the fit of the latter.To perform the fit, we have implemented a particle swarm op-timisation algorithm (Kennedy & Eberhart 1995). In this algorithm,a pack of particles efficiently searches the minimum of a functionin a given parameter space. Each particle communicates with theothers at each step, and they are attracted both to their local and theswarm global minima, with a relative strength that can be tuned.The velocity and position of the particles are updated in every step,depending solely on the swarm status in the previous step. For ourapplication, we use a swarm of particles and iterations, find-ing that an average of 100-150 steps are enough to converge to theglobal minimum.In Fig. 7 we present the main result of this paper. We showthe best BCM fits at three different redshifts, z = 0 , , and . Wehave marked with a grey dotted line the scales where we estimatethe shotnoise amplitude to be approximately 30% of the clusteringamplitude: k ≈ h Mpc − at z = 2 , and k ≈ h Mpc − at z = 1 .We remind the reader that, in the cosmology rescaling process, thebox of the simulations can vary of length, and so the shotnoiselevel can be slightly different. Due to the significant contributionof shotnoise, results at small scales and high redshifts should beinterpreted carefully.Dashed lines show the results when fitting only the powerspectrum measurements. In this case, we recover the accuracy of found in A20 in the power spectrum at all scales and redshifts.However, the baryonic impact on the bispectrum can be over- orunder-estimated by up to 20%. In contrast, when fitting the power spectrum and bispectrum together , the accuracy of the power spec-trum slightly degrades, but we obtain significantly better agreementwith the bispectrum.For all redshifts and hydrodynamical simulations considered,we obtain joint fits that are − accurate for the power spectrumand for the bispectrum. We note that the worse performance isobtained for the simulations with the most extreme feedback e.g.Illustris and BAHAMAS high-AGN. In the case of BAHAMAS,which is arguably the most realistic simulation for our purposes,the fits describe simultaneously the baryonic effects on the powerspectrum and bispectrum to better than .These results are achieved considering the late-time reaccretedgas fixed to zero. In particular, we note that the bump around k ≈ − h Mpc − in the bispectrum of Illustris TNG-300 and EAGLEare correctly reproduced by the model, despite the absence of thereaccreted gas component. To understand which BCM componentcauses this enhancement of the bispectrum at small scales, we isolatethe impact of each baryon component to the clustering, similarly towhat done in § 4. This analysis, reported in Fig. 8, clearly showthat the bound gas causes an enhancement at small scales in bothpower spectrum and bispectrum. Furthermore, the back-reaction ofthe gas overdensity to the dark matter adds the necessary power toreproduce correctly the measurements of the Illustris TNG-300. Wecan conclude that, by simply assuming the gas as a double power-law, the model has enough flexibility (over the range of scales weconsidered) to explain the “bump” in the bispectrum measured inIllustris TNG-300 and EAGLE.We have repeated the fits letting free the corresponding massfraction, M r , however, this did not result in noticeably improvedfits. This finding is consistent with the hypothesis that, within theaccuracy of our model and simulated data, and over the scalesconsidered, the reaccreted gas is not necessary to reproduce theclustering of the hydrodynamical simulations analysed.Finally, one could wonder what is the smallest number of freeparameters necessary to produce accurate results. In A20 it is shownthat with only 4 parameters it is possible to fit the power spectrumat , and arguably the 7-8 parameters used here are degenerate,and effectively recastable into a model with a smaller parameter set.We leave the exploration of degeneracies between parameters andthe finding of a minimal-model parameter set for a future work. As already pointed out in previous works (Chisari et al. 2019; Aricòet al. 2020), despite the baryonic parameters do not have a specificredshift dependence, when fitting the clustering at different redshiftsthey show a clear evolution. To quantify the inaccuracies obtainedby fixing the baryonic parameters, we apply the baryonificationalgorithm to our simulation to snapshots that correspond to highredshifts ( z = 1 , ), using the best-fitting parameters at z = 0 found in § 5.In Fig. 9 we show the power spectra and reduced bispectraobtained. The error in the power spectrum is in most of the casesbelow , for extreme models around 7-10 % . On the other hand,the reduced bispectrum can be off of − . These errors mustbe taken into account when fixing at face value a set of baryonicparameters in multiple redshifts. We have so far analysed, for simplicity, only the equilateral con-figuration of the reduced bispectrum. In this section, we explore
MNRAS , 1–14 (2020) G. Aricò et al.
Figure 10.
Upper panel:
Measurements of the reduced squeezed bispectraat z = 0 in EAGLE, Illustris, Illustris TNG-300, BAHAMAS, Bahamaslow-AGN and BAHAMAS high-AGN (symbols), according to the legend.The shaded band show our prediction obtained by fitting power spectrumand reduced equilateral bispectrum inthe two most extreme models, EAGLEand BAHAMAS high-AGN. Lower panel:
Difference between the ratios ofreduced squeezed bispectra predicted by our baryon correction model andmeasured in the hydrodynamical simulations. The dashed lines show themodel parameters fitting only the power spectra of the hydrodynamical sim-ulations, whereas the solid lines show the model with parameters constrainedby fitting both power spectra and equilateral reduced bispectrum.
Figure 11.
Upper panel : Relation between the baryonic impact on the powerspectrum at k = 1 h Mpc − , defined as ∆ P ( k ) /P ( k ) , and the halo baryonfraction in haloes with a mass [ × , × ] M (cid:12) . The black dashedline displays the fit provided by van Daalen et al. (2019), with the grey andlight grey shaded bands marking a 1% and 2% deviation, respectively. Forcomparison, the colored stars emply the baryon fraction measured directly inthe hydrodynamical simulations by van Daalen et al. (2019). The colouredsymbols indicate the measurements of our baryonified simulations whenusing the best-fitting values calibrated against the power spectra (circles),and both power spectra and bispectra (diamonds). Lower panel:
Same as theupper panel, but for the reduced bispectrum. In this case, the dashed linerepresents a simple linear regression of the symbols displayed. the baryonification performance for the “squeezed” configuration,which measures the correlation between points on isosceles trian-gles with one side much smaller than the other two in k -space, sothat k (cid:28) k = k . The squeezed bispectrum might be seen asa “conditional” two-point correlation which quantifies the depen- dence of small-scale nonlinearities on the large-scale backgroundoverdensity.It has been shown that, in some cases, the baryonic effecton the squeezed bispectrum can be directly related to the powerspectrum at small scales, when considering a k long enough tonot be affected by baryonic physics (Barreira et al. 2019; Fore-man et al. 2019). Specifically, Barreira et al. (2019) have measuredthe “power spectrum response functions” in the Illustris TNG-300,using the separate universe approach, finding that they are largelyunaffected by baryonic physics. This suggests that the information inthe squeezed bispectrum is already contained in the power spectrum,and thus knowing the latter we can predict the former. However, asshown in Foreman et al. (2019), the analytical predictions givenfrom the power spectrum response function are not always in agree-ment with the hydrodynamical simulations, e.g. BAHAMAS. Thiscould be a hint that, in some cases, the response function are notfully specified by the power spectra.Here, we take a somewhat agnostic approach, and test if wecan predict correctly the squeezed bispectrum starting from theinformation contained in the power spectrum and the equilateralbispectrum. To do so, we apply to our gravity-only simulations aBCM with the parameters that reproduce both the power spectrumand reduced equilateral bispectrum for a given hydrodynamicalsimulations. Then, we measure the reduced squeezed bispectrum( k = k > k ∼ . h − Mpc ) and compare it with those mea-sured directly in the hydrodynamical simulations.In Fig. 10 we show the results obtained at z = 0 . First, wecan notice that, as for the case of equilateral configurations, whenconsidering baryon physics the reduced squeezed bispectrum is en-hanced with respect to the gravity-only one. However, the baryoniceffects in the squeezed bispectrum are smaller than those in theequilateral configuration – spanning a range of − , against a − measured in the equilateral configuration.We also see that our predictions for the squeezed reduced bis-pectrum agree very well with the simulation measurements, reach-ing a ∼ accuracy in all cases. This further supports the ideathat the modifications to the density field in the barionification isaccurately capturing the three-dimensional distortions induced bybaryons, and not simply fitting an effective distortion in the powerspectrum.For comparison, we also display in Fig. 10 the predictionswhen tuning our model using only the power spectrum. As for theequilateral bispectrum, the impact of baryons is not captured veryaccurately, with discrepancies generally within (EAGLE andBAHAMAS low-AGN ≤ ) to up to (Bahamas high-AGN). Recently, it has been shown that there is a tight correlation betweenbaryonic effects on the power spectrum and the baryonic fractioninside haloes of M ≈ M (cid:12) (van Daalen et al. 2019). The bestfits of the baryon correction model has been shown to be able toaccurately recover such correlation, even if for large power spec-trum suppression, which correspond to very strong AGN feedback,it tends to underestimate the baryon fraction measured in hydrody-namical simulation (Aricò et al. 2020).We now explore whether adding the information on the bis-pectrum the baryonic halo fractions become more constrained, andadditionally, whether there is a relation between reduced bispec-trum and baryon fraction, analogous to the one found for the powerspectrum.In Fig. 11 we show how, indeed, in the case of the BAHAMAS MNRAS000
Same as theupper panel, but for the reduced bispectrum. In this case, the dashed linerepresents a simple linear regression of the symbols displayed. the baryonification performance for the “squeezed” configuration,which measures the correlation between points on isosceles trian-gles with one side much smaller than the other two in k -space, sothat k (cid:28) k = k . The squeezed bispectrum might be seen asa “conditional” two-point correlation which quantifies the depen- dence of small-scale nonlinearities on the large-scale backgroundoverdensity.It has been shown that, in some cases, the baryonic effecton the squeezed bispectrum can be directly related to the powerspectrum at small scales, when considering a k long enough tonot be affected by baryonic physics (Barreira et al. 2019; Fore-man et al. 2019). Specifically, Barreira et al. (2019) have measuredthe “power spectrum response functions” in the Illustris TNG-300,using the separate universe approach, finding that they are largelyunaffected by baryonic physics. This suggests that the information inthe squeezed bispectrum is already contained in the power spectrum,and thus knowing the latter we can predict the former. However, asshown in Foreman et al. (2019), the analytical predictions givenfrom the power spectrum response function are not always in agree-ment with the hydrodynamical simulations, e.g. BAHAMAS. Thiscould be a hint that, in some cases, the response function are notfully specified by the power spectra.Here, we take a somewhat agnostic approach, and test if wecan predict correctly the squeezed bispectrum starting from theinformation contained in the power spectrum and the equilateralbispectrum. To do so, we apply to our gravity-only simulations aBCM with the parameters that reproduce both the power spectrumand reduced equilateral bispectrum for a given hydrodynamicalsimulations. Then, we measure the reduced squeezed bispectrum( k = k > k ∼ . h − Mpc ) and compare it with those mea-sured directly in the hydrodynamical simulations.In Fig. 10 we show the results obtained at z = 0 . First, wecan notice that, as for the case of equilateral configurations, whenconsidering baryon physics the reduced squeezed bispectrum is en-hanced with respect to the gravity-only one. However, the baryoniceffects in the squeezed bispectrum are smaller than those in theequilateral configuration – spanning a range of − , against a − measured in the equilateral configuration.We also see that our predictions for the squeezed reduced bis-pectrum agree very well with the simulation measurements, reach-ing a ∼ accuracy in all cases. This further supports the ideathat the modifications to the density field in the barionification isaccurately capturing the three-dimensional distortions induced bybaryons, and not simply fitting an effective distortion in the powerspectrum.For comparison, we also display in Fig. 10 the predictionswhen tuning our model using only the power spectrum. As for theequilateral bispectrum, the impact of baryons is not captured veryaccurately, with discrepancies generally within (EAGLE andBAHAMAS low-AGN ≤ ) to up to (Bahamas high-AGN). Recently, it has been shown that there is a tight correlation betweenbaryonic effects on the power spectrum and the baryonic fractioninside haloes of M ≈ M (cid:12) (van Daalen et al. 2019). The bestfits of the baryon correction model has been shown to be able toaccurately recover such correlation, even if for large power spec-trum suppression, which correspond to very strong AGN feedback,it tends to underestimate the baryon fraction measured in hydrody-namical simulation (Aricò et al. 2020).We now explore whether adding the information on the bis-pectrum the baryonic halo fractions become more constrained, andadditionally, whether there is a relation between reduced bispec-trum and baryon fraction, analogous to the one found for the powerspectrum.In Fig. 11 we show how, indeed, in the case of the BAHAMAS MNRAS000 , 1–14 (2020) atter bispectrum with baryons high-AGN simulation, the fit of the bispectrum marginally improvesthe baryon fraction estimation. The fact that both the best-fittingparameter set can accurately reproduce the power spectra of thehydrodynamical simulations, but predict slightly different baryonfractions, can be a hint of some degeneracies between parameterswhich is broken when including the bispectrum information. Forthe Illustris, the opposite is true: the gas fraction in clusters differsmore from its true value. This likely points to the fact that someof the baryonification assumptions somewhat break for extremefeedback scenarios. This could be related to gas fractions that arenot monotonic with halo mass, or that these events affect gas beyondthe boundaries of a halo (a process not included in our model). Onthe other hand, we note that the gas fractions in these simulationsare in clear tension with observations which prefer values ≈ . .Regardless of the simulation, the bottom panel of Fig. 11 showsthat the baryonic effects on the reduced equilateral bispectrum cor-relate with the baryonic fraction: the smaller the baryon fraction,the larger the bispectrum enhancement. Remarkably, the predictionfrom our model, when fitted with a simple linear regression, showsa trend as tight as the one found in the power spectrum ( ). Tohave an idea of the predictions from hydrodynamical simulations,we infer the baryon fractions from the power spectra measurementsusing the fitting function provided by van Daalen et al. (2019),and combine them to the measurement of the reduced bispectraenhancement. By doing so, we find that all the predictions are stillwithin , except for BAHAMAS high-AGN, which is slightly offbut still well within .It would be very interesting to extend the analysis of van Daalenet al. (2019) to the bispectrum, to check if including a vast num-ber of hydrodynamical simulations the relation still holds with alow scatter. Nevertheless, we stress that, a priori, the BCM doesnot predict a tight relation between baryon fraction and clustering.In A20 (Fig. 8) it was shown that the baryon fraction-clusteringrelation is more relaxed when considering the full BCM parame-ters ranges. Interestingly, it seems that the calibration and subgridphysics with which hydrodynamical simulations are run, translatesinto constraints and degeneracies of the BCM parameters, and thusconstrain the baryon fraction. In this paper, we have used a combination of cosmology scaling andbaryonification algorithms, to reproduce with a negligible compu-tational time the density fields of various hydrodynamical simula-tions, up to very small and non-linear scales ( k = 5 h Mpc − ) andfor two- and three-point statistics.Below we summarise our main findings: • Baryonic physics causes an enhancement in the reduced equi-lateral bispectrum at all the scales considered, roughly monotoni-cally with the strength of the feedback mechanisms; • It is possible to simultaneously reproduce the baryonic effectson the power spectrum and on reduced bispectrum (with and − precision, respectively), as measured in EAGLE, Illustris,Illustris TNG-300, and three different AGN implementations ofBAHAMAS, • In contrast, a baryon model tuned to only reproduce the powerspectrum, can lead to up to ∼ discrepancies in the reducedbispectrum; • We find that a double power-law gas density profile is flexibleenough to reproduce the bump at small scales measured in thebispectrum of some hydrodynamical simulations (see Foreman et al. 2019). It appears thus that an additional modelling of gas overdensityat relatively small scales is superfluous. • The model parameters that best fit the power spectrum andequilateral bispectrum also predict changes to the squeezed config-urations at the ∼ level. • The baryon parameters are not redshift independent; ignoringtheir time dependence results in a inaccuracy in the powerspectrum, and − in the reduced bispectrum, up to z = 2 ; • Analysing the best-fitting models to the hydrodynamical sim-ulations, we find a correlation between baryonic effects on the bis-pectrum and baryon fraction inside haloes, similar to the one for thepower spectrum found in van Daalen et al. (2019),Overall, our results support the physical soundness (as wellas our specific numerical implementation) of baryonification algo-rithms. This also encourages its use not only in spherically-averaged2-point statistics, but also in cross-correlations and in other statisticssuch as peak counts.The next generation surveys will produce a huge amount ofdata, which is only partially interpretable with the current theoreticalmodels. This paper is a contribution to the effort to overpass modelsbased on only gravitational interactions, and fully exploit the dataup to higher-order statistics. We anticipate that our approach will bea valid tool for a fast production of mock density fields, accurate tovery small scales and statistics of order higher than 2-point, usefulfor pipeline validation, blind comparisons or for direct exploitationof the data, e.g. marginalising over baryonic effects.
ACKNOWLEDGEMENTS
The authors acknowledge the support of the E.R.C. grant 716151(BACCO). C.H.-M. acknowledges support from the Spanish Min-istry of Economy and Competitiveness (MINECO) through theprojects AYA2015-66211-C2-2 and PGC2018-097585-B-C21. SCacknowledges the support of the “Juan de la Cierva Formación” fel-lowship (FJCI-2017-33816). We thank Simon Foreman for makingpublic the power spectra and bispectra measurements used in thiswork, as well as the code bskit . We acknowledge the Illustris, Illus-tris TNG, BAHAMAS, and EAGLE teams, for providing/makingpublic the data of their hydrodynamical simulations. We are grate-ful to Alex Barreira and Simon Foreman for carefully reading thedraft, and helping us to improve the manuscript with their preciousfeedback. We thank Yetli Rosas-Guevara for providing us with themeasurements of the mass density profiles in the Illustris TNG-300simulation. G.A. thank Lurdes Ondaro Mallea and Marcos PellejeroIbañez for useful discussions.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Aihara H., et al., 2018, PASJ, 70, S4Angulo R. E., Hilbert S., 2015, MNRAS, 448, 364Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., FrenkC. S., 2012, MNRAS, 426, 2046MNRAS , 1–14 (2020) G. Aricò et al.
Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M.,Stücker J., 2020, arXiv e-prints, p. arXiv:2004.06245Aricò G., Angulo R. E., Hernández-Monteagudo C., Contreras S., ZennaroM., Pellejero-Ibañez M., Rosas-Guevara Y., 2020, MNRAS, 495, 4800Barreira A., Nelson D., Pillepich A., Springel V., Schmidt F., Pakmor R.,Hernquist L., Vogelsberger M., 2019, MNRAS, 488, 2079Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57Benitez N., et al., 2014, preprint, ( arXiv:1403.5237 )Chisari N. E., et al., 2019, arXiv e-prints,Colombi S., Jaffe A., Novikov D., Pichon C., 2009, MNRAS, 393, 511Contreras S., Zennaro R. E. A. M., Aricó G., Pellejero-Ibañez M., 2020,arXiv e-prints, p. arXiv:2001.03176Crain R. A., et al., 2015, MNRAS, 450, 1937DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036Foreman S., Coulton W., Villaescusa-Navarro F., Barreira A., 2019, arXive-prints, p. arXiv:1910.03597Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018,AJ, 156, 160Hellwing W. A., Schaller M., Frenk C. S., Theuns T., Schaye J., Bower R. G.,Crain R. A., 2016, MNRAS, 461, L11Jenkins A., et al., 1998, ApJ, 499, 20Kennedy J., Eberhart R., 1995, in Proceedings of ICNN’95 - InternationalConference on Neural Networks. pp 1942–1948 vol.4Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, AstronomyLetters, 44, 8Laureijs R., et al., 2011, preprint, ( arXiv:1110.3193 )Marinacci F., et al., 2018, MNRAS, 480, 5113McAlpine S., et al., 2016, Astronomy and Computing, 15, 72McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465,2936McCarthy I. G., Bird S., Schaye J., Harnois-Deraps J., Font A. S., vanWaerbeke L., 2018, MNRAS, 476, 2999Mead A. J., Peacock J. A., 2014a, MNRAS, 440, 1233Mead A. J., Peacock J. A., 2014b, MNRAS, 445, 3453Mead A. J., Peacock J. A., Lombriser L., Li B., 2015, MNRAS, 452, 4203Naiman J. P., et al., 2018, MNRAS, 477, 1206Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Nelson D., et al., 2018, MNRAS, 475, 624Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2Pillepich A., et al., 2018, MNRAS, 475, 648Planck Collaboration et al., 2018, arXiv e-prints,Renneby M., Hilbert S., Angulo R. E., 2018, MNRAS, 479, 1100Ruiz A. N., Padilla N. D., Domínguez M. J., Cora S. A., 2011, MNRAS,418, 2422Schaye J., et al., 2015, MNRAS, 446, 521Schneider A., Teyssier R., 2015, J. Cosmology Astropart. Phys., 12, 049Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., AmaraA., Refregier A., 2019, J. Cosmology Astropart. Phys., 2019, 020Schneider A., Stoira N., Refregier A., Weiss A. J., Knabenhans M., StadelJ., Teyssier R., 2020, J. Cosmology Astropart. Phys., 2020, 019Scoccimarro R., 2000, ApJ, 544, 597Sefusatti E., Komatsu E., 2007, Phys. Rev. D, 76, 083004Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, MN-RAS, 460, 3624Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F.,Nelson D., Hernquist L., 2015, MNRAS, 452, 575Springel V., 2005, MNRAS, 364, 1105Springel V., et al., 2018, MNRAS, 475, 676The EAGLE team 2017, arXiv e-prints, p. arXiv:1706.09899Troxel M. A., et al., 2018, Phys. Rev. D, 98, 043528Verde L., Treu T., Riess A. G., 2019, Nature Astronomy, 3, 891Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L.,2013, MNRAS, 436, 3031Vogelsberger M., et al., 2014a, MNRAS, 444, 1518Vogelsberger M., et al., 2014b, Nature, 509, 177Watson D. F., Conroy C., 2013, ApJ, 772, 139Wong K. C., et al., 2020, MNRAS,
Figure A1.
Upper panel:
Baryon suppression of the matter powerspectrum (left) and reduced bispectrum (right), defined as S ( k ) ≡ T ( k ) BCM /T ( k ) GrO for T ( k ) = P ( k ) , Q ( k ) , at z = 0 . Solid, dashed,dashed-dotted and dotted lines are computed with simulations of box side , , and h − Mpc and 1536 , 768 , 384 and 192 parti-cles, respectively. Colors are referred to different halo mass bins, expressedin decimal logarithm of h − M (cid:12) , with which the baryon corrections havebeen computed, according to the legend. Lower panel:
Difference in sup-pression between a paired and fixed simulation and a single realisation, forthe four different volumes specified in the legend of the upper panel.
Figure A2.
Baryon suppression of the matter power spectrum (left) andreduced bispectrum (right), defined as S ( k ) ≡ T ( k ) BCM /T ( k ) GrO for T ( k ) = { P ( k ) , Q ( k ) } , at z = 0 . Solid, dashed, dashed-dotted and dottedlines are computed with simulations of box side h − Mpc and 768 ,512 , 384 and 256 particles, respectively. Colours are referred to differenthalo mass bins, expressed in decimal logarithm of h − M (cid:12) , with which thebaryon corrections have been computed, according to the legend.Zennaro M., Angulo R. E., Aricò G., Contreras S., Pellejero-Ibáñez M.,2019, MNRAS, 489, 5938van Daalen M. P., McCarthy I. G., Schaye J., 2019, arXiv e-prints, APPENDIX A: CONVERGENCE TEST
In this Appendix, we show the tests we have performed to assure thatthe baryonic effects on the clustering measurements have converged.First, we test the convergence with the simulation box size.For this, we have used our suite of simulations with h − Mpc , h − Mpc , h − Mpc , and h − Mpc of box side, with N = 192 , N = 384 , N = 768 and N = 1536 cubic particles, re- MNRAS000
In this Appendix, we show the tests we have performed to assure thatthe baryonic effects on the clustering measurements have converged.First, we test the convergence with the simulation box size.For this, we have used our suite of simulations with h − Mpc , h − Mpc , h − Mpc , and h − Mpc of box side, with N = 192 , N = 384 , N = 768 and N = 1536 cubic particles, re- MNRAS000 , 1–14 (2020) atter bispectrum with baryons spectively. All the simulations have same force and mass resolution,and share the same initial conditions. For each different volume wehave run two simulations with fixed amplitude and shifted phasesas reported in §2.In Fig. A1 we show the suppression S ( k ) , defined as the ra-tio between baryonified and gravity-only matter power spectra andequilateral bispectra, measured in the four different boxes at z = 0 .We have also separated the contribution to the clustering ofdifferent halo masses, to get more insight on the origin of the dis-crepancies between the different boxes. As expected, we note that theboxsize does not affect sensibly haloes of M ≤ h − M (cid:12) . How-ever, the abundance of massive haloes ( M = 10 − h − M (cid:12) )varies consistently among the various boxes, leading to discrepan-cies in the baryonic effects that are still within in the powerspectrum, but slightly higher in the bispectrum ( − ).In the bottom panels of Fig. A1 we display the impact of using a paired and fixed simulation against a single realisation. Also in thiscase, the biggest impact is found in the bispectrum, with a maximumof ≈ . bias when using a h − Mpc box, whereas we detecta maximum of ≈ in the power spectrum. In the analysis, wemake use of a single realisation of the h − Mpc box, which isshown to be converged within .We have performed also a mass resolution test, by using foursimulation with the same box size, h − Mpc , and different num-ber of particles: N = 256 , N = 384 , N = 512 , N = 768 particles. Also in this case, we split the contribution of differenthalo mass bins. As shown in Fig. A2, we have found that resolu-tion effects are larger in large haloes, and their impact in the powerspectrum and bispectrum is within ≈ . APPENDIX B: FOLDING OF THE PARTICLEDISTRIBUTION
Measuring the three-point clustering with the classical Fourier es-timators can be very expensive in terms of memory and CPU, es-pecially when using covering larger dynamical ranges. In fact, it iseasy to see that, being k Ny = πN g /L box the Nyquist frequency ofthe grid, increasingly large number of grid points N g are requiredto get a given accuracy at a fixed wavenumber, when using pro-gressively larger simulation boxes L box . Additionally, when using“interlacing” to suppress aliasing, the number of grids used must bedoubled (Sefusatti et al. 2016).However, since our measurements are limited by discretenessnoise (and not cosmic variance), we can obtain accurate estimates ofFourier statistics on small scales by folding the density field (Jenkinset al. 1998; Colombi et al. 2009). The idea is to fold the particle dis-tribution by re-applying the periodic boundary conditions assuminga new boxsize L (cid:48) = L/f , where we call f the number of foldings.If L (cid:48) is large enough to assure that the modes inside the new boxare uncorrelated, we can measure in principle the clustering froma new effective fundamental wavenumber k (cid:48) f = 2 π/L (cid:48) up to a neweffective Nyquist frequency, given by k (cid:48) Ny = πN g /L (cid:48) . For instance,by folding the box 4 times, we will get to Nyquist frequency 4 timeshigher.In Fig. A3 we apply this technique to our h − Mpc simulation, folding the particles up to 6 times, and reaching a k Ny ≈ h Mpc − with a and a grid for the power spec-trum and the bispectrum, respectively. Even using a TSC schemeon interlaced grids, we note that it is safer to use the measurementsup to k = k (cid:48) Ny / in the bispectrum. Also, the measurement of thelargest modes of the folded box are noisy because they are sparsely Figure A3.
Left panels:
Matter power spectrum measured in a (inter-laced) mesh, folding the box up to 6 times, following the technique explainedin the text (coloured dots). For comparison, the matter power spectrum mea-sured with a mesh and not folding the box is plotted as a black solidlines. The equivalent Nyquist frequencies for each folded box is plotted asa solid line. In the bottom panel, we display the ratio between the powerspectrum measure with a and a mesh. Right panels:
Similarly tothe left panels, we display the measured bispectra using a mesh andthe folding of the box, and compare with a not-folded mesh. Figure A4.
Left panels:
Ratio of baryonified and GrO matter power spectra S ( k ) , measured in a (interlaced) mesh, folding the particles up to 6times, following the technique explained in the text (coloured dots). Forcomparison, the measurements with a mesh is plotted as a black solidlines. The equivalent Nyquist frequencies for each folded box is plotted asa coloured solid line. In the bottom panel, we display the difference of theratios ∆ S ( k ) measured with a and a mesh. Right panels:
Similarto the left panels, but displaying the matter bispectrum instead of the powerspectrum.MNRAS , 1–14 (2020) G. Aricò et al. sampled; for this reason, it is convenient to discard these modes,taking for instance wavenumbers k > k (cid:48) f .Using these precautions, we show in Fig. A4 that using thistechnique we can achieve an accuracy well within in the estima-tion of the ratios, using a small fraction − of the computationalresources. Although it is common to use the folding technique tocompute power spectra, to our knowledge, this is the first time it hasbeen shown to be accurate for bispectrum measurements. This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000