Semianalytic methods can generate baryoncorrected fields from Nbody simulations baryonification and are rapidly becoming a ubiquitous tool in modeling structure formation on nonlinear scales. We extend this formalism to consistently model the weak lensing and thermal SunyaevZeldovich tSZ fields directly on the fullsky with an emphasis on higherorder correlations. We use the auto and cross Nthorder moments with N in 2 3 4 as a summary statistic of the lensing and tSZ fields and show that our model can jointly fit these statistics measured in IllustrisTNG to within measurement uncertainties for scales above gtrsim 1 rm Mpc and across multiple redshifts. The model predictions change only minimally when including additional information from secondary halo properties such as halo concentration and ellipticity. Each individual moment is dependent on halos of different mass ranges and has different sensitivities to the model parameters. A simulationbased forecast on the ULAGAM simulation suite shows that the combination of all moments measured from current and upcoming lensing and tSZ surveys can jointly constrain cosmology and baryons to high precision. The lensing and tSZ field are sensitive to different combinations of the baryonification parameters with degeneracy directions that are often orthogonal and the combination of the two fields leads to significantly better constraints on both cosmology and astrophysics. Our pipeline for maplevel baryonification is publicly available at httpsgithub.comDhayaaAnbajaganeBaryonForge.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Maplevel baryonification Efficient modelling of higherorder correlations in the weak lensing and thermal SunyaevZeldovich fields
{'No Header': 'Preprint typeset using L A T E X style openjournal v. 09/06/15', 'Map-level baryonification: Efficient modelling of higher-order correlations in the weak lensing and thermal Sunyaev-Zeldovich fields': 'Dhayaa Anbajagane ( ) 1 , 2 , Shivam Pandey 3 , 4 , and Chihway Chang 1 , 2  \n- 1 Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA\n- 2 Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA\n- 3 Department of Physics, Columbia University, 538 West 120th Street, New York, NY, USA 10027, USA and\n- 4 Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA  \nVersion December 11, 2024', 'ABSTRACT': "Semi-analytic methods can generate baryon-corrected fields from N-body simulations ('baryonification') and are rapidly becoming a ubiquitous tool in modeling structure formation on non-linear scales. We extend this formalism to consistently model the weak lensing and thermal Sunyaev-Zeldovich (tSZ) fields directly on the full-sky, with an emphasis on higher-order correlations. We use the auto- and cross𝑁 th-order moments, with 𝑁 ∈ { 2 , 3 , 4 } , as a summary statistic of the lensing and tSZ fields, and show that our model can jointly fit these statistics measured in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG to within measurement uncertainties, for scales above ≳ 1 Mpc and across multiple redshifts. The model predictions change only minimally when including additional information from secondary halo properties, such as halo concentration and ellipticity. Each individual moment is dependent on halos of different mass ranges and has different sensitivities to the model parameters. A simulationbased forecast on the U/l.pc/a.pc/g.pc/a.pc/m.pc suite shows that the combination of all moments, measured from current and upcoming lensing and tSZ surveys, can jointly constrain cosmology and baryons to high precision. The lensing and tSZ field are sensitive to different combinations of the baryonification parameters, with degeneracy directions that are often orthogonal, and the combination of the two fields leads to significantly better constraints on both cosmology and astrophysics. Our pipeline for map-level baryonification is publicly available at https://github.com/DhayaaAnbajagane/BaryonForge .", '1 INTRODUCTION': "Measurements of weak lensing, which is a direct probe of the cosmic density field, have provided some of the best constraints on the properties of the Universe, both on the initial conditions and on their subsequent evolution (Abbott et al. 2022; Asgari et al. 2021; Secco et al. 2022a; Amon et al. 2022; More et al. 2023; Omori et al. 2023; Chang et al. 2023; Madhavacheril et al. 2024). These constraints are expected to significantly improve with the advent of surveys like the Dark Energy Spectroscopic Instrument DESI Collaboration et al. (2016), Vera C. Rubin Observatory (The LSST Dark Energy Science Collaboration et al. 2018), the Euclid mission (Racca et al. 2016), the Roman space telescope (Spergel et al. 2015) and more, all of which will provide vastly richer datasets and therefore enable measurements of the small-scale density field  \nwith improved statistical precision.  \nHowever, an ongoing limitation in using such small-scale measurements is modeling the impact of baryons ('baryonic imprints') on the density field. Baryons, which constitute ≈ 15% of the total matter in our Universe, with the rest being accounted for by dark matter (DM), experience a vast set of astrophysical process across a range of scales. Such processes alter the distribution and thermodynamics of these baryons, and also the DM phase space via the gravitational coupling between the two components ( e.g. , Gnedin et al. 2004; Abadi et al. 2010; Duffy et al. 2010; Anbajagane et al. 2022a; Shao & Anbajagane 2024). The signatures of these processes can be prevalent even on quasi-linear scales (roughly a few times a halo radius) and will be a significant component of any smallscale signal (i.e. the density near/within a halo). Their presence can bias any cosmological analyses that do not account for this effect ( e.g. , Secco et al. 2022a; Amon & Efstathiou  \n2022; Gatti et al. 2022; Anbajagane et al. 2023). This limitation in scales is further accentuated by the fact that many novel cosmological signatures are imprinted onto these small scales. Thus, efficient models of these baryonic imprints have the potential to not only improve the constraining power for current models of interest, such as Λ CDM and 𝑤 CDM, but also open the opportunity to constrain extended models like those of primordial non-gaussianities (Coulton et al. 2022; Anbajagane et al. 2024b; Jung et al. 2023; Goldstein et al. 2024).  \nThe most precise model of baryonic imprints comes from studying the formation of structure in various hydrodynamical simulations which model the evolution of gas, stars, black holes etc. in addition to the DM ( e.g. , McCarthy et al. 2017; Springel et al. 2018). The simulations have a finite resolution scale that is larger than many relevant astrophysical processes (such as supernovae explosions, or supermassive black hole accretion, to name a few) and must therefore implement these processes as 'sub-grid' physics that approximate their effects (see Vogelsberger et al. 2020, for a review). The subgrid models of each simulation vary in the equations that are solved, and in how the equations are parameterized. Studies of these simulations, however, have showcased the variety of predictions manifesting from these different, but often equally realistic, choices in the models of galaxy formation ( e.g. , Anbajagane et al. 2020; Lim et al. 2021; Lee et al. 2022; Cui et al. 2022; Stiskalek et al. 2022; Anbajagane et al. 2022a,b; Shao et al. 2023; Gebhardt et al. 2023). The type and amplitude of these differences vary significantly depending on the exact properties - such as the gas mass/temperature/pressure, the stellar mass/age/metallicity etc. - being studied and the regime in mass, redshift, environment etc. being explored. Thus, while the simulations are valuable in capturing the range of possible baryonic imprints, one cannot use them to robustly model these imprints in the density field. Furthermore, studies on the thermodynamic properties of gas also find differences between the measurements from data and the predictions from these hydrodynamic simulations ( e.g. , Hill et al. 2018; Amodeo et al. 2021; Pandey et al. 2022; Anbajagane et al. 2022c).  \nGiven the current challenges in modelling baryons in simulations, it is advantageous to instead use a phenomenological model. Such a model does not require specifying any sub-grid physics and instead requires inputs that are observables, or are quantities directly obtained from observables, of the different baryonic matter components. Baryonification is a method through which N-body simulations - which evolve collisionless matter (dark matter) under gravity and accurately capture all non-linear evolution - are modified through phenomonological models such that the new density field mimics one that has baryonic imprints (Schneider & Teyssier 2015; Schneider et al. 2019). This method has been used extensively to model the matter power spectrum with baryonic imprints, and has  \nbeen shown to have the flexibility needed to capture the behaviors ranging across many different simulations (Schneider et al. 2019; Giri & Schneider 2021; Aric'o et al. 2021b). It has already been used to analyse 2-point correlation functions measured in widefield surveys (Chen et al. 2023; Aric'o et al. 2023; Bigwood et al. 2024).  \nThe fiducial baryonification method is performed by displacing the positions of particles in the N-body simulations to alter the overall density distribution in the volume. This procedure has been computationally tractable for building emulators that translate the model parameters to prediction a given statistic in real space, e.g. , the 3D power spectrum such as done in Giri & Schneider (2021); Aric'o et al. (2021b). These emulators can then be used in standard analysis pipelines to predict, for example, the weak lensing correlation functions ( e.g. , Aric'o et al. 2023). However, surveys are rapidly adopting simulation-based analyses as a way to accurately forward model systematic effects, and also to incorporate more information from the field ( e.g. , Gatti et al. 2023, 2024b; Jeffrey et al. 2024; Harnois-Deraps et al. 2024; Cheng et al. 2024). Some work has been done to show the baryonification model generates fields that accurately model baryonic imprints on higher-order statistics, such as the density field bispectrum (Aric'o et al. 2021a) and the peaks of the weak lensing field (Lee et al. 2023). These higher-order statistics extend upon the simpler, two-point statistics - which describe the correlation between two points in space, e.g. , the power spectrum - by also extracting correlations between three or more points ( e.g. , Fluri et al. 2019; Gatti et al. 2020; Zurcher et al. 2021; Fluri et al. 2022; Euclid Collaboration et al. 2023; Anbajagane et al. 2023; Gatti et al. 2024a; Jeffrey et al. 2024).  \nThe use of baryonification in simulation-based analyses of surveys is limited by the computational expense of the method. Survey datasets frequently require lightcone maps covering a large sky fraction for their predictions (upcoming surveys will need ≈ 15 , 000 deg 2 coverage) and these maps are constructed from O( 10 2 ) 3D snapshots/boxes in a simulation ( e.g. , Kacprzak et al. 2023; Anbajagane et al. 2024b; Jeffrey et al. 2024). Applying baryonification on each box is an computationally expensive procedure - and also a memory intensive one, since it requires storing the particles of each snapshot - given analyses often require O( 10 3 -10 4 ) such simulations. An attractive alternative then is to perform baryonification directly on pixels of 2D maps. This has been explored by Fluri et al. (2019) for performing simulation-based inference of weak lensing data from the Kilo-Degree Survey (Kuijken et al. 2015). An aspect of our work focuses on updating this method, by robustly accounting for the pixel window functions and line-of-sight projections and then validating its performance on simulations.  \nWe have so far only discussed the total density field,  \nwhereas other probes - often more closely linked to the baryonic matter components - are actively used to constrain astrophysics and cosmology. Of these, one of the most commonly used and actively studied ones is the thermal SunyaevZeldovich (tSZ) field (Sunyaev & Zeldovich 1972), observed using millimeter surveys such as the South Pole Telescope (Carlstrom et al. 2011; Benson et al. 2014), the Atacama Cosmology Telescope (Fowler et al. 2007; Thornton et al. 2016), and the Planck satellite (Planck Collaboration et al. 2020). The tSZ effect is the inverse Compton scattering of CMB photons with energetic electrons along the line of sight. It probes the integrated gas pressure and therefore the thermodynamics of baryons in and around the halo (see Carlstrom et al. 2002; Mroczkowski et al. 2019, for reviews). As a result, it is advantageous to include the tSZ probe in analyses, both for its cosmological information - as it, like the total matter field, is sensitive to the initial conditions and their subsequent evolution - and even more so for its information on the gas thermodynamics. However, pursuing this direction requires a consistent, common model that predicts the impact of astrophysical processes on both the distribution of matter and the thermodynamics of baryons. Studies of baryonification have thus far been exclusively on the density field (Schneider et al. 2019; Giri & Schneider 2021; Aric'o et al. 2021b; Lee et al. 2023; Aric'o et al. 2023), with the inclusion of tSZ only recently being explored (Aric'o & Angulo 2024). /one.sup  \nPandey et al. (2024, henceforth, P24) show that the halo profiles used in the existing baryonification model can be extended, with simple modifications, to also predict the gas pressure profiles (and gas number density profile) in addition to the total matter density profiles. To et al. (2024, see their Figure 1) have further demonstrated that this model can consistently predict the tSZ-halo mass scaling relation alongside the density two-point correlations. Thus, the original baryonification model of Schneider et al. (2019, henceforth, S19) can be built on to now also predict the tSZ observables (as well as other observables that trace the gas number density, such as the kinematic SZ effect or the X-ray flux). While P24 focus on using the halo model approach to predict two-point measurements of weak lensing and the tSZ, the model remains to be explored and validated at the map-level for different summary statistics that combine information from the lensing and tSZ fields.  \nThis work explores exactly this aspect, by both building and validating a map-level baryonification model that can be used to predict higher-order correlations. We do this in three ways: (i) we first validate the map-level baryonification method in modelling the density and tSZ fields, using the moments of the fields up to fourth order as our summary statistic, (ii) then we test the sensitivity of our predictions to the in-  \nclusion of secondary halo properties such as concentration and ellipticity, as well as additional features in the input halo profiles such as cosmological shocks and non-thermal pressure, and (iii) finally, we showcase the ability of existing and upcoming surveys to constrain these baryonic effects, either separately or in a joint analysis with cosmological parameters varied as well.  \nOf particular note is that our modelling pipeline is built atop the open-source C/o.pc/r.pc/e.pc C/o.pc/s.pc/m.pc/o.pc/l.pc/o.pc/g.pc/y.pc L/i.pc/b.pc/r.pc/a.pc/r.pc/y.pc suite (CCL, Chisari et al. 2019) developed for the Rubin LSST Dark Energy Science Collaboration (DESC). This choice both allows our methods to be easily accessible/usable by the community, by borrowing the many user accesibility features built for the suite, and more importantly allows the halo profiles used in this method to also be propagated through other halo model tools already built with CCL.  \nWe organize this paper as follows: Section 2 describes the baryonification model used in this work. Section 3 describes the simulations, the summary statistics used to model and validate our method, and the quantification of measurement uncertainty and model best-fit. Section 4 describes a subset of the tests we perform to validate our modelling pipeline and characterize its behavior. Finally, section 5 details the power of current and upcoming surveys in constraining the parameters associated with these baryonic effects. We conclude in Section 6.  \nA number of appendices provide further details on the methods we used: the lensing and tSZ forward modelling approach in Appendix A, the sensitivity of the predictions to additional model and methodology choices in Appendix B and C, the redshift evolution of the baryonification predictions in Appendix D, the dependence between Fisher information and baryonification model-complexity in Appendix E, and finally, some salient computation details on the method, include the runtimes, in Appendix F.", '2 BARYONIFICATION': 'We first describe the halo model used in this work in Section 2.1. Then, we detail the procedures used to include baryon signatures into the different large-scale structure fields we consider: the matter density field (Section 2.2) and the tSZ field (Section 2.3). Our model follows from that presented in P24, which builds on the density field-only baryonification model of S19 and Giri & Schneider (2021). Following the CCL convention, all distance scales used in defining the profiles are written in comoving Mpc units, and all masses as in M ⊙ . We will frequently use 𝑀 200c, which is a spherical overdensity mass defined as the mass contained within a halo-centric sphere of average density ⟨ 𝜌 ⟩ = 200 𝜌 𝑐 ( 𝑧 ) , with 𝜌 𝑐 ( 𝑧 ) being the critical density of the universe at redshift 𝑧 . The radius of  \nthe sphere is 𝑅 200c.  \nA number of salient computational details about the method are described in Appendix F, including the characteristic runtimes of the pipeline. The baryonified lensing and tSZ maps for a simulation /two.sup can be generated on the order of minutes on an Intel broadwell chip with 40 cores; see Appendix F for more details.', '2.1 Dark matter baryon (DMB) halo model': 'The DMBhalomodelisasumovermultiple different components. Each component is parameterized in a different way, with the parameterization informed by prior work in simulations and/or observations. We now describe each component below:  \nThe dark matter profile is modelled by a simple NFW form (Navarro et al. 1997) with an additional truncation term that scales as 𝑟 -4 ,  \n𝜌 NFW ( 𝑟 ) = 𝜌 0 GLYPH<18> 𝑟 𝑟 𝑠 GLYPH<19> -1 GLYPH<18> 1 + 𝑟 𝑟 𝑠 GLYPH<19> -2 GLYPH<18> 1 + 𝑟 2 𝑟 2 𝑡 GLYPH<19> -2 , (1)  \nwhere 𝜌 0 is a normalization coefficient set by requiring that the mass enclosed within 𝑅 200c is 𝑀 200c,  \n𝜌 0 = 𝑀 200c GLYPH<20> ∫ 𝑅 200c 0 4 𝜋𝑟 2 𝜌 NFW ( 𝑟 ) 𝑑𝑟 GLYPH<21> -1 . (2)  \nNote that this profile is implicitly a function of 𝑟 𝑠 , or alternatively the halo concentration 𝑐 200c = 𝑅 200c / 𝑟 𝑠 . The truncation term is necessary as the baryonification model requires as estimate of the total enclosed mass at 𝑟 → ∞ and this is a divergent quantity for the standard NFW profile as the profile scales as 𝑟 -3 at 𝑟 ≫ 𝑟 𝑠 . We follow S19 and set 𝑟 𝑡 = 4 𝑅 200c (see Table 1).  \nThetwohalo profile , which describes the extended matter distribution around a halo contributed to by its neighboring structures, is defined similar to S19, as  \n𝜌 2h = 𝜌 𝑚 ( 𝑧 ) GLYPH<20> 1 + 𝜉 𝑚𝑚 𝑏 ( 𝜈 200c ) GLYPH<21> , (3)  \nwith 𝜌 𝑚 ( 𝑧 ) being the comoving matter density at redshift 𝑧 , 𝜉 𝑚𝑚 the linear matter correlation function, and 𝑏 ( 𝜈 200c ) is the halo bias given by  \n𝑏 ( 𝜈 200c ) = 1 + 𝑞𝜈 2 200c -1 𝛿 𝑐 + 2 𝑝 𝛿 𝑐 ( 1 + 𝑞𝜈 2 200c ) 𝑝 . (4)  \nHere, 𝛿 𝑐 = 1 . 686 / 𝐷 ( 𝑧 ) , with 𝐷 being the growth factor normalized to 𝐷 ( 𝑧 = 0 ) = 1, and 𝜈 200c is the peak height of  \nthe halo, and is computed as,  \n𝜈 200c = 𝛿 𝑐 / 𝜎 ( 𝑅 200c ) , (5)  \nwith 𝜎 ( 𝑅 200c ) being the root-mean square of the density field smoothed with a tophat of radius 𝑅 200c. We fix 𝑞 = 0 . 707 and 𝑝 = 0 . 3, following S19 which in turn is based on the findings of Sheth & Tormen (1999).  \nThe NFW profile and the two halo profile can be summed to predict the total matter distribution in the dark matter-only (DMO) case. Now, we describe the profiles of the baryonic components.  \nThe stellar profile model follows Mohammed et al. (2014); Schneider & Teyssier (2015) and is a simple powerlaw with an exponential cutoff,  \n𝜌 star = 𝑀 star 4 𝜋 3 / 2 𝑅 ℎ 1 𝑟 2 exp GLYPH<20> -GLYPH<18> 𝑟 2 𝑅 ℎ GLYPH<19> 2 GLYPH<21> , (6)  \nwhere 𝑅 ℎ = 𝜖 ℎ 𝑅 200c is the half-light radius of the stellar halo, and 𝑀 star is the normalization of the profile obtained as,  \n𝑀 star = 𝑓 cga 𝑀 tot , (7)  \nwhere we have defined a theoretical total halo mass, 𝑀 tot ( 𝑟 → ∞) , through the expression  \n𝑀 𝑋 ( 𝑟 ) = ∫ 𝑟 0 d 𝑟 4 𝜋𝑟 2 𝜌 𝑋 ( 𝑟 ) . (8)  \nThe central galaxy fraction, 𝑓 cga, is then given by  \n𝑓 cga = 2 𝐴 GLYPH<20> GLYPH<18> 𝑀 200c 𝑀 1 GLYPH<19> 𝜏 cga + GLYPH<18> 𝑀 200c 𝑀 1 GLYPH<19> 𝜂 cga GLYPH<21> -1 . (9)  \nHere, the two power-law exponents are given by 𝜏 cga = 𝜏 + 𝜏 𝛿 and 𝜂 cga = 𝜂 + 𝜂 𝛿 , where 𝜏 and 𝜂 correspond to the scaling of the total stellar fraction, and 𝜏 𝛿 and 𝜂 𝛿 are a correction to convert the total stellar fraction scaling to the central galaxy stellar fraction scaling. The factor of 2 ensures that for 𝑀 200c = 𝑀 1, we have 𝑓 cga = 𝐴 . The original baryonification model of S19 used a different version of Equation (9), which included 𝜂 and 𝜂 𝛿 but did not have 𝜏 and 𝜏 𝛿 . These latter parameters control the slope of the stellar-mass halo-mass relation for low-mass halos and is included in this work. /three.sup In practice, the 𝜏 and 𝜏 𝛿 parameters have no impact in our work as the signal for both weak lensing and tSZ is generated by halos that are above the 𝑀 1 = 3 × 10 11 M ⊙ mass scale (see Table 1 for the fiducial values of different parameters), and such halos are sensitive only to 𝜂, 𝜂 𝛿 and not 𝜏, 𝜏 𝛿 . However, we include this parameter here as it is a simple extension that generalizes the model for smaller mass halos.  \n/three.sup S19 note that analyses of weak lensing only use halos above 𝑀 200c > 10 12 M ⊙ , for which their single power-law model is accurate and adequate. We extend this to the double power-law model traditionally used in observational and theoretical analyses of galaxies ( e.g. , Behroozi et al. 2019, see their Figure 9). This extension allows the baryonification model to be used consistently with halos of lower masses.  \nAs discussed before, the quantity 𝑀 tot does not diverge now due to our addition of a truncation radius, 𝑟 𝑡 , in Equation (1). Without this factor, 𝑀 tot = 𝑀 NFW ( 𝑟 → ∞) → ∞ . This is relevant as we use this total mass, 𝑀 tot, alongside the relevant mass fractions discussed above, to define the total mass of stars/gas mass in a halo. The baryonification model assumes that 𝑀 tot for a given halo is the same in Universes with and without the presence of baryons. This assumption is valid as the redistribution of matter is mostly localized around the halo, and for sufficiently large radii (we use 𝑟 → ∞ in definition 𝑀 tot) the enclosed mass is constant even after this redistribution process ( e.g. , Ayromlou et al. 2022; Gebhardt et al. 2023).  \nThe gas profile is a modified version of the generalized NFW profile (Nagai et al. 2007, see their Equation A1), with two distinct length scales, /four.sup  \n𝜌 gas = 𝜌 gas , 0 GLYPH<18> 1 + 𝑟 𝑅 co GLYPH<19> -𝛽 GLYPH<18> 1 + GLYPH<18> 𝑟 𝑅 ej GLYPH<19> 𝛾 GLYPH<19> -𝛿 -𝛽 𝛾 . (10)  \nHere, 𝑅 co = 𝜃 co 𝑅 200c is the length scale of the gas core, and 𝑅 ej = 𝜃 ej 𝑅 200c is that of the ejected gas. The scaling 𝛽 , which controls the slope between 𝑅 co < 𝑟 < 𝑅 ej, has an additional mass dependence as,  \n𝛽 = 3 ( 𝑀 200c / 𝑀 𝑐 ) 𝜇 𝛽 1 + ( 𝑀 200c / 𝑀 𝑐 ) 𝜇 𝛽 , (11)  \nwhich asymptotes to 3 for 𝑀 ≫ 𝑀 𝑐 , and 0 for 𝑀 ≪ 𝑀 𝑐 . Equation (11) is the updated parameterization from Giri & Schneider (2021) and guarantees positivitity of 𝛽 for all masses. The other power-law indices, 𝛾 and 𝛿 , control the slopes at 𝑟 ∼ 𝑅 ej and 𝑟 ≫ 𝑅 ej, respectively.  \nThe normalization, 𝜌 gas , 0 is set in an analagous manner to the stellar profile normalization of Equation (7), and is given as  \n𝜌 gas , 0 = 𝑓 gas 𝑀 tot GLYPH<20> ∫ ∞ 0 4 𝜋𝑟 2 𝜌 gas ( 𝑟 ) 𝑑𝑟 GLYPH<21> -1 , (12)  \nwith 𝑓 gas derived from the baryon and stellar mass fraction,  \n𝑓 gas = 𝑓 b -𝑓 star , (13)  \n𝑓 b = Ω b / Ω m , (14)  \n𝑓 star = 2 𝐴 GLYPH<20> GLYPH<18> 𝑀 200c 𝑀 1 GLYPH<19> 𝜏 + GLYPH<18> 𝑀 200c 𝑀 1 GLYPH<19> 𝜂 GLYPH<21> -1 , (15)  \nThe baryon fraction within an isolated, toy-model halo is simply the cosmological baryon fraction, and once we specify the stellar fraction - which contains in it the rich physics of star formation - we can obtain the gas fraction of the halo. Note that 𝑓 star denotes the fraction of all stellar matter, including those in satellite galaxies, and we will use this fact in Equation  \n(20). Next, following P24, we add additional mass and redshift dependence to multiple gas profile parameters through the following promotion for each parameter,  \n𝑋 → 𝑋 GLYPH<18> 𝑀 200c 𝑀 𝑋 GLYPH<19> 𝜇 𝑋 ( 1 + 𝑧 ) 𝜈 𝑋 ( 𝑐 200c ) 𝛼 𝑋 . (16)  \nwhere 𝑀 𝑋 is a pivot mass, 𝜇 𝑋 , 𝜈 𝑋 , 𝛼 𝑋 control the scaling with mass, redshift, and halo concentration. Note that in practice, most of these additional scaling parameters are not included in the final model; see Section 2.5 and Table 1 for more details.  \nThe collisionless matter profile is the final component of the baryonification density profile model. This profile constitutes both the DM as well as the subhalos/galaxies in the halo, as both components are collisionless and interact only gravitationally (and not hydrodynamically) with the gas distribution. We follow previous baryonification models in ignoring any explicit modelling of gas within individual subhalos/galaxies. For the observables we wish to model - the tSZ field, and the baryonic imprints in the density field - the gas distribution in substructure is irrelevant compared to the gas of the host halo. The collisionless matter profile is modelled by accounting for the adiabatic, gravitational contraction/expansion of the original matter distribution (Blumenthal et al. 1986; Gnedin et al. 2004) - where this original distribution is taken to be an NFW profile; see Equation (1) - due to the presence of baryonic components, namely gas and stars.  \nThis adiabatic relaxation is performed via the transformation 𝑟 → 𝑟 / 𝜁 , where 𝜁 is obtained by solving  \n𝜁 ( 𝑟 ) -1 = 𝑎 GLYPH<20> GLYPH<18> 𝑀 𝑖 ( 𝑟 ) 𝑀 𝑓 ( 𝑟 ) GLYPH<19> 𝑛 -1 GLYPH<21> , (17)  \nwhere 𝑎 and 𝑛 are phenomenological parameters. For 𝑎 = 𝑛 = 1, the expression corresponds to exact angular momentum conservation. We instead take 𝑎 = 0 . 3 and 𝑛 = 0 . 2 (see Table 1), which provides a better match to simulations (Abadi et al. 2010). The deviation from 𝑎 = 𝑛 = 1 arises because galaxy growth is not an instantaneous process (Gnedin et al. 2004; Abadi et al. 2010).  \nThe initial (final) mass distribution 𝑀 𝑖 ( 𝑀 𝑓 ) is given by,  \n𝑀 𝑖 = 𝑀 nfw ( 𝑟 ) (18)  \n𝑀 𝑓 = 𝑓 CLM 𝑀 nfw ( 𝑟 ) + 𝑀 gas ( 𝑟𝜁 ) + 𝑀 star ( 𝑟𝜁 ) . (19)  \nWesolve for 𝜁 iteratively, and find that in practice it converges ( < 1% differences) in under 10 iterations. Upon obtaining a solution, we compute the final 𝜌 CLM profile as,  \n𝜌 CLM ( 𝑟 ) = 𝑓 CLM 4 𝜋𝑟 2 𝑑 𝑑𝑟 𝑀 nfw ( 𝑟 ) , (20)  \nwith the definition 𝑓 CLM = ( Ω m -Ω 𝑏 )/ Ω m +( 𝑓 star -𝑓 cga ) . The latter two fractions are defined in Equation (9) and (15), and their difference 𝑓 star -𝑓 cga corresponds to the mass fraction of  \nstars in satellite galaxies. This quantity, 𝑓 CLM, corresponds to the fraction of collisionless matter in the total matter distribution.  \nWe have now described all the profile components needed for modelling the baryonic imprints. We can write the total density distribution in the DMO and DMB halo models as  \n𝜌 DMO = 𝜌 NFW + 𝜌 2h (21)  \n𝜌 DMB = 𝜌 CLM + 𝜌 gas + 𝜌 star + 𝜌 2h (22)  \nFigure 1 shows examples of the different profiles we use in this work. We see that the collisionless matter profile (CLM) is amplified at small scales and suppressed at larger scales, when compared to the dark matter NFW profile, as is expected from the process of adiabatic relaxation. The DMB model is further amplified at small scales due to the presence of the stellar component. We have verified that our pipeline reproduces the profiles shown in Figure 1 of S19.', '2.2 Matter density field': "Thenovelty of the baryonification model is in taking results of an N-body simulation - which includes accurate predictions for the gravitationally induced correlations of non-linear structure - and perturbing it to induce baryon-like features in the density field. Such an approach requires a displacement function, Δ 𝑑 ( 𝑟 ) , which specifies the distance a particle must be perturbed from its current location.  \nTo define this function, we first compute the enclosed mass corresponding to the density distributions in Equations (22) and (22). Then, the displacement function is given by  \nΔ 𝑑 ( 𝑟 ) = 𝑀 -1 DMB ( 𝑀 DMO ( 𝑟 )) -𝑟, (23)  \nwhere 𝑀 -1 DMB ( 𝑀 ) is the functional inversion of enclosed massradius function /five.sup and outputs a radius given a mass, and 𝑀 DMO ( 𝑟 ) is the enclosed mass in the DMO case. For a given halo, we evaluate this function at the location of every particle within some maximum halo-centric distance (nominally set to 20 𝑅 200c), and perform a radial offset to the particle position with respect to the halo. In practice, we do not immediately offset the particle locations. Instead, we loop over all halos in the volume, and accumulate the offsets per particle. Once the loop is finished, we now shift the particles by the total, accumulated offset. In this way, the algorithm and its predictions are invariant to changes in the ordering of halos in the loop. The offset positions obey periodic boundary conditions of the simulation box; for example, any pixel that is offset beyond the left edge of the simulation volume will reappear on the right edge of the volume.  \nThus far, we have used the formalism of S19. However, as discussed previously, we require a faster method for applying this displacement function to large suites of full-sky simulations. Fluri et al. (2019) introduced the shell baryonification method, where the displacement is performed on individual map pixels rather than on individual matter particles. In this case, the displacement is written as,  \nΔ 𝑑 ( 𝑟 p ) = 𝑀 -1 DMB , p ( 𝑀 DMO , p ( 𝑟 p )) -𝑟 p , (24)  \nwhere we have substituted the 3D distance 𝑟 with the projected distance 𝑟 p. The enclosed masses are all now masses enclosed within the projected radius,  \n𝑀 𝑋,𝑝 ( 𝑟 p ) = ∫ 𝐿 p / 2 0 2 𝑑𝑙 ∫ 𝑟 p 0 𝑑𝑟 p2 𝜋𝑟𝜌 𝑋 GLYPH<18> √︃ 𝑙 2 + 𝑟 2 p GLYPH<19> , (25)  \nwhere we now integrate over the line-of-sight distance 𝑙 . The factor of 2 in the integral over 𝑑𝑙 arises from the symmetry in 𝑙 used to set integration limits to 0 < 𝑙 < 𝐿 p / 2 instead of -𝐿 p / 2 < 𝑙 < 𝐿 p / 2. The choice of 𝐿 p / 2 as the limit, where 𝐿 p is the thickness of the shell (or simulation box, when considering a 3D snapshot), is different from Fluri et al. (2019), where the integration limits were set to 50 𝑟 p instead. We motivate the former choice by recognizing that the mass in a given map pixel is a line-of-sight integration overing the distance -𝐿 p / 2 < 𝑙 < 𝐿 p / 2. Thus, a similar choice must be made in the baryonification model.  \nThis choice of projection scale can also be motivated by taking the limit 𝐿 p →∞ , where 𝑀 DMB , p ≈ 𝑀 DMO , p since the large-scale matter distribution (which is much less affected by baryons than the small-scale distribution, and is modelled by the two-halo term in our work) dominates the mass in the pixel. In this limit, the displacement function is simply Δ 𝑑 ( 𝑟 p ) ≈ 0. Bysetting 𝐿 p as the integration distance, our model predictions are formally consistent in these extreme limits. Appendix B shows the changes in the predicted displacement due to variations in the assumed projection scale. We note that the projection integral in Equation (25) assumes the halo is spherically symmetric. We test the anisotropic baryonification corrections in Section 4.1.3, but the baryonification model in those cases will still use the integral in Equation (25).  \nSimilar to the case with the particle baryonification, we loop over all halos in a given shell and accumulate offsets to the pixel locations. Once the loop is done, we shift the pixel locations by the sum of the accumulated offsets. For each offset parent pixel, we use the H/e.pc/a.pc/l.pc/p.pc/y.pc interpolation routine get interp weights to (i) determine the nearest four H/e.pc/a.pc/l.pcP/i.pc/x.pc child pixels that the parent pixel contributes to, and (ii) determine the contribution the parent pixel makes to each child pixel, which is determined by the weights provided by the function. Thus, we regrid the offset pixel back onto the H/e.pc/a.pc/l.pcP/i.pc/x.pc grid using the neighboring pixels and corresponding  \nFigure 1. An illustration of the different profiles used in the baryonification pipeline. The top left shows the dark-matter only (DMO) and dark matter baryon (DMB) models, including all the components that constitute each profile. The top right, bottom left, and bottom right, show the fiducial gas, star, and pressure profile, respectively. In addition to the fiducial parameter choice, listed in Table 1, we show a handful of variations (see legends) to illustrate the impact of certain parameters on the profiles.  \n<!-- image -->  \ncontributions as determined by the H/e.pc/a.pc/l.pcP/i.pc/x.pc routines. This procedure is guaranteed to preserve the total mass/density in the shell.  \nIn Section 4, we perform a number of tests to evaluate the model's accuracy and these tests use 2D rectilinear simulation maps rather than curved sky maps. In this case, the baryonification-based mass reassignment is done similar to the H/e.pc/a.pc/l.pcP/i.pc/x.pc case: for each offset parent pixel, we find its nearest four child pixels, calculate the area overlap with each child pixel, and re-assign the parent pixel's values to the four child pixels.  \nWehave verified that our particle baryonification pipeline, operating on the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG DMO simulation, can reproduce the 3D power spectrum results from Figure 2 of S19. We have also checked - using the shell baryonification pipeline - that the angular power spectrum of the baryon-corrected density fields shows the same scale-dependent variation with model parameters as in that Figure; for example, using 𝜃 ej = 2 leads to an excess of power on intermediate scales. /six.sup", '2.3 Thermal Sunyaev-Zeldovich (tSZ) effect': "The tSZ effect is sourced by the thermal pressure of the electrons in the cosmic distribution of gas. Given we already have a model for the total density profile, in Equation (22), and for the gas density profile, in Equation (10), we can now assume the halo is in hydrostatic equilibrium and trivially obtain an estimate of the gas pressure,  \n𝑑𝑃 𝑑𝑟 = -𝐺𝑀 1 -halo DMB ( 𝑟 ) 𝑟 2 𝜌 gas ( 𝑟 ) . (26)  \nThis expression can be integrated to obtain the pressure, 𝑃 = -∫ 𝑟 ∞ 𝑑𝑃 𝑑𝑟 𝑑𝑟 , using the boundary condition of 𝑃 ( 𝑟 →∞) = 0. Note that Equation (26) uses the enclosed mass corresponding to only the one-halo term, denoted above as 𝑀 1 -halo DMB ( 𝑟 ) , and ignores the two-halo contribution. This is so our pressure profile corresponds to only the one-halo component. The tSZ field is generated by pasting tSZ profiles around all halos in the simulation. Under such a method, the two-halo term of the tSZ field is naturally modelled as the clustering of halo positions includes this two-halo component. Thus, we avoid double-counting the two-halo term by explicitly ignoring it in Equation (26).  \nEquation (26) is derived using hydrostatic equilibrium, which assumes that the gravitational pull on the gas is balanced entirely by the thermal pressure - random motions of the ions  \n-of the gas. While this is generally true for the gas closer to the core of the halo, it is not accurate for the gas at the outskirts ( e.g. , Nelson et al. 2014). Thus, Equation (26) only represents the total pressure required to support the gas, and we must estimate the fraction contributed separately by the thermal and the 'non-thermal' pressure, where the latter primarily constitutes turbulent motions of the gas (i.e. the gas' velocity dispersion).  \nWe follow P24 by including in our predictions a parameteric model, based on Shaw et al. (2010), of the non-thermal pressure,  \n𝑓 nt = 𝛼 nt 𝑓 𝑧 GLYPH<18> 𝑟 𝑅 200c GLYPH<19> 𝛾 nt , (27)  \n𝑓 𝑧 = min [( 1 + 𝑧 ) 𝜈 nt , ( 𝑓 max -1 ) tanh ( 𝜈 nt 𝑧 ) + 1 ] , (28)  \n𝑓 max = 6 -𝛾 nt / 𝛼 nt , (29)  \nwhere 𝛼 nt is the amplitude of non-thermal pressure, 𝛾 nt its scaling with radius, and 𝜈 nt its scaling with redshift. We set 𝛾 nt = 0 . 5 and 𝜈 nt = 0 . 3 following P24. The quantities 𝑓 𝑧 and 𝑓 max are defined so 𝑓 nt is between 0 and 1 for 𝑟 < 6 𝑅 200c. We also explicitly clip its value between 0 < 𝑓 nt < 1. The non-thermal pressure is varied in all our predictions. We also show in Section C2 the predictions of this model compared to others in the literature and show the model is adequately flexible.  \nAric'o & Angulo (2024) also introduce an extension of their baryonification method (Aric'o et al. 2021b) to include a tSZ field consistent with the baryon-corrected density field. They first assign 'synthetic' gas and star particles to their simulations using the DM particle dataset. The pressure profile is then modelled using the prescription of Komatsu & Seljak (2002). That prescription connects the gas pressure and density by assuming the gas is a polytrope ( 𝑃 ∝ 𝜌 Γ ), which is allowed under the condition of hydrostatic equilibrium, and then by assuming the total matter distribution follows an NFW-like profile (Komatsu & Seljak 2001, see their Equation 7). Note that Γ is furthermore assumed to be independent of scale, whereas it has been found to be a function of radius in both simulations ( e.g. , Battaglia et al. 2012, see their Figure 4) and observations (Capelo et al. 2012); it does, however, depend explicitly on halo concentration and therefore halo mass. Finally, a particle at a given radial distance, 𝑟 , is assigned a value given by the pressure profile evaluated at that distance.  \nIn comparison, our work predicts the pressure field directly on the sky, avoiding the use of particles to reduce computational cost and relying instead on the halo positions. Our pressure profile is obtained using the assumptions of hydrostatic equilibrium as well, but we do not assume an NFW profile for the total matter distribution. Instead, we directly use 𝑀 dmb ( 𝑟 ) predicted by the baryonification model; see Equation (22). Our model also does not assume the gas is a polytrope or that  \nit has a scale-independent polytropic index. This allows our predictions to be more easily generalized to both the outskirts of halos and to lower-mass halos. Our profile painting method can capture anisotropic features in the field via the use of elliptical profiles (see Section 4.1.3), but it does not capture more fine-grained morphological features, as would be possible in the particle-based method of Aric'o & Angulo (2024).", '2.4 Pixelization': "Our predictions are pixelized maps rather than a discrete, particle/point. There is therefore a window function needed due to this the pixelization, and must be accounted for during the theory modelling. We do so by modifying the profiles as,  \n𝑋 ( 𝑟 ) → GLYPH<20> 𝑊 𝑝 ★ 𝑋 GLYPH<21> ( 𝑟 ) , (30)  \nwhere [ 𝑊 pix ★ 𝑋 ] represents a convolution between the pixel window function, 𝑊 pix, and the radial profile, 𝑋 . This convolution is performed entirely in Fourier space. In the case of 2D maps, which have square pixels, we assume a circular tophat with an area equivalent to that of the square pixel. /seven.sup We have tested that this approximation makes a negligible difference in the convolved profiles: the results are identical for scales smaller/larger than the pixel-scale, and differ at the percent level around the pixel-scale. For the full-sky maps, we approximate the Healpix pixel window function by a Gaussian window function which we verifiy is accurate to better than 1%,  \n𝑊 pix ( ℓ ) = exp GLYPH<20> -𝜎 2 2 ℓ ( 1 + ℓ ) GLYPH<21> , (31)  \nwith 𝜎 = 𝜃 pix /( 4 √ ln 2 ) , 𝜃 pix the pixel resolution scale obtained from the nside2resol function of H/e.pc/a.pc/l.pcP/y.pc, and ℓ is the angular multipole. This analytic function is advantageous as it can be evaluated at any ℓ -which is necessary for performing FFTlog-based Fourier transforms to high numerical precision - whereas the default Healpix window function is only computed/provided up to an ℓ max set by the Healpix map's resolution. /eight.sup  \nFigure 2 shows an example of the pixelization's impact on the density profiles and therefore the displacement function. For larger pixelizations, we smooth over baryonic imprints and so the displacement function approaches zero. It is at zero for 𝐿 = 10 Mpc, which is about 20 times larger than the radius 𝑅 200c = 0 . 45 Mpc for the halo mass we show. There is a pileup of mass (identified as the 'bump' in the profiles) at the scale of pixel, where this scale is shown in vertical dotted lines. The  \nFigure 2. Thedisplacement function (top) and the DMB profile (bottom) corresponding to different pixel/map resolutions (pixel length scale denoted in the legend). The black dotted line is the result from no convolutions, and it is consistent with the result from convolving the predictions with a vanishingly small pixel. The pixel window function smooths the mass on small scales, thereby leading to (i) nearly zero displacement at small scales, and (ii) a pileup of mass at the scale of the pixel (denoted by vertical lines).  \n<!-- image -->  \nblack dotted line in Figure 2 shows the unsmoothed result, and there is numerical consistency between it and the result from using a very small pixel window function. We have verified that the FFTlog calculations, /nine.sup which are used to convert the profiles to Fourier space and perform the convolution, preserve the profiles down to the 10 -5 % level, which is adequate for our modelling. We use these pixel convolutions for all steps in our work. The pixel scales vary depending on the exact simulation/analysis being performed, and we will denote the exact values where relevant.", '2.5 Model parameterization': 'The profile models in Section 2.1 and Section 2.3 has a number of free parameters, but not all of them need to be varied. Table 1 summarizes the parameters of the model, and then the priors we use when performing any model comparisons/fits in the discussions to follow. The key parameters to vary are those associated with the gas profile. Note that this includes 𝜂 , given the stellar fraction directly controls the gas fraction as shown in Equation (13). The rest of the parameters are fixed to values discussed in S19. The parameter 𝜏 , which we introduce in this work in Equation (15), is set to 𝜏 = -1 . 5 which is the 1 𝜎 upper bound from the findings of Moster et al. (2013, see their Table 1).  \nIn Equation (16) we modified the gas profile parameters to include additional dependence on mass, redshift, and halo concentration. We set these additional scalings to zero for most parameters; the notable exceptions are the scale radii, 𝜃 ej and 𝑀 𝑐 , where we follow P24 in varying the relevant mass- and redshift-dependence parameters associated these quantities. We also follow P24 in setting the concentrationdependence parameter 𝛼 𝑋 = 0. The concentration dependence of baryonic feedback has been recently measured in many works ( e.g. , Shao et al. 2023; Shao & Anbajagane 2024; Gebhardt et al. 2023; Pandey et al. 2024) but given our model parameterization is already flexible we do not include this additional dependence.  \nThe halo concentration is formally a free parameter in the model. In practice, for a given halo, we use the halo concentration measured directly in the simulation when available and otherwise assume the simulation-calibrated relation of Diemer & Joyce (2019) - alongside the halo mass, redshift, and assumed cosmology of the simulation - to assign the halo a 𝑐 200c value.', '3 SIMULATIONS & SUMMARY STATISTICS': 'Weuse three distinct simulation suites and one set of summary statistics to validate our model predictions against the simulation measurements. The I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG suite (Section 3.1) and the Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc suite (Section 3.2) are used to validate the baryonification model. The U/l.pc/a.pc/g.pc/a.pc/m.pc suite (Section 3.3) is used to simulate surveys as part of our Fisher forecast. All analyses of the fields are done after summarizing the fields into a data vector containing moments (from 2nd to 4th order) of the lensing and tSZ field (Section 3.4). All maps made from the simulations are either a 2D rectilinear grid of 512 2 pixels or a H/e.pc/a.pc/l.pcP/i.pc/x.pc curved-sky map of NSIDE = 1024.', '3.1 I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulations': "The I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG /one.sup/zero.sup suite (Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018) is a set of high-resolution hydrodynamical simulations that self-consistently evolve a wide set of astrophysical processes (Weinberger et al. 2017; Pillepich et al. 2018a). The simulation volume spans a box with side 𝐿 box = 205 Mpc / ℎ . At such volumes, the particle mass resolution is high but at the cost of poor number statistics for the larger halos. The latter limits our accuracy in validating the pressure field model, as the pressure field depends more heavily on the most massive halos. The suite contains simulations of three resolution levels, and we use the third, TNG300-3, variant. This run has 625 3 particles each of dark matter and  \nTable 1. The input parameters to the profiles of the baryonification model, described in Section 2.1 and 2.3. For each parameter, we list its fiducial value, the prior used when fitting it (left blank if the parameter is fixed in our analysis) and a description of the parameter's physical meaning. The prior ranges are similar to, but generally broader than, the ones defined in P24, see their Table 1. Any mass, redshift, and concentration-dependent parameters - as defined by Equation (16) - are fixed to 0 if they are not shown above.  \ngas, with a DM particle mass resolution of 4 . 5 × 10 9 M ⊙ . While higher resolution runs are also available, we utilize only this low-resolution run as its resolution more closely mimics that of the larger simulations used for modelling full-sky lensing and tSZ fields. The pixel-scale of the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG 2D maps is 0 . 6 Mpc.  \nEach simulation in the suite also has a DMO counterpart, that shares the exact same initial conditions. The pair of simulations can be used to make measurements whose uncertainties are cosmic variance suppressed (see Section 3.5 for details) and therefore allow more precise validations than otherwise possible. We use the public data release, as introduced in Nelson et al. (2019). There also exist other, public simulation suites that span much large volumes than I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG -such as M/a.pc/g.pc/n.pc/e.pc/t.pc/i.pc/c.pc/u.pc/m.pc (Hirschmann et al. 2014) and B/a.pc/h.pc/a.pc/m.pc/a.pc/s.pc (McCarthy et al. 2017) - and would therefore be better suited for our tests of the tSZ field and lensing field cross-correlations. These simulations also span a wide range of baryonic feedback prescriptions, some of which are stronger compared to the relatively weaker feedback in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. However, unlike the case with I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG, the quantities  \nwe require from these simulations (the particle snapshots) are not publicly accessible. For these reasons, we use I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG for our initial validation here, but highlight the need for a more precise validation - using a larger simulation before employing this model on survey data.", '3.2 Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc simulations': 'Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc /one.sup/one.sup is a large suite of DMO simulations spanning a wide range of cosmological parameters (Villaescusa-Navarro et al. 2020). We use the high-resolution runs generated the fiducial cosmology which is based on Planck Collaboration et al. (2016). These run have 1024 3 particles spanning a large volume ( 𝐿 box = 1000 Mpc / ℎ ), and therefore has much poorer mass/spatial resolution than I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG with each DM particle of mass ≈ 10 11 M ⊙ , but also have a statistical sample of halos of higher masses. In our work, we use these simulations for tests that (i) use the halo concentration, as we can utilize the 𝑐 200c estimates from the R/o.pc/c.pc/k.pc/s.pc/t.pc/a.pc/r.pc catalogs (Behroozi et al. 2013) provided with these simulations, and  \n(ii) for checking sensitivities to higher masses that are not testable in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. The pixel-scale of the 2D maps is 3 Mpc.', '3.3 U/l.pc/a.pc/g.pc/a.pc/m.pc simulations': "The U/l.pc/a.pc/g.pc/a.pc/m.pc /one.sup/two.sup suite is a set of full-sky N-body simulations designed for widefield survey analyses, and was introduced in Anbajagane et al. (2024b). The simulations are based on the Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc and Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc-P/n.pc/g.pc suites (Villaescusa-Navarro et al. 2020; Coulton et al. 2022) and complement them by enabling Fisher forecasts of observables from wide-field surveys. In particular, they vary four cosmology parameters, and four inflationary parameters as well. In this work, we use these simulations to perform forecasts for weak lensing and tSZ surveys (Section 5).  \nThe simulations resolve 512 3 particles in a volume with 𝐿 box = 1 Gpc / ℎ , following the same configuration as the fiducial (not high-resolution) Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc run. However, the halo catalogs provide friends-of-friends (FoF) mass, estimated with a linking length of 𝑏 = 0 . 2 in units of the mean interparticle separation, instead of 𝑀 200c. This catalog is used only in our forecast analysis (Section 5), and we make the approximation of treating 𝑀 fof masses as though they were 𝑀 200c ones. FoF linking lengths of 𝑏 = 0 . 2 will generally correspond to halo overdensities of 200 𝜌 𝑐 , but this statement depends on the concentration of the specific halo (More et al. 2011, see their Figure 4). Even with this caveat, we can still meaningfully study the degeneracy directions of the different baryonification parameters, as well as the change in relative constraining power across different analysis configurations.  \nHalos of 𝑀 fof > 10 14 M ⊙ are adequately resolved by the simulation as these halos contain at least 100 particles. /one.sup/three.sup Anbajagane et al. (2024b, see their Figure 12 and 13) show that the halo mass function and halo bias of this catalog agreed with those from higher-resolution simulations. This mass limit is also adequate for studying tSZ auto-correlations: we show in Section 4.3, using the high-resolution Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc suite, the sensitivity of our observable (the moments of the fields) to different mass scales. For the scales of interest to us, the tSZ-only measurements are sensitive to halos above 𝑀 200c > 10 14 . 5 M ⊙ . For the lensing-only measurements, we have already shown in Anbajagane et al. (2024b, see their Figure 10 and 11) that the density and lensing field are accurately predicted in these simulations.  \nModelling the cross-correlation of the lensing and tSZ field through these simulations incurs one caveat. Pandey et al.  \n(2022, see their Figure 2) compute the mass-dependence of the lensing and tSZ real-space, two-point cross-correlation function. For the smallest scales we measure ( 𝜃 ≈ 3 ' ) they show the peak contributing halo mass is 𝑀 200c ≈ 10 14 M ⊙ but halos below that mass-scale also contribute a significant amount to the total signal. Our map will therefore underestimate this signal since we only paint pressure profiles around halos with a minimum mass of 10 14 M ⊙ . Higher-order correlations push the peak contributing halo mass to higher masses, therefore making any missing contribution from lower mass halos less relevant. Regardless, we stress that the absence of lower mass halos in our analysis will reduce the impact of baryons on our simulated signal (relative to their impact when using all halos) and therefore will weaken our Fisher information on baryonic processes. We do not risk quoting artificially tighter bounds due to this caveat. In summary, the U/l.pc/a.pc/g.pc/a.pc/m.pc suite still has sufficient, though not ideal, characteristics for performing our joint analysis of lensing and tSZ fields.", '3.4 Field moments statistic': "The baryonification model has been validated for the density field on both the power spectrum (Schneider et al. 2019; Aric'o et al. 2023), as well as the bispectrum (Aric'o et al. 2021a) and the weak lensing peaks (Lee et al. 2023). It has also been validated for the tSZ power spectrum (Aric'o & Angulo 2024). In this work, we focus on the moments of the field which have been used extensively on weak lensing data ( e.g. , Van Waerbeke et al. 2013; Chang et al. 2018; Peel et al. 2018; Petri et al. 2015; Gatti et al. 2022, 2023, 2024b,a), and are a computationally efficient estimators of higher-order information in the fields.  \nThese moments are computed as,  \n⟨ 𝐴 ( 1 ) 𝐴 ( 2 ) . . . 𝐴 ( 𝑁 ) ⟩( 𝜃 ) = 1 𝑁 pix -1 𝑁 pix ∑︁ 𝑖 = 1 𝐴 ( 1 ) 𝑖 𝐴 ( 2 ) 𝑖 . . . 𝐴 ( 𝑁 ) 𝑖 , (32)  \nfor some set of fields 𝐴 ( 1 ) , 𝐴 ( 2 ) , . . . , 𝐴 ( 𝑁 ) , where all fields are smoothed on some scale, 𝜃 . We follow previous works ( e.g. , Gatti et al. 2020; Anbajagane et al. 2023) in performing the smoothing using a tophat-filter, which is defined in harmonic space as  \n𝐵 ( ℓ ) = 2 𝐽 1 ( ℓ𝜃 ) ℓ𝜃 (33)  \nwhere 𝐽 1 ( 𝑥 ) is a Bessel function of the first order. In this work, we will consider two fields: the density field and the tSZ field. For 𝑁 = 2, the moments capture the same information as a power spectrum, and this has been checked extensively in many recent works ( e.g. , Anbajagane et al. 2023; Gatti et al. 2024b).  \nAn ideal validation of our model would also use other summary statistics - such as the poly-spectra, cumulative  \ndistribution functions, the wavelet harmonics, etc. - but we choose a single statistics, the moments, for simplicity in presenting and interpreting our results. Different statistics will probe somewhat similar types of information, even if the exact definitions vary significantly. Thus, validating our approach for one set of statistics serves as a positive sign for other statistics as well, though an explicit validation must always be performed before employing the baryonification method, with that given statistic, on data. One particular note is that the isotropic shape of the filter in Equation (32) means the moments are an angle-averaged quantity. Correlation functions for 𝑁 > 2 points depend on the exact configuration (or shape) of the points being correlated; for correlations of three points, the shapes are different triangles. Statistics computed on isotropically smoothed fields, such as the moments defined here, average over all these configurations/shapes. Thus, the validation of baryonification for these moments does not guarantee the validation of shape-specific correlations.", '3.5 Estimation of measurement uncertainty and model best fit': "One key validation done in this work is to compare the predictions of the baryonification model to the fields from I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. The robustness of the match is quantified in relation to the uncertainty on the simulation measurement. There are two different contributions to this uncertainty: (i) cosmic/sample variance, with quantifies the intrinsic uncertainty of our measurement due to using only a finite volume/sample of the Universe, and (ii) shot noise, which is the poission uncertainty from using a finite number of discrete objects.  \nIn the case of I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG, our key measurement is the fractional difference between the moments measured on the true fields and those measured on the predicted 'baryonification-based' fields. The latter fields are derived using products from the dark matter-only simulation in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. The DMO and hydrodynamical simulations share the same initial conditions, and therefore our predicted fields have the same initial conditions of the hydrodynamical simulations. In this case, the uncertainty on the fractional difference is cosmic variance-suppressed. Note that it is suppressed , not cancelled , as the late-time density distributions in the two simulations still vary slightly due to differences in the dynamical equations being solved in each (the DMO simulation does not solve any hydrodynamical equations).  \nThe uncertainty is estimated as follows: we define the ratio 𝑅 ≡ 𝑋 dmo / 𝑋 hydro, where 𝑋 is the 𝑁 th -order moment of the density field and/or tSZ field. The uncertainties on 𝑅 are computed through a leave-one-out jackknife resampling of the 2D fields, with 64 patches. In each sampling, we remove a contiguous, square patch of the field and compute 𝑅 using the remaining area. In the case of the density field, this  \nis straightforward as both DMO and hydrodynamical simulations produce density fields. However, this process cannot be repeated as is for the tSZ field as the DMO simulation does not have an associated tSZ field. We instead create a mock field through the following transformation,  \n𝑦 dmo = iFFT GLYPH<18> FFT ( 𝑦 hydro ) × FFT ( 𝛿 dmo ) FFT ( 𝛿 hydro ) GLYPH<19> . (34)  \nEquation (34) uses the tSZ field of the hydrodynamical simulation and modifies it in Fourier space to now correspond to the density field of the DMO simulation. This field is then used in estimating 𝑅 for moments of the tSZ field. This new 'artificial' tSZ field accounts for differences in the exact realizations of the hydrodynamical and DMO simulations, which requires information about fourier-mode phases, and also the differences in the clustering of matter on small scales, which requires information about the fourier-mode amplitudes. Note that there are many other reasonable ways to estimate a 'artificial' tSZ field, such as the combination 𝑦 dmo = iFFT ( FFT ( 𝛿 dmo ) × √︁ 𝑃 𝑦𝑦, hydro / 𝑃 𝛿𝛿, hydro ) . We find these alternatives result in larger covariances for the ratio 𝑅 , andwetherefore use the method described above since a tighter (or smaller) covariance estimate results in a more stringent test of our method.  \nDue to the small volume probed by I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG, not all jackknife patches are representative volumes of the Universe (for example, only two out of 64 patches contains a halo with 𝑀 200c > 10 15 M ⊙ ), and therefore the jackknife provides a biased estimate of the covariance. However, this estimate is still useful in interpreting the amplitude of deviations between the baryonification model and I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG.  \nWe find the best fit baryonification model for the simulation measurement by generating predictions for a large number of points in parameter space and performing a simple minimum𝜒 2 search using the fractional residuals: given as 𝜒 2 = ( 𝑀 / 𝑠 -1 ) 𝑇 C -1 ( 𝑀 / 𝑠 -1 ) , where we compare simulations, 𝑠 , with the model, 𝑀 . During this procedure, we set the off-diagonal terms of the covariance matrix, C , to zero. Our numerical estimate of this matrix is noisy, and given the comparatively long datavectors (120 to 240 datapoints) we can incur numerical artifacts/errors during matrix inversion. Increasing the number of jackknife resamplings will alleviate this issue at the cost of further biasing our covariance estimates. /one.sup/four.sup For the purpose of checking that the baryonification model jointly fits the simulation measurements, this choice still enables a meaningful analysis while avoiding numerical artifacts. However, this highlights that a more precise validation on larger simulations - for example, by comparing the  \nbaryonification predictions to full-sky simulations like M/i.pc/l.pc/l.pc/e.pc/n.pc/i.pc/u.pc/m.pcTNG (Pakmor et al. 2023) or the F/l.pc/a.pc/m.pc/i.pc/n.pc/g.pc/o.pc project (Schaye et al. 2023) - will be needed to more robustly characterize this model.", '4 MODEL VALIDATION': 'A primary focus of this work is validating the baryonification model for the density and tSZ fields up to 4th order, by using auto/cross moments of the field as our summary statistic. Here, we will use the density field as an analog of the lensing field as our validation is done with the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulation for which we cannot create full-sky maps (such maps require simulations, like those of U/l.pc/a.pc/g.pc/a.pc/m.pc suite, which have at least twenty to thirty times larger volumes). Validating the model for the density field will automatically validate it for the lensing field, as the latter derives directly from the former: see Equation (A1). /one.sup/five.sup Wewill use the full-sky approach for our forecasting results in Section 5.  \nWe split our validation steps into two categories: (i) methodology and (ii) model. The former encapsulates all decisions related to how we move the particles or paint the tSZ field (including how we use the simulation data as inputs). The latter includes changes in the profile model that can then be propogated into the rest of the pipeline. For brevity, we discuss only the former case here. The latter analysis, which is detailed in Appendix C, shows that our profile models are adequate in their flexibility/parameterization.', '4.1.1 Comparing Particle and Shell baryonification': "An critical choice in this work is modelling the signatures of baryon astrophysics using just 2D fields and halo catalogs. As mentioned earlier, the original baryonification methods focused explicitly on 3D particle snapshots (S19, Aric'o et al. 2023) though some work has employed a 2D approach for weak lensing analyses (Fluri et al. 2019). In the 3D approach, denoted as 'particle' baryonifiction, one applies the displacement function to the particles in the 3D snapshot. This snapshot can then be compressed into a lightcone, which is the relevant dataproduct for lensing analyses. In the 2D approach, denoted 'shell' baryonification, we apply the displacement function to the 2D pixels. The latter is computationally cheaper and the discretization step in converting particles to maps significantly reduces the memory footprint.  \nThis is an increasingly important consideration given the community is already producing simulation suites nearing O( 10 4 ) individual simulations ( e.g. , Villaescusa-Navarro et al. 2020; Coulton et al. 2022; Kacprzak et al. 2023; Anbajagane et al. 2024b).  \nFor our tSZ profile painting method, we can interchangeably use the 3D and 2D approach as the difference is simply in what step the line-of-sight integration is done. /one.sup/six.sup In the first method, we paint profiles on a 3D grid (we do not use particle information when painting profiles) and then integrate along line-of-sight by collapsing the map along one direction. In the second, the integration is done within the definition of the projected, 2D profile, before any painting happens. Thus, the tSZ field does not suffer accuracy drops between using 3D fields and 2D fields.  \nOn the other hand, the density field predication can vary between the particle and shell baryonification methods. We test this explicitly by performing the two on the 𝑧 = 0 snapshot from the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulation; in this case 'shell' baryonification is done on the 2D rectilinear map and not the full-sky maps. We first perform both types of baryonification for 32 points in the parameter prior space (defined in Table 1), drawn using a Sobol sequence (Sobol 1967). We use only 32 models given the computational expense of the particle baryonification procedure, which takes O( 10 3 ) longer than the shell-based method. We will use a significantly finer sampling (1024 samples) in Section 8 when comparing our shell baryonification predictions to I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG measurements.  \nOnce we have made the map for a given model, we take the fractional difference of the moments measured on the two maps from particle/shell baryonification, and then compute the average absolute ratio over all 32 models. Figure 3 shows this average difference for the moments of the density field. In practice, the offsets are within the 1-2% level, which is a negligible error. The error in the third moments increases at around the pixel-scale of ≈ 1 Mpc, and the same for the fourth moments increases below ≲ 2 Mpc. These differences are still within the cosmic variance-suppressed uncertainties of the measurement (see gray bands in Figure 8). We also show the 1 𝜎 spread of these absolute, fractional differences and find them to also be around the 1 -2% level for the scales mentioned above.  \nIn general, a difference between the two methods can occur in two different ways: (i) the geometry of the offsets are fundamentally different in the particle and shell baryonification methods, as in the former some particles will be displaced just along the line-of-sight and have no impact on the final, 2D  \nFigure 3. The absolute fractional difference in the baryonification predictions, averaged over 32 parameter choices, between the 2D shell method and the 3D particle method (see Section 2.2 for details on both methods). The dashed lines are the standard deviation of the absolute fractional differences. The 32 models are a Sobol sequence drawn from a prior defined in Table 1. The bias is within 2% across all scales and increases for higher-order statistics for scales below ≲ 2 Mpc. We deem this a negligible systematic; see text for more details.  \n<!-- image -->  \nfield, whereas the projected displacements are always perpendicular to the line-of-sight; and (ii) the presence of interlopers, which are halos that overlap in 2D space but are distant in 3D space, as the offsets predicted by the smaller halo will impact the matter distribution of the larger halo given both distributions will contribute to the pixel being displaced. Formally, this latter effect is indeed accounted for through the two halo profile in Equation (3). Both effects are alleviated when the projection is done over smaller distances, i.e. if the shells are thin. The tests above are done with the projection scale set to 𝐿 = ≈ 300 Mpc whereas in full-sky simulations - such as the U/l.pc/a.pc/g.pc/a.pc/m.pc we perform shell baryonification on in Section 5 the shell thickness is ≈ 60 -100 Mpc. Therefore, Figure 3 is a more strict test of the method.  \nFollowing the discussions in Section 2.2, the model for the displacement function and projected displacement function should be equivalent as both are derived from the same input profile and are physical models of the 3D and projected density distribution, respectively. Figure 3 shows our predictions obey this condition down to at least ≈ 1 Mpc (2 Mpc) for up to thirdorder (fourth-order) in the density field.", '4.1.2 Dependence on concentration': "To first order, halos can be defined using just their mass and redshift. Once these two properties are given, many other halo properties can be approximately derived through simple scaling arguments. However, there are a number of secondary halo properties that add significant information about the halo. Of these, the halo concentration plays a key role and has been studied extensively as a result.  \nIn our current baryonification model, which does re-  \nFigure 4. The impact of concentration on the displacement function for halos of two different masses. There is a strongly scale-dependent effect from varying 𝑐 200c . Note that we use a large range of 𝑐 200c values (beyond our physical expectations for halos of these masses) for illustrative purposes. In general, lower 𝑐 200c halos require stronger contraction (more negative displacements) on small scales, and the physical scale of maximum displacement is also shifted to slightly larger scales. This prediction is made using the fiducial values of the baryonification parameters (see Table 1).  \n<!-- image -->  \nquire an input 𝑐 200c per halo, we simply use a precomputed concentration-mass relation - in our case the model of Diemer & Joyce (2019) - and assign each halo a 𝑐 200c. This is then used in the rest of the pipeline. This is effectively a mass-only model as the halo is still defined using just its mass; a halo of a given mass will also have a fixed concentration value. However, at fixed halo mass there is a ∼ 40% scatter in the concentration values ( e.g. , Wechsler et al. 2002; Diemer & Kravtsov 2015; Diemer & Joyce 2019; Ragagnin et al. 2019; Anbajagane et al. 2022a). Furthermore, the concentration value of a halo depends on the halo's cosmological context, i.e. its entire accretion history, its environment and more ( e.g. , Wechsler et al. 2002; Mansfield & Kravtsov 2020). Thus, the baryonification model's predictions for the cosmological correlations can be more accurate if we self-consistently include the concentration of each halo as measured in the simulation.  \nFirst, we show in Figure 4 the impact of 𝑐 200c on the displacement function for halos of the Milky Way mass-scale and the cluster mass-scale. In both cases, lower 𝑐 200c values lead to a more negative displacement at small scales. A lower 𝑐 200c implies a more diffuse halo and therby requires a stronger displacement value to transform the mass distribution into the more concentrated form that is generated by the presence of baryonic matter. The second effect of a lower 𝑐 200c is  \nFigure 5. The fractional difference in the baryonification predictions, for the 2nd, 3rd and 4th moments of the density and tSZ field, when including/ignoring the concentration information. The difference is below 1% for almost all correlations, and is most amplified only for the cross-correlations where varying the concentration will generate more spatially correlated features across the density and tSZ field.  \n<!-- image -->  \nincreasing the radial scale corresponding to the maximum displacement. This is because for a more diffuse halo, the gas distribution starts dominating the dark matter profile at larger radii. This analysis uses the fiducial baryonification parameters as defined in Table 1.  \nWe then test the impact of this difference on the statistics of interest. Note that this test is done on Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc given those simulations have a readily available concentration measurement, and therefore the range of scales we probe here is larger than that of the other tests. We generate 32 models across the parameter prior and show the median, fractional deviation in the statistics. Figure 5 shows there is a ≈ 2%offset in the measured moments down to 5 Mpc, which is the smallest scale we can probe in Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc. In summary, while the 𝑐 200c variation does have an impact on the statistics, it is not significant above > 5 Mpc and is completely negligible above > 10 Mpc. It is likely the effect is more pronounced towards smaller scales and for higher orders of the fields. We find it is best to use the simulation-measured concentration when available.  \nFigure 4 also shows that the difference in the displacement function is minor when varying the concentration by nearly two orders of magnitude. This implies the model is insensitive to using other halo concentration-mass relations ( e.g. , Child et al. 2018; Ishiyama et al. 2021; Anbajagane et al. 2022a; Shao et al. 2023; Sorini et al. 2024) which deviate, to varying  \nlevels, from the predictions of Diemer & Joyce (2019).  \nOne can also trivially amplify the impact of concentration on the measurements by setting 𝛼 𝑋 ≠ 0 (see Equation (16) and Section 2.1) to enforce a relationship between 𝑐 200c and the gas profile parameters. Thus our test above is primarily showing that the default, fiducial model is not sensitive to the choice of concentration. While P24 find values of 𝛼 𝜃 ej that allow the model predictions to match the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulations, more detailed simulation studies are required to understand the correlation between the gas parameters and the halo concentration. Such a correlation could be important for the baryonification technique as the halos' concentration values are correlated with their spatial distribution, an effect commonly referred to as 'assembly bias' ( e.g. , Zentner et al. 2014). Note that some aspect of this phenomenon can already be included in the current baryonification method, through using the concentration measured for each individual halo in the simulation (as opposed to assuming a concentration-mass relation, which will not include the effects of assembly bias).", '4.1.3 Dependence on ellipticity': "Another critical secondary property of the halo is its shape. To good approximation, halos are quasi-spherical objects. However, they are not exactly spherical, and their ellipticities/orientations are sourced by the cosmological process of structure formation; meaning the shape and orientation of halos exhibit cosmologically sourced, spatial correlations. Such correlations are closely connected to the 'intrinsic alignments' effect studied in the weak lensing community (see Troxel & Ishak 2015, for a review).  \nIt is therefore important to determine the impact of the ellipticity correlations in the baryonification pipeline, for modelling the statistics of interest to us. The implementation of ellipticities in our model is trivial. The halo-centric distance to a pixel is now scaled as 𝑥 2 + 𝑦 2 → 𝑥 2 + 𝑦 2 / 𝑞 2 where 𝑞 ≤ 1 is the ratio of the minor axis to the major axis. /one.sup/seven.sup Under this formalism, the major axis is assumed to be of length 𝑅 200c. In practice, this scaling is done by transforming the coordinate system - rotating the pixel grid so the major/minor axis of the halo is along the x-axis and y-axis of the grid - and then rescaling the 𝑦 distance by 𝑞 as denoted above. The relevant quantity, a displacement or gas pressure, is computed at the location given by this rescaled distance. We perform this test with I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG, where we computed the orientations of halos using the same pipeline as Anbajagane et al. (2022a); this work computes the mass inertia tensor in 3D to compute the 3D ellipticity. In our case, we take only the 2 × 2 subset of the tensor corresponding to the x and y directions, and  \nFigure 6. Different maps of the gas pressure: from the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG hydrodynamic run (top left), from the baryonification model of the DMO run assuming circular profiles (top right), elliptical profiles (bottom left), and 'extreme' elliptical profiles (bottom right). The former elliptical model uses ellipticies computed directly in the simulation, whereas the latter is simply a scaled, more-elliptical version of the former, 𝑞 ex = 0 . 6 𝑞 TNG . We can clearly see the correlation between the halo orientations and the large-scale density field: the alignments are along filaments or point in the direction of the nodes.  \n<!-- image -->  \ncompute a 2D ellipticity using the new tensor.  \nFigure 6 shows a visual comparison of theoretical maps made with and without the orientation information. One can clearly see the correlation of the ellipticities with the largescale density field. The alignments are often in the direction of filaments or are pointing towards nodes in the field. This is most prevalent in 'Ellip. Ex.' map, which is an extreme version of the ellipticity, given as 𝑞 ex = 0 . 6 𝑞 TNG, and shown purely for illustrative purposes.  \nFigure 7 then shows the change in the moments of the fields due to the inclusion/exclusion of the halo ellipticity in our modelling. The change in the statistics is below percentlevel in almost all cases. In summary, the ellipticity of halos is not necessary for modelling the moments of the field. Other statistics - particularly those with anisotropic information, such as the bispectrum which is a function of angle, and not just isotropic information like the moments - may  \nFigure 7. The fractional difference in the baryonification model predictions when excluding/including ellipticity in the model. We show all moments of the density and tSZ field, up to the 4th moment. The impact is sub-percent for all moments, and generally grows towards small-scales.  \n<!-- image -->  \nhave stronger sensitivity to this choice and therefore must be separately tested. Similarly, methods that involve the explicit detection/characterization of filamentary structures, such as topological measures ( e.g. , Heydenreich et al. 2021, 2022; Euclid Collaboration et al. 2023), could have a stronger sensitivity to the choice of including/ignoring ellipticity.", '4.2 Comparison to hydro simulation': "Wenowcheckthatour baryonification model can jointly fit measurements from hydrodynamical simulations. This is done by computing the fractional difference between the simulation measurements and our best-fit baryonification predictions for all moments of the density and tSZ field from 2nd to 4th order - and comparing the difference against the measurement uncertainty described in Section 3.5.  \nFigure 8 compares our predictions with measurements made on I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. Our predictions use the DMO simulation whereas the measurements are made on the full hydrodynamical simulation. Each subplot shows the fractional difference for different moments. The gray-bands denote the 1 , 2 , 3 𝜎 uncertainty bounds computed using the method in Section 3.5. The baryonification model is accurate to within the cosmic variance-suppressed uncertainties of the moments; the fraction differences are always within the 1-2 𝜎 bounds (gray bands). Note that the percent-level agreement for the large-scales of the density auto-correlations (at any order) is expected as the baryonic imprints are negligible on such scales.  \nFigure 8. The fractional difference between the moments of the baryonified DMO simulations and of the TNG hydrodynamical simulation. All moments are computed on 2D maps at 𝑧 = 0. The gray bands show the 1, 2, 3 𝜎 uncertainties. The dotted black lines denote 10% fractional errors. Cosmic variance is highly suppressed given both maps derive from the same initial conditions. The joint fit (yellow line) is within 1-2 𝜎 for all different moments. We a-priori expect some irreducible deviations due to the low statistics of massive halos in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG; see Section 4.2 for details. For the panel with ⟨ 𝑦𝑦 ⟩ , the 'Individual Fit' line is underneath the 'Joint Fit' line.  \n<!-- image -->  \nWe a-priori expect to poorly match some features on small-scales, and particularly for any moments that involves the tSZ field. The signal for the tSZ field comes from the most massive clusters, but the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulation has only a handful of such clusters ( e.g. , only two clusters above 𝑀 200c > 10 15 M ⊙ ). For a larger simulation, with a larger sample of clusters contributing to the signal, any individual characteristics of the cluster would average out over the population, and the baryonification prediction - which by construction represents only the average cluster - would match the simulation measurement better. When the simulation contains only a few clusters, these individual characteristics can have a more dominant role in the emergent signal. For example, the most massive cluster in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG - which will dominate the total tSZ signal, especially for the higher-order moments (see Figure 9) - has undergone a recent merger. Our model will not account for this characteristic. However, a significantly larger simulation would have enough massive systems such that the individual merger history of a single object would not have a notable impact on the observed field.  \nFigure 8 also shows the results from fitting only the density  \nfield, or only the tSZ field. The density field is not a sensitive probe of all the different baryonic physics, and therefore a fit to the density alone does not correctly estimate the amplitude of the tSZ field. We also show the result from fitting a single moment at a time. Unsurprisingly, this fit is the best of all, given the large degrees of freedom of the model compared to number of datapoints used to constrain it. Note that in some moments, even the invidiual fit is in the 2 𝜎 regime, which indicates the presence of the individualistic features in the few most massive clusters (which are not captured in the model). We have separately used the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pc simulations (Vogelsberger et al. 2014), a predecessor to the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG suite we use here, and verified that our model predictions are not systematically low on small-scales in those analyses. This further indicates that the offsets seen in Figure 8 is a stochastic fluctuation, induced by the individual histories of the handful of the largest clusters in the simulation.  \nUnder the current statistical uncertainties, we find that our model is able to jointly fit the simulation data at a fractional accuracy that is at worst between 5% to ≈ 50%, depending on the exact moment being considered, with the fractional  \ndifference increasing for moments with larger measurement uncertainties (for example, ⟨ 𝑦 3 ⟩ or ⟨ 𝑦 4 ⟩ , whose signal is dominated by just a few halos in the TNG volume). Given the caveats related to TNG that we have discussed throughout the text above, a more robust analysis is needed to calibrate the baryonification model accuracy to much better precision. Such an analysis would require a much larger simulation; for example B/a.pc/h.pc/a.pc/m.pc/a.pc/s.pc (McCarthy et al. 2017) or M/i.pc/l.pc/l.pc/e.pc/n.pc/i.pc/u.pc/m.pc-T/n.pc/g.pc (Pakmor et al. 2023). We leave these efforts to future work.  \nWhile the validation of the baryonification method for the density two-point statistics (or 2nd moments) has been done to the percent-level ( e.g. , Giri & Schneider 2021; Aric'o et al. 2021b), the model for the other eleven measurements shown in Figure 8 can still be useful even if the model accuracy is larger than 1%: the noise level in the higher-order moments increases significantly compared to the 2nd moments ( e.g. , Anbajagane et al. 2023, see their Figure 6, for the lensing-only moments), while the tSZ field is generally more noise-dominated than the lensing field. These two factors add leniency to the validation requirements of these other eleven moments. The actual requirements for using baryonifcation is survey-specific and must be determined explicitly for each survey and choice of statistic.  \nAppendix D also shows that the baryonification model can jointly fit all moments measured across multiple redshifts in the I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulations.", '4.3 Mass dependence': 'In Section 4.2 we show that our model is flexible enough to jointly fit the simulation measurements within their statistical uncertainties. We now use this model to check the massdependence of the predicted signal for the different moments. This is analogous to similar plots made for the two-point function in Pandey et al. (2022, see their Figure 2) or To et al. (2024, see their Figure 3), and lend insight into the origin of the signal in different measurements. In general, we expect higher-order moments to be more sensitive to more massive halos than the 2nd moments, and similarly, for the tSZ field to be more sensitive to more massive halos than the lensing field. The former is because taking higher powers of a variable amplify the tails of its distribution, and the density field has a strongly positive-skewed tail - occupied by massive halos -whose presence/importance is amplified by taking higher powers of the field. The latter is because weak lensing depends linearly on the density field and therefore mass, whereas the tSZ signal follows ∝ 𝑀 5 / 3 ( e.g. , Lee et al. 2022, see their Figure 1 for a comparison across simulations) and therefore is sensitive to higher masses than lensing.  \nWe estimate the mass-dependence of the moments using the Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc simulations, and in particular using their high-  \nresolution runs, where halos of 10 13 M ⊙ are resolved by 100 particles. The contributions from halos of a given mass, 𝑀 1, is computed through the difference 𝑋 ( > log 10 𝑀 1 ) -𝑋 ( > log 10 𝑀 2 ≡ log 10 𝑀 1 -𝛿 ) , where 𝛿 = 0 . 25 and 𝑋 is a Nthorder moment measured on a density and/or tSZ field, with the field generated using only halos of mass greater than 𝑀 1 or 𝑀 2. We use the baryonification model to make fields for different choices of 𝑀 and take simple numerical differences to estimate the derivative of the measurement with mass, d ln 𝑋 / d ln 𝑀 . We approximate d ln 𝑋 ≈ Δ 𝑋 / 𝑋 fid where 𝑋 fid is the fiducial signal from using all halos. Thus d ln 𝑋 / d ln 𝑀 represents the fractional change - compared to the fiducial case of using all available halos - in the moments of the fields due to halos of a given mass.  \nFor the tSZ field, the change d ln 𝑋 / d ln 𝑀 represents the contribution of a given halo mass to the total tSZ signal, but for the density field it represents the contribution to the baryonic imprints in the field and not the density field as a whole. Note that for the density field, the 𝑋 fid used in the expression Δ 𝑋 / 𝑋 fid is the total density field signal (and not just the total baryonic imprints signal). For this reason, the fractional differences for the density field will have a lower amplitude than those for the tSZ field. For all estimates, we have used 10 high-resolution simulations and averaged the derivatives across realizations. The mass-dependence will change based on the chosen baryonification parameters; for example, changing the 𝜇 𝑋 parameters in Table 1 can modify these derivatives. For our estimate of the derivative, we will use the fiducial parameter values listed in Table 1.  \nFigure 9 shows these derivatives across different halo masses, for five different apertures used in computing the moments. As expected, the ⟨ 𝛿𝛿 ⟩ correlation has the lowest mass-dependence, as all other correlations have derivatives that peak towards higher masses. The derivatives of the density auto-correlations are also negative as the presence of halos (and namely, the astrophysical processes within them) generally suppresses power on these scales, while for the tSZ field the derivative will be positive as the presence of these halos increases the signal. The derivatives of moments that include the tSZ field are orders-of-magnitude larger than those that include the density field alone, showcasing the dramatically higher sensitivity of the tSZ field to baryonic physics. As discussed previously, the mass-dependence of a measurement peaks at higher masses when we increase the order of the moment or if we include the tSZ field.  \nThe mass-dependence in Figure 9 highlights the difficulty in validating our method using simulations that span only sub- Gpc volumes. If we focus on the 2nd moments of the tSZ field, the minimum contributing mass is 𝑀 200c ≈ 10 14 . 5 M ⊙ , whereas for the 4th moments this mass scale is 𝑀 200c ≈ 10 15 . 2 M ⊙ . The I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG simulations have  \nFigure 9. The derivative of the baryonification signal, with respect to the halo mass, for the moments of the density field and tSZ field (different panels, 𝑋 = ⟨ 𝛿𝛿 ⟩ , ⟨ 𝛿𝑦 ⟩ , . . . ) measured at different smoothing scales (colors). The y-axis represents the fractional contribution, from halos of a given mass 𝑀 , to the total signal. For the density field, the signal is the baryonic imprints in the field rather than the density field itself. See text for details on the measurement. The mass-dependence of higher-order moments, or any moments including the tSZ field, is increased to higher masses when compared to the 2nd moments of the density field. The exact values of the derivatives depends on the baryonification parameters used, and we use the fiducial values shown in Table 1.  \n<!-- image -->  \n/circledot  \n/circledot  \n/circledot  \n/circledot  \naround O( 100 ) halos at the former mass scale, but only a single halo at the latter mass scale. We require hydrodynamical simulations of Gpc-scale volumes in order to validate these higher-order moments to good accuracy. However, we reiterate that in practical use-cases, the accuracy requirements for these higher-order moments - especially those of the tSZ field - will be more lenient given the larger measurement uncertainties of such moments.', '5 FORECASTING': "Wenowforecast the power of current and future surveys in constraining these baryonification parameters, and study the parameters' different degeneracy directions. This is done by forward modelling the observables from surveys; namely, the lensing convergence field and the thermal Sunyaev-Zeldovich field. For brevity, we only briefly describe the procedure below and provide a detailed discussion in Appendix A.  \nThe lensing field is modelled using the same pipeline from Anbajagane et al. (2023, 2024b). The maps are generated from the density fields of the U/l.pc/a.pc/g.pc/a.pc/m.pc simulations (Section 3.3), and then processed to contain all the observational effects found in the data, including the relevant survey mask/area and source galaxy redshift distribution. The two lensing surveys we focus on are: (i) the Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005), which is an optical imaging survey of 5,000 deg 2 of the southern sky, and is cur-  \nrently the largest precision photometric dataset for cosmology. The Year 3 data products and cosmology results are available (Sevilla-Noarbe et al. 2021; Abbott et al. 2022), while the legacy Year 6 dataset is not yet available at the time of publishing this work. We focus on this Year 6 dataset. We then consider (ii) the Rubin Observatory Legacy Survey of Space and Time (LSST), which is a 14,000 deg 2 survey that probes higher redshifts, and is the successor to current weak lensing surveys. We forecast only for the final, year 10 dataset. The characteristics of the surveys (shape noise, and 𝑛 ( 𝑧 ) distribution parameters) are given in Table A1. More details can be found in Appendix A1  \nNext, the tSZ field is modelled by following the approach of Raghunathan et al. (2022a,b). We include simple Gaussian models for a number of 'foreground' componenets, such as the cosmic infrared background, radio background etc., as well as the thermal and atmospheric noise. We focus on two CMB surveys: (i) SPT-3G (Benson et al. 2014) is a survey of the southern sky with considerably deeper (lower noise) observations than other CMB surveys. The entire footprint covers O( 10 4 ) deg 2 of the sky with different noise levels. We focus on the main 1500 deg 2 footprint (i.e. excluding any wide/summer fields) that also overlaps completely with DES Y6. (ii) The Simons Observatory (SO, Ade et al. 2019) is an upcoming survey of ≈ 20 , 000 deg 2 of the sky, and has nearly complete overlap with the LSST Y10 footprint. For both surveys, we assume the analysis is supplemented by Planck data like in ex-  \nting surveys (Madhavacheril et al. 2020; Bleem et al. 2022). The characteristics of the surveys (noise properties, beam, frequency bands etc.) are given in Table A2. Further details are provided in Appendix A2.  \nWe note that all our analyses use relatively 'medium'resolution data-products: the U/l.pc/a.pc/g.pc/a.pc/m.pc simulations are run with 512 3 particles, and produce maps of NSIDE = 1024. These choices set a limit on the scales we can analyse. Simulations of higher resolution enable access to even smaller scales, which would place tighter constraints on the baryonic processes discussed in this work. Therefore all our constraints below could be improved further by pushing the analysis to even smaller scales. In this work, we are limited by the resolution of the simulations available for a Fisher forecast.  \nThese results should also not be over-interpreted or generalized since our analysis explicitly limits the tSZ measurements to 𝜃 > 3 ' , whereas the tSZ power spectrum has the highest signal-to-noise at even smaller scales, ℓ ∼ 3000 (Raghunathan et al. 2022b, see their Figure 4). Our tSZ maps are generated at NSIDE = 1024 - the same choice as that of the lensing field - whereas the maps from current SPT datasets are made at NSIDE = 8192 (Bleem et al. 2022). Thus, there is more usable information in the tSZ field that we have not accessed in our datavectors given we limit our maps (and therefore our analysis) to a coarser resolution. This is done primarily in the interest of minimizing the computation cost, as generating a tSZ map of NSIDE = 8192 would take around 64 × longer than generating one of NSIDE = 1024. Note that for the lensing analysis, the choice of NSIDE = 1024 is realistic, though it may be possible to access smaller scales ( NSIDE = 2048 , 4096) - if the lensing systematics are better controlled/modelled - which would improve our constraints. Theconstraints from both fields could also improve/deteriorate depending on the summary statistics being considered; for examples in the lensing field, see Euclid Collaboration et al. (2023); Anbajagane et al. (2023); Gatti et al. (2024a). We have not explored this direction, and have limited ourselves to the moments in this work, as they have already been used to constrain cosmology (Gatti et al. 2022, 2024a).  \nWe now present the Fisher information for a number of cases, each changing the statistics being measured, the fields being used, and/or the parameters being varied. The key goal of these estimates is to understand the degeneracy directions of the different parameters under different analysis choices, and also the relative constraining power of surveys. The Fisher information is estimated as,  \n𝑭 𝑖 𝑗 = ∑︁ 𝑚,𝑛 𝑑 e 𝑋 𝑚 𝑑𝜃 𝑖 GLYPH<0> C -1 GLYPH<1> 𝑚𝑛 𝑑 e 𝑋 𝑛 𝑑𝜃 𝑗 , (35)  \nwhere 𝑑 e 𝑋 𝑚 𝑑𝜃 𝑖 is the mean derivative of point 𝑚 in data vector 𝑋 with respect to parameter 𝜃 𝑖 , where the mean is computed us-  \nns of the surveys for the analysis of DES and SPT, and 300 realizations for that of LSST and SO. We compute this derivative around the fiducial parameter values, as detailed in Table 1. Then, C -1 is the inverse of the numerically estimated covariance matrix and includes the Hartlap correction factor (Hartlap et al. 2007),  \nC -1 → 𝑁 sims -𝑁 data -2 𝑁 sims -1 C -1 . (36)  \nThe Hartlap factor for all analyses in this work is ≳ 0 . 9. The datavector, 𝑋 , consists of the moments of the field, computed on 10 scales ranging from 3 . 2 ' < 𝜃 < 200 ' . These are the same choices as previous works on the moments of the lensing field (Gatti et al. 2020; Anbajagane et al. 2023; Gatti et al. 2024a) and we adopt it for the tSZ field as well. The covariance matrix is estimated by generating 6000 to 8000 realizations of the surveys, with the number varying depending on the survey being analysed. The O( 10 4 ) realizations are enabled by the 2000 independent simulations available in the U/l.pc/a.pc/g.pc/a.pc/m.pc simulation suite. These simulations are run at the fiducial cosmology, chosen to be the best fit values of Planck Collaboration et al. (2016), and the derived baryonification products use the fiducial values of the model parameters (see Table 1).  \nThe Fisher information estimated in Equation (35) can be artificially increased by numerical noise as this noise can break parameter degeneracies; see Coulton et al. (2023, their Appendix A) for a discussion of this. We quantify this amplification in our estimates by varying the number of realizations used to compute the derivatives and the covariance. The Fisher information changes by ≲ 2% if we estimate the covariance using half the number of realizations. Doing the same exercise on the derivative estimates change the Fisher information by ≲ 10%. Our discussion below focus on constraints improving by more than factors of 2, and are therefore unaffected by any such numerical uncertainties at the 10% level.  \nFigure 10 shows the constraints for different combinations of moments, as measured on mock maps of DES Y6 and SPT 3G; the numbers are listed in Table 2 as well. In addition to the baryonification parameters, we also include the cosmological parameters Ω m, which is the matter energy fraction at 𝑧 = 0, and 𝜎 8, which is the root-mean squared fluctuations of the density field at 𝑧 = 0, smoothed on a 8 Mpc / ℎ scale. We include, and therefore marginalize over, all baryonification parameters that are listed with a prior in Table 1. We however only show a subset of them in this plot for visibility.  \nThe inclusion of higher-order moments (blue and yellow contours) greatly improves the constraints over using the 2nd momentsalone(purple contours). This is expected as baryonic imprints are tied to the astrophysical processes within halos and are therefore stronger on non-linear scales. These scales are also where the field is more non-Gaussian, and parameter constraints improve from including measurements of any  \nFigure 10. Fisher constraints on the baryonification/cosmology parameters (see Table 1 and Section 2) from moments of the lensing and tSZ fields constructed with DES Y6 and SPT 3G specifications, respectively. The contours show 1 𝜎 and 2 𝜎 constraints from using the 2nd/3rd/4th moments. The dashed green/peach lines along the diagonal panels show the 1D posteriors from using only the 3rd/4th moments, respectively. We show only a subset of parameters for the sake of brevity; the full constraints are listed in Table 1. The inclusion of higher-order information significantly increases constraining power.  \n<!-- image -->  \nhigher-order correlations at these scales. Aric'o et al. (2021a) already show that a range of baryonification parameter choices can accurately fit a matter power spectrum from simulations but only a subset of those can also jointly fit the matter bispectrum, indicating that the addition of higher-order information constrains the parameter space further. We expand on this by using up to 4th order in the fields, and also by using both the lensing and tSZ fields. Figure 10 also shows the marginalized constraints from using only the 3rd moments (4th moments), presented as a green (peach) dashed line in the 1D posteriors  \nalong the diagonal panels of the triangle plot. /one.sup/eight.sup In many cases, the 'only 3rd moments' and 'only 4th moments' constraints are more constraining than the 2nd moments alone (though this is not the case for the cosmology parameters).  \nFor most parameters, the largest relative improvement is between using 2nd moments alone and using the combination of 2nd and 3rd moments. Including the 4th moments to this  \nTable 2. The constraints presented in Figure 10, 11, and 12 for (i) different statistics measured on both the lensing and tSZ field (left columns), (ii) for different combinations of fields using the 2nd and 3rd moments (middle columns), and (iii) for different combinations of surveys using the 2nd and 3rd moments of the lensing and tSZ fields (right columns), respectively. The inclusion of higher-order statistics improves all constraints by factors of 3 to 4. Including the tSZ field in the analysis similarly improves constraints by up to factors of 10. The constraints from LSST Y10 and SO are factors of 4-5 better than that of DES Y6 and SPT-3G. The parameters are partitioned visually into cosmological parameters, baryonification parameters shown in the figures of this work, and remaining parameters that are not presented but are marginalized over. In analysis (ii), we do not use 𝛼 nth due to it having no impact on the 𝜅 field; see text in Section 5 for details.  \ncombination improves constraints as well but the relative gain is smaller. This is similar to the lensing-only analysis of Anbajagane et al. (2023), where the 4th and 5th moments did not improve the constraining power. The one exception we find is 𝜃 ej, which improves by nearly a factor of three when including the 4th moment.  \nFigure 10 indicates that a number of degeneracy directions are broken by including higher-order information. It also shows that the cosmological parameters, Ω m and 𝜎 8, can be constrained even when marginalizing over astrophysical systematics. Note that these analyses use both the lensing field and the tSZ field. Limiting the analysis to the lensing field alone deteriorates the constraining power significantly (see Figure 11 below). In general, this result highlights the utility in jointly constraining cosmology and astrophysics using a combination of the lensing field and tSZ field. P24 have already shown this to be the case when using only two-point statistics in the data vector (see their Figure 8).  \nWhile the above analysis decomposes the information into different orders, we now decompose it by field. In Figure 11 we split the Fisher information by convergence field and tSZ field. In the former case, we use all auto/cross-correlations between all tomographic bins of the lensing data, while in the latter we use only auto-correlations of the singular tSZ field. We then show the constraints from moments that have at least one power of the lensing field and one of the tSZ field, and finally the constraints from all moments of the two fields. While the above analysis also used the 4th order moments, we will now use the combination of the 2nd and 3rd as our fiducial datavector, given this measurement has already been used (for the lensing field) to constrain cosmology ( e.g. , Gatti et al. 2022, 2024a). This analysis also does not vary the  \nFigure 11. Fisher constraints on the baryonification and cosmological parameters from different combinations of the convergence and tSZ field as measured by DES Y6 and SPT 3G. The datavector uses the 2nd and 3rd moments. The convergence field and tSZ field have different sensitivities to the parameters, and the degeneracy directions are often orthogonal. As a result the combination of fields is significantly more constraining. The 1D marginal distributions of Ω 𝑚 and 𝜎 8 also show a dotted blue line, which is the lensing-only constraints without marginalizing over any baryonification parameters.  \n<!-- image -->  \nµ  \nnon-thermal pressure parameter 𝛼 nt, see Equation (27), since this parameter has no impact on the 𝜅 field, leading to zero Fisher information and therefore a singular Fisher matrix. We  \nchoose to remove this parameter for three variants (lensingonly, tSZ-only, cross correlations-only) in this analysis - and not just for the 𝜅 -only analysis where the issue of a singular Fisher matrix arises - to preserve consistency across analysis setups and therefore simplify the comparisons of the different constraints.  \nFigure 11 shows that the lensing and tSZ field have different sensitivities to combinations of parameters; in many panels, the degeneracy directions of the convergence-only and tSZonly case are nearly orthogonal. As a result, any analysis using combinations of the convergence and tSZ field is far more constraining than the individual fields. Another interesting note is that the lensing field is able to better constrain combinations of the redshift parameters, 𝜈 𝜃 ej and 𝜈 𝑀 c , as shown in Table 1. This is because the lensing survey has more redshift information via the availability of tomographic bins, whereas the tSZ field is a singular field integrated over all redshifts. This result highlights the highly complementary nature of combining the lensing and tSZ fields when constraining baryonification models; P24 (see their Figure 8) have shown the same for the two-point functions. It is also interesting that the crosscorrelations alone can provide a relatively good constraint on the baryonification parameters (as well as cosmology). Such cross-correlations are more robust measurements than their auto-correlation counterparts, as additive systematics in each field (that are uncorrelated across the different fields) will not be present in the measurement.  \nThe field-by-field comparison of Figure 11 also highlights the improvement in cosmology constraints ( Ω m and 𝜎 8) due to the addition of the tSZ field. The cross-correlation of the lensing and tSZ fields is a significant factor in this improvement. The results also indicate that the tSZ-only analysis is unable to jointly constrain baryonification models alongside cosmology. However, this result should not be over-interpreted or generalized since our analysis explicitly limits the tSZ measurements to 𝜃 > 3 ' , whereas the tSZ power spectrum has the highest signal-to-noise at even smaller scales, ℓ ∼ 3000 (Raghunathan et al. 2022b, see their Figure 4). Our tSZ maps are generated at NSIDE = 1024 - the same choice as that of the lensing field - whereas the maps from current SPT datasets are made at NSIDE = 8192 (Bleem et al. 2022). Thus, there is more usable information in the tSZ field that we have not accessed in our datavectors given we limit our maps (and therefore our analysis) to a coarser resolution. This is done primarily in the interest of minimizing the computation cost, as generating a tSZ map of NSIDE = 8192 would take around 64 × longer than generating one of NSIDE = 1024. For the lensing analysis, the choice of NSIDE = 1024 is realistic, with some potential to extend to higher resolution if the systematics can be controlled robustly.  \nThe 1D marginal distributions in Figure 11 also shows a  \nµ  \nFigure 12. Fisher constraints of the baryonification parameters from different survey datasets. The datavector includes only the 2nd and 3rd moments, which contain the bulk of the signal-to-noise. The constraints from LSST Y10 and SO is between 3 to 5 times better than that from DES Y6 and SPT-3G. We have marginalized over all other parameters in the analysis. The constraints are also listed in Table 2.  \n<!-- image -->  \ndotted blue line, which are the results from a lensing-only analysis (using 2nd and 3rd moments) but with all baryonification parameters fixed. The constraints on cosmology 𝜎 ( Ω 𝑚 ) = 0 . 008 and 𝜎 ( 𝜎 8 ) = 0 . 013 - are comparable to those from using all 2nd and 3rd auto/cross-moments of the lensing and tSZ field after marginalizing over all baryonification parameters. Note that this compares lensing-only with a lensing plus tSZ analysis. In Appendix E, we consistently use a lensing plus tSZ analysis, with all baryonification parameters either fixed or varied in different stages, and find the constraints on cosmology degrades by ≈ 60% if we vary all baryonification parameters.  \nFinally, Figure 12 compares constraints from current surveys (DES Y6 and SPT-3G) with upcoming ones (LSST Y10 and SO). We find the parameter degeneracies in the upcoming surveys are broadly similar to those in the current surveys. The extended redshift range probed by LSST Y10 (compared to DES Y6) has not changed the degeneracy directions. The parameter constraints themselves improve by factors of 3 to 5 depending on the specific parameter. While these results are simplistic in their accounting of systematic uncertainties - we have not marginalized over intrinsic alignments for the lensing field, and we include no systematic model to be marginalized over for the tSZ field - they show the gain that could be realized from the final survey datasets of the next decade, relative  \nto constraints from ongoing surveys.  \nAll analyses above, summarised in Table 2, show that a numberofanalysis setups result in parameter uncertainties that are of the same order as the fiducial parameter values listed in Table 1: (i) the 𝜈 𝑋 parameters, which characterize the redshift evolution of 𝜃 ej and 𝑀 c, (ii) 𝛾 which is the gas slope at 𝑅 ∼ 𝑅 ej, and (iii) the mass-dependence parameters 𝜇 𝑋 . This may indicate that robust predictions of the lensing and tSZ moments can be obtained while fixing these parameters. However, a robust validation of this statement requires fitting this model (with and without varying the parameters mentioned here) to mock surveys derived from different hydrodynamical simulations. Appendix D also presents results that indicate the 𝜈 𝑋 parameters may not be necessary for jointly modelling the moments measurements across redshift. Bigwood et al. (2024), who use the baryonification model of Giri & Schneider (2021) to analyse data from DES and the Atacama Cosmology Telescope (ACT), also find qualitatively similar behavior to us for constraints on their mass-dependence parameter 𝜇 (see their Table B2).", '6 CONCLUSIONS': 'Cosmology is in a golden age of data-driven science, and this is set to only improve over the upcoming decade. For example, measurements of weak lensing will extend to larger fractions of the sky, and to much higher statistical precision (Spergel et al. 2015; The LSST Dark Energy Science Collaboration et al. 2018; Racca et al. 2016). Such measurements will access a broader range of scales than was previously available. However, the actual use of these measurements - namely those on smaller, more non-linear scales - for obtaining constraints is a non-trivial task due to uncertainties in modelling the physics of baryons on these scales. Addressing this issue is tantamount to better constraining Λ CDM cosmology ( e.g. , Secco et al. 2022b; Amon et al. 2022; Amon & Efstathiou 2022) and also to using lensing to extract information about other cosmological phenomenon, such as inflation (Anbajagane et al. 2024b; Goldstein et al. 2024).  \nWe introduce here a map-level baryonification technique that jointly generates baryon-corrected density fields (and thereby lensing fields) as well as tSZ fields, using N-body simulations as an input. This is done in a computationally efficient manner by working directly with full-sky maps from simulations rather than with three-dimensional particle snapshots. We perform various validations on this technique and characterize its performance on 2D simulation fields, as well as in mock analyses of wide-sky datasets. We perform all validations by using the moments of the fields, from 2nd to 4th order, as our summary statistics. We include all autoand cross-correlations at these orders. Our main results are as follows:  \n- · The baryonification method can be applied directly to 2D fields (rather than 3D snapshots) in a robust manner. After accounting for the pixel window function and the line-of-sight projection scale, the residuals between the 2D and 3D version is under 2% for most physical scales (Figure 3).\n- · Themodelpredictions are fairly insensitive to secondary halo properties, such as halo concentration and ellipticity. (Figure 5 and 7).\n- · The model is adequately flexible and can jointly fit all moments (up to 4th order) of the total matter density and gas pressure fields, as measured in I/l.pc/l.pc/u.pc/s.pc/t.pc/r.pc/i.pc/s.pcTNG. The residuals are within the uncertainties in all cases (Figure 8). The model can also jointly fit across multiple redshifts (Figure D1).\n- · The moments have significant differences in the halo mass their signal is most sensitive to; the total mass range spans nearly two orders of magnitude across the different moments, highlighting their complementary nature (Figure 9).\n- · A Fisher forecast of current/upcoming surveys shows the inclusion of higher-order information dramatically improves the constraining power on the baryonification model (Figure 10).\n- · The tSZ and lensing fields are more constraining for specific subsets of parameters, and the parameter degeneracy directions in the different fields are nearly orthogonal, leading to the combination of fields providing significantly better constraints (Figure 11). The lensing data is also more informative on redshift-dependent parameters due to the tomographic redshift information of the dataset.  \nThis work provides a first, initial confirmation that the baryonification method is viable for jointly modelling the lensing and tSZ fields directly on the sky, across multiple orders in the fields. We stress that this is not a final and conclusive validation of the map-level baryonification technique. As was pointed out before, a more thorough validation will require full-sky hydrodynamical simulations, such as M/i.pc/l.pc/l.pc/e.pc/n.pc/i.pc/u.pc/m.pc-TNG (Pakmor et al. 2023) and/or the F/l.pc/a.pc/m.pc/i.pc/n.pc/g.pc/o.pc project (Schaye et al. 2023), and a focus on other popular summary statistics used for analyses of weak lensing and tSZ fields.  \nOnce the baryonification method has been more precisely validated, a limiting step in applying it to data will be the forward model of the tSZ field and its foregrounds, done consistently across multiple orders . Omori (2022) have provided a precisely validated prescription for converting an N-body  \nsimulation into realistic tSZ fields (with all foregrounds consistently modelled, and then also verified at the two-point level), but this requires high fidelity simulation products which are often not available for the large simulation suites used for forward modelling the lensing field (Kacprzak et al. 2023; Anbajagane et al. 2024b; Jeffrey et al. 2024; Gatti et al. 2024a). Explorations of more approximate tSZ forward models will be needed in order to analyse the tSZ field with the same simulation-based forward modelling approach used for the lensing field.  \nThe baryonification method has rapidly become a important tool in characterizing the two-point correlation functions measured by lensing surveys. The method, however, is a general prescription on how to add baryonic fields to N-body simulations through an approximate, phenomological approach. In this work, we use this formalism to model the lensing and tSZ field. One could similarly use this prescription to model other fields that are observables of gas physics, such as the kinematic SZ effect, which has already been measured using different combinations of CMB and galaxy surveys ( e.g. , Mallaby-Kay et al. 2023; Hadzhiyska et al. 2024), and the X-ray luminosity, measured for a large fraction of the sky by the eROSITA survey (Bulbul et al. 2024). In fact, some of these observables have already been used to constrain baryonification models (Schneider et al. 2019, 2022; Grandis et al. 2024; Bigwood et al. 2024). Thus, we have only probed some aspects of the full potential provided by the baryonification formalism.  \nThe current state of widefield survey science is rapidly advancing towards surveys that cover O( 1 ) fractions of the sky, and have overlapping areas of O( 10 4 ) deg 2 . In this datarich age, the advent of multi-survey, multi-wavelength science provides powerful ways to extract robust physical inferences at ever-increasing precision. The baryonification method provides one way of performing joint-analyses across different surveys in a rapid, flexible way, and is complementary to other approaches such as dedicated hydrodynamical simulations. Our prescription for the map-level baryonification pipeline is documented and made publicly available. We hope this enables easy explorations/extensions of the baryonification pipeline that build up the synergies of the growing multisurvey landscape.', 'ACKNOWLEDGEMENTS': 'DAissupported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE 1746045. SP is supported by the Simons Collaboration on Learning the Universe. CC is supported by NSF grant AST-2306166.  \nWe thank Marco Gatti and Yuuki Omori for discussions during early stages of this work, Srini Raghunathan for guid-  \nance on the tSZ forward models used in Raghunathan et al. (2022a,b), and Kayla Kornoelje for details on the linear combination method used by SPT. We also thank the referee for their many helpful comments that greatly enhanced our discussion and results. The conception and pursuit of this work arose from many conversations between DA and Eric Baxter. The authors are grateful to Eric for his collaboration and friendship over the many years, and for his efforts in making our research communities vibrant avenues for discussion and discovery. He will be missed.  \nAll analysis in this work was enabled greatly by the following software: P/a.pc/n.pc/d.pc/a.pc/s.pc (McKinney 2011), N/u.pc/m.pcP/y.pc (Van der Walt et al. 2011), S/c.pc/i.pcP/y.pc (Virtanen et al. 2020), and M/a.pc/t.pc/p.pc/l.pc/o.pc/t.pc/l.pc/i.pc/b.pc (Hunter 2007). We have also used the Astrophysics Data Service (ADS) and arXiv preprint repository extensively during this project and the writing of the paper.', 'DATA AVAILABILITY': 'All simulation datasets used in this work are publicly available, and can be found at the links provided in the relevant subsections of Section 3.  \nThe baryonification pipeline is publicly available at https://github.com/DhayaaAnbajagane/ BaryonForge .', 'REFERENCES': "Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435  \nAbbott T. M. C., et al., 2022, Phys. Rev. D, 105, 023520  \nAde P., et al., 2019, J. Cosmology Astropart. Phys., 2019, 056  \nAmodeo S., et al., 2021, Phys. Rev. D, 103, 063514  \nAmon A., Efstathiou G., 2022, MNRAS, 516, 5355  \nAmon A., et al., 2022, Phys. Rev. D, 105, 023514  \nAnbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686  \nAnbajagane D., Evrard A. E., Farahi A., 2022a, MNRAS, 509, 3441  \nAnbajagane D., et al., 2022b, MNRAS, 510, 2980  \nAnbajagane D., et al., 2022c, MNRAS, 514, 1645  \nAnbajagane D., et al., 2023, MNRAS, 526, 5530  \nAnbajagane D., et al., 2024a, MNRAS, 527, 9378  \nAnbajagane D., Chang C., Lee H., Gatti M., 2024b, J. Cosmology Astropart. Phys., 2024, 062  \nAric'o G., Angulo R. E., 2024, arXiv e-prints, p. arXiv:2406.01672  \n- Aric'o G., Angulo R. E., Hern'andez-Monteagudo C., Contreras S., Zennaro M., 2021a, MNRAS, 503, 3596\n- Aric'o G., Angulo R. E., Contreras S., Ondaro-Mallea L., Pellejero-Iba˜nez M., Zennaro M., 2021b, MNRAS, 506, 4070  \nAric'o G., Angulo R. E., Zennaro M., Contreras S., Chen A.,  \nHern'andez-Monteagudo C., 2023, arXiv e-prints, p. arXiv:2303.05537  \nAsgari M., et al., 2021, A&A, 645, A104  \nAung H., Nagai D., Lau E. T., 2020, arXiv e-prints, p. arXiv:2012.00977  \n- Ayromlou M., Nelson D., Pillepich A., 2022, arXiv e-prints, p. arXiv:2211.07659  \nBattaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 75  \nPandey S., Salcido J., To C.-H., Hill J. C., Anbajagane D., Baxter E. J., McCarthy I. G., 2024, arXiv e-prints, p. arXiv:2401.18072  \nPeel A., Pettorino V., Giocoli C., Starck J.-L., Baldi M., 2018, A&A, 619, A38  \nPetri A., Liu J., Haiman Z., May M., Hui L., Kratochvil J. M., 2015, Phys. Rev. D, 91, 103511  \nPillepich A., et al., 2018a, MNRAS, 473, 4077  \nPillepich A., et al., 2018b, MNRAS, 475, 648  \nPlanck Collaboration et al., 2016, A&A, 594, A13  \nPlanck Collaboration et al., 2020, A&A, 641, A1  \nPratt C. T., Qu Z., Bregman J. N., 2021, ApJ, 920, 104  \nQuilis V., M J., Ib´a˜nez ., S´aez D., 1998, ApJ, 502, 518  \nRacca G. D., et al., 2016, in MacEwen H. A., Fazio G. G., Lystrup M., Batalha N., Siegler N., Tong E. C., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave. p. 99040O ( arXiv:1610.05508 ), doi:10.1117/12.2230762  \nRagagnin A., Dolag K., Moscardini L., Biviano A., D'Onofrio M., 2019, MNRAS, 486, 4001  \nRaghunathan S., et al., 2022a, ApJ, 926, 172  \nRaghunathan S., et al., 2022b, ApJ, 928, 16  \nRankine W. J. M., 1870, Philosophical Transactions of the Royal Society of London, 160, 277  \nSchaye J., et al., 2023, MNRAS, 526, 4978  \nSchneider A., Teyssier R., 2015, J. Cosmology Astropart. Phys., 2015, 049  \nSchneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmology Astropart. Phys., 2019, 020  \nSchneider A., Giri S. K., Amodeo S., Refregier A., 2022, MNRAS, 514, 3802  \nSecco L. F., et al., 2022a, Phys. Rev. D, 105, 023515  \nSecco L. F., et al., 2022b, Phys. Rev. D, 105, 023515  \nSevilla-Noarbe I., et al., 2021, ApJS, 254, 24  \nShao M., Anbajagane D., 2024, The Open Journal of Astrophysics, 7, 29  \nShao M. J., Anbajagane D., Chang C., 2023, MNRAS, 523, 3258  \nShaw L. D., Nagai D., Bhattacharya S., Lau E. T., 2010, ApJ, 725, 1452  \nSheth R. K., Tormen G., 1999, MNRAS, 308, 119  \nSobol I. M., 1967, USSR Computational Mathematics and Mathematical Physics, 7, 86  \nSorini D., Bose S., Pakmor R., Hernquist L., Springel V., Hadzhiyska B., Hern'andez-Aguayo C., Kannan R., 2024, arXiv e-prints, p. arXiv:2409.01758  \nSpergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757  \nSpringel V., et al., 2018, MNRAS, 475, 676  \nStiskalek R., Bartlett D. J., Desmond H., Anbajagane D., 2022, MNRAS, 514, 4026  \nSunyaev R. A., Zeldovich Y. B., 1972, Comments on Astrophysics and Space Physics, 4, 173  \nTeyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195  \nThe Dark Energy Survey Collaboration 2005, arXiv e-prints, pp astro-ph/0510346  \nThe LSST Dark Energy Science Collaboration et al., 2018, arXiv e-prints, p. arXiv:1809.01669  \nThornton R. J., et al., 2016, ApJS, 227, 21  \nTo C.-H., Pandey S., Krause E., Dalal N., Anbajagane D., Weinberg D. H., 2024, arXiv e-prints, p. arXiv:2402.00110  \nTroxel M. A., Ishak M., 2015, Phys. Rep., 558, 1  \nVan Waerbeke L., et al., 2013, MNRAS, 433, 3373  \nVan der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22  \nVillaescusa-Navarro F., et al., 2020, ApJS, 250, 2  \nVirtanen P., et al., 2020, Nature Methods, 17, 261  \nVogelsberger M., et al., 2014, Nature, 509, 177  \nVogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42  \nWechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002,  \nApJ, 568, 52 Weinberger R., et al., 2017, MNRAS, 465, 3291 Zentner A. R., Hearin A. P., van den Bosch F. C., 2014, MNRAS, 443, 3044 Zhang Z. J., Chang C., Larsen P., Secco L. F., Zuntz J., LSST Dark Energy Science Collaboration 2022, MNRAS, 514, 2181 Zhu Z., et al., 2021, in Galaxy Cluster Formation II (GCF 2021) - Virtual  \nWorkshop. p. 56 ( arXiv:2107.01663 ), doi:10.5281/zenodo.5013757  \nZurcher D., Fluri J., Sgier R., Kacprzak T., Refregier A., 2021, J. Cosmology Astropart. Phys., 2021, 028  \nZurcher D., et al., 2022, MNRAS, 511, 2075", 'A FORWARDMODELLING PIPELINE': "Akeyresult of our work is the simulation-based forecast of current and future surveys (Section 5). The forecast is done by forward modeling the observables from these surveys; namely, the lensing convergence field (Section A1) and the thermal Sunyaev-Zeldovich field (Section A2). This appendix details the exact procedures used to create our forward model for each field.  \nThroughout, we use maps with NSIDE = 1024. The angular scale of the pixels (3 . 4 ' ) correspond to physical scales of 0 . 4 < 𝑟 < 7 Mpc for the redshift range spanned by the simulation, 0 . 1 < 𝑧 < 3 . 5.", 'A1 Weak lensing': "The lensing convergence field, 𝜅 , is a line-of-sight integral of the density field  \n𝜅 ( ˆ n , 𝑧 𝑠 ) = 3 2 𝐻 2 0 Ω m 𝑐 2 ∫ 𝑧 𝑠 0 𝛿 ( ˆ n , 𝑧 𝑗 ) 𝜒 𝑗 ( 𝜒 𝑠 -𝜒 𝑗 ) 𝑎 ( 𝑧 𝑗 ) 𝜒 𝑠 𝑑𝑧 𝑗 𝑑𝜒 𝑑𝑧 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝑧 𝑗 , (A1)  \nwhere 𝑧 𝑠 is the redshift of the 'source' plane/galaxies being lensed, ˆ n is the pointing direction on the sky, 𝛿 is the overdensity field, 𝜒 is the comoving distance from an observer to a given redshift, 𝑎 is the scale factor, 𝐻 0 is the Hubble constant, Ω m is the matter energy density fraction at 𝑧 = 0, and 𝑐 is the speed of light. Here, we have used the shorthand 𝜒 ( 𝑧 𝑠 ) ≡ 𝜒 𝑠 and 𝜒 ( 𝑧 𝑗 ) ≡ 𝜒 𝑗 .  \nWe use the same forward-modelling pipeline from Anbajagane et al. (2023, 2024b), whose description we reproduce below for completeness. The lensing maps are generated from the density fields of the U/l.pc/a.pc/g.pc/a.pc/m.pc simulations, and then processed to include all the observational effects found in the data. These procedures have been utilized in many analyses/forecasts of weak lensing data (Fluri et al. 2019; Zurcher  \net al. 2021; Fluri et al. 2022; Gatti et al. 2022; Zurcher et al. 2022; Gatti et al. 2023; Anbajagane et al. 2023, 2024b; Gatti et al. 2024b,a; Jeffrey et al. 2024). As mentioned in Section 5, we focus on two surveys: (i) the Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005), Year 6 dataset, and the (ii) the Rubin Observatory Legacy Survey of Space and Time (LSST), year 10 dataset. We detail below our forward modeling procedure:  \nConstructing lensing convergence shells. The model starts with lightcone shells of the particle counts, which are two-dimensional, HEALP/i.pc/x.pc maps of the (projected) particle counts at different redshifts. We convert the particle count map to an overdensity map as  \n𝛿 𝑖 = 𝑁 𝑖 p /⟨ 𝑁 p ⟩ -1 , (A2)  \nwhere 𝑁 𝑖 p is the number of particles in pixel 𝑖 , and the average is computed over all pixels in the shell. The density shells can then be converted into the convergence 𝜅 using Equation A1, after converting the integral over redshift into a discrete sum over lightcone shells. Note that our baryonification step is performed separately on each individual density shell, using the halos in just that shell. The baryonified density shells are then used in the summation of Equation (A1).  \nSource galaxy redshift distributions. We perform a weighted average of the convergence shells to construct the convergence within the different tomographic bins of each survey. The weights used in this averaging are the source galaxy redshift distribution, 𝑛 ( 𝑧 ) , of the chosen bin and survey,  \n𝜅 𝐴 ( ˆ n ) = 𝑁 steps ∑︁ 𝑗 = 1 𝑛 𝐴 ( 𝑧 𝑗 ) 𝜅 ( ˆ n , 𝑧 𝑗 ) Δ 𝑧, (A3)  \nwhere 𝜅 𝐴 is the true convergence of a tomographic bin, 𝐴 . The 𝑛 ( 𝑧 ) is obtained from the following: for DES Year 6 and LSST Year 10 we use the same 𝑛 ( 𝑧 ) as that used in Zhang et al. (2022, see their Table 1). The LSST modeling in that work follows the baseline analysis choices of The LSST Dark Energy Science Collaboration et al. (2018, see their Appendix D2.1). The redshift distributions for DES Y6, and LSST Y1 and Y10 are parameterized as,  \n𝑑𝑁 𝑑𝑧 ∝ 𝑧 2 exp GLYPH<20> -GLYPH<18> 𝑧 𝑧 0 GLYPH<19> 𝛼 GLYPH<21> , (A4)  \nwith parameters given in Table A1. Once the 𝑛 ( 𝑧 ) of the full survey is defined, we split it into 4 (5) tomographic bins for DES Y6 (Y10) of equal number density. Each bin is then convolved with a Gaussian of width given by the photometric redshift uncertainty, also quoted in Table A1. The DES Y6 distribution is non-zero only between 0 . 2 < 𝑧 < 1 . 3, following Zhang et al. (2022, see their Table 1). The LSST distributions are cut at 𝑧 < 3 . 5. See Figure 1 of Anbajagane et al. (2024b) for the exact distributions we use in this work.  \nTableA1. Theredshift distribution and source galaxy number density assumed for the upcoming surveys. All numbers are taken from Zhang et al. (2022, see their Table 1).  \nConstructing lensing shear shells. Weak lensing surveys measure galaxy shapes, which primarily trace the shear field, 𝛾 , and not the convergence field, 𝜅 . However, the shear and convergence field can be transformed into each other using the Kaiser-Squires (KS) transform (Kaiser & Squires 1993), implemented in harmonic space as  \n𝛾 ℓ𝑚 𝐸 + 𝑖𝛾 ℓ𝑚 𝐵 = -√︄ ( ℓ + 2 )( ℓ -1 ) ℓ ( ℓ + 1 ) GLYPH<18> 𝜅 ℓ𝑚 𝐸 + 𝑖𝜅 ℓ𝑚 𝐵 GLYPH<19> , (A5)  \nwhere 𝑋 { 𝐸,𝐵 } are the E-mode and B-mode (or Q and U polarizations, in H/e.pc/a.pc/l.pcP/i.pc/x.pc notation) of the field.  \nShape noise. Upon generating the two shear fields, we add shape noise in real space. The forward-modelled field includes Gaussian shape noise with a standard deviation given as  \n𝜎 𝛾 = 𝜎 𝑒 √︁ 𝑛 gal 𝐴 pix , (A6)  \nwhere 𝑛 gal is the source galaxy number density, and 𝐴 pix is the pixel area for a given map resolution. All maps in this work use NSIDE = 1024, corresponding to a pixel resolution of 3 . 2 ' . The per-galaxy shape noise is taken to be 𝜎 𝑒 = 0 . 26.  \nSurvey mask & mass map construction. The noisy shear field - which is the sum of the true shear fields and the shape noise fields - is then masked according to the survey footprint. For DES Y6 we use the provided survey mask in the year 3 data release (Sevilla-Noarbe et al. 2021). For LSST Y10, we divide the full-sky into three equal-area cutouts of roughly 14,000 deg 2 each, which is the expected area coverage (The LSST Dark Energy Science Collaboration et al. 2018, see Appendix C1). The noisy shear maps are converted to convergence using the relation in Equation A5. We only use the resulting E-mode field, 𝜅 𝐸 , for our analyses. This follows the same procedures used in Chang et al. (2018); Jeffrey et al. (2021).  \nTo summarize, lensing convergence maps are constructed from the raw particle number count maps. The 𝑛 ( 𝑧 ) distributions for a given survey are used to obtain the convergence map in a given tomographic bin. This convergence map is converted to shear maps, the relevant shape noise is added, the relevant survey mask is applied, and then the noisy shear maps are converted back to a noisy convergence map. The set of procedures listed above is the standard approach for forward modeling the lensing field ( e.g. , Zurcher et al. 2021, 2022; Gatti et al. 2022; Anbajagane et al. 2023, 2024b). Thus, our final convergence maps will be an accurate representation of  \nthe survey data. We have intentionally ignored the modelling of intrinsic alignments (Troxel & Ishak 2015; Lamman et al. 2023) - which is a systematic effect considered in all the models cited above - as this work focuses on a pure test of the degeneracy directions between cosmology and the baryonification paraemeters. We will explore, in future works, the relationship between the baryonification model and other systematics.", 'A2 Thermal Sunayev-Zeldovich effect': "Therawmeasurements from a CMB survey are the temperature or power measured at different frequencies. At a given frequency, the observations contain contributions from a wide variety of astrophysical and cosmological sources, of which one is the tSZ field. This tSZ field can be filtered out of these maps through a 'Linear Combination' (LC) procedure, which is a weighted sum of frequency maps performed in harmonic space ( e.g. , Madhavacheril et al. 2020; Bleem et al. 2022). The exact method used to calculate the weights results in a variety of LC methods, such as the internal Linear combination (ILC), Needlet ILC etc. We use the LC algorithm of Bleem et al. (2022), which uses a theoretically estimated covariance matrix for the angular power spectra of the frequency maps. This estimate requires models for the other, 'foreground' fields that contribute power to the different frequency maps. For this, we follow the approach of Raghunathan et al. (2022a,b), who base their model on the measurements of George et al. (2015) from the South Pole Telescope (SPT) data.  \nAs mentioned above in Section 5, we consider two CMB surveys: SPT-3G (Benson et al. 2014) and the Simons Observatory (SO, Ade et al. 2019). For both surveys, we assume their analyses are supplemented by Planck data, following the procedures used in existing surveys (Madhavacheril et al. 2020; Bleem et al. 2022). As a result, we forward model the Planck maps in this work. We follow Bleem et al. (2022) in using data from only the high-frequency instrument (HFI) of Planck . We also additionally use the 545 GHz band (which is not used in Bleem et al. (2022)), but as we will show later in Figure A1, the weights for the map is such that it provides no information to the final data product. The characteristics of the surveys (noise properties, beam, frequency bands etc.) are given in Table A2. We now detail our forward modelling procedure for the tSZ field:  \nTrue tSZ field. Our simulated tSZ emission is obtained by pasting pressure profiles onto a H/e.pc/a.pc/l.pcP/i.pc/x.pc map. We use all halos in the lightcone to do so. The pasted profiles are already integrated over the line of sight. We therefore only need to change the units of the map through the simple rescaling,  \n𝑦 ( ˆ n ) = 𝜎 𝑇 𝑘 𝐵 𝑚 𝑒 𝑐 2 𝑃 𝑒 ( ˆ n ) (A7)  \nTable A2. The noise and beam properties of the different surveys in different frequency bands (in units of GHz). For Planck , we list the 100, 143, and 217 properties under the 90, 150, 220 columns for brevity; our actual analysis makes Planck frequency maps at 100, 143, and 217 GHz. The 𝛼 knee and ℓ knee of SO are constant across all bands. For SPT-3G, 𝑁 atm = 𝑁 white , and for Planck 𝑁 atm = 0. All noise estimates are in micro-kelvin-arcminute units. The beam full-width half-max, 𝜃 FWHM , is given in arcminutes. The values for SPT-3G and SO are taken from Raghunathan et al. (2022b, see their Table 1 and 2) while the values for Planck are from Madhavacheril et al. (2020, see their Table 1). See references therein for how these values are derived.  \nwhere 𝑃 𝑒 ( ˆ n ) is the projected electron pressure in direction ˆ n , 𝜎 𝑇 is the Thomson cross-section, 𝑘 𝐵 is the Boltzmann constant, and 𝑚 𝑒 𝑐 2 is the rest energy of an electron. As mentioned above, the raw dataproducts of a CMB survey are the temperature maps at different frequencies. We can convert the tSZ signal into these frequencies via,  \nΔ 𝑇 ( 𝜈, ˆ n ) = 𝑓 SZ ( 𝜈 ) 𝑇 cmb 𝑦 ( ˆ n ) (A8)  \nwhere the frequency response of the tSZ signal, 𝑓 SZ, is given as,  \n𝑓 SZ GLYPH<18> 𝑥 ≡ ℎ𝜈 𝑘 𝐵 𝑇 cmb GLYPH<19> = 𝑥 coth ( 𝑥 / 2 ) -4 . (A9)  \nHere, ℎ is the Planck constant and 𝑇 cmb = 2 . 7255 𝐾 is the temperature of the CMB.  \nForegrounds. As we discussed previously, the temperature map at a given frequency has contributions from multiple fields, normally called 'foregrounds'. We follow Raghunathan et al. (2022b) in modelling the CMB, radio, and infrared foregrounds. The kinematic SZ emission is ignored as it is subdominant, relative to these foregrounds, by more than an order of magnitude and is thus negligible for our work (George et al. 2015, see their Figure 3). We generate Gaussian realizations of all these foregrounds, by modelling their harmonic power spectra and using the synfast routine in H/e.pc/a.pc/l.pcP/y.pc to make full-sky realizations. The lensed 𝐶 CMB ℓ is modelled using C/a.pc/m.pc/b.pc run at the best-fit cosmology from Planck (Planck Collaboration et al. 2016).  \nIn differential temperature units ( Δ 𝑇 ), the CMB template has no frequency dependence. The radio and infrared foregrounds follow George et al. (2015), and can be written as  \n𝐶 Radio ℓ,𝜈 1 ,𝜈 2 = 𝐴 Radio 𝐺 15 𝐹 Radio 𝜈 1 ,𝜈 2 GLYPH<18> ℓ 3000 GLYPH<19> 2 GLYPH<18> 2 𝜋 ℓ ( ℓ + 1 ) GLYPH<19> (A10)  \n𝐶 CIB , poiss . ℓ,𝜈 1 ,𝜈 2 = 𝐴 CIB , poiss . 𝐺 15 𝐹 CIB , poiss . 𝜈 1 ,𝜈 2 GLYPH<18> ℓ 3000 GLYPH<19> 2 GLYPH<18> 2 𝜋 ℓ ( ℓ + 1 ) GLYPH<19> (A11)  \n𝐶 CIB , clus . ℓ,𝜈 1 ,𝜈 2 = 𝐴 CIB , clus . 𝐺 15 𝐹 CIB , clus . 𝜈 1 ,𝜈 2 GLYPH<18> ℓ 3000 GLYPH<19> 0 . 8 GLYPH<18> 2 𝜋 ℓ ( ℓ + 1 ) GLYPH<19> (A12)  \nThe radio background is modelled as having a poisson spatial distribution (and is therefore essentially constant in ℓ ), while the CIB is modelled as a combination of a poisson-distributed component and a clustered component. The frequency response function, 𝐹 𝜈 1 ,𝜈 2 is defined as,  \n𝜖 𝜈 1 ,𝜈 2 = GLYPH<18> 𝐵 ' ( 𝜈 0 ) 𝐵 ' ( 𝜈 0 ) 𝐵 ' ( 𝜈 1 ) 𝐵 ' ( 𝜈 2 ) GLYPH<19> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> T = Tcmb (A13)  \n𝐹 Radio 𝜈 1 ,𝜈 2 = 𝜖 𝜈 1 ,𝜈 2 GLYPH<18> 𝜈 1 𝜈 2 𝜈 2 0 GLYPH<19> -𝛼 (A14)  \n𝐹 CIB , x 𝜈 1 ,𝜈 2 = 𝜖 𝜈 1 ,𝜈 2 GLYPH<18> 𝜈 1 𝜈 2 𝜈 2 0 GLYPH<19> 𝛽 𝑥 GLYPH<18> 𝐵 ( 𝜈 1 ) 𝐵 ( 𝜈 2 ) 𝐵 ( 𝜈 0 ) 2 GLYPH<19> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> T = Tcib (A15)  \nwhere 𝐵 is the blackbody function, 𝐵 ' ≡ 𝑑𝐵 / 𝑑𝑇 its derivative with temperature, and 𝑥 ∈ { poiss . , clus . } . The function 𝜖 𝜈 1 ,𝜈 2 converts spectral energies to differential temperature ( Δ 𝑇 ) units. The radio component is modelled with a power-law frequency dependence, while the CIB is a modified blackbody where the modification is to the low-frequency, Rayleigh-Jeans regime of the distribution. For the CIB model in Equation (A15), we use 𝛽 poiss . = 1 . 505, 𝛽 clus . = 2 . 510, 𝑇 cib = 20 𝐾 (George et al. 2015). For the radio foreground, we use 𝛼 = -0 . 6, following Raghunathan et al. (2022b).  \nNote that our goal in including foregrounds is only to obtain a reasonably accurate noise term for the tSZ maps. As a result, we follow Raghunathan et al. (2022a,b) in modelling these foregrounds as uncorrelated Gaussian fields, whereas in actuality these fields will be correlated. Such correlations were found to be negligible for analyses of clusters (Raghunathan et al. 2022a, see their Section 3.6), which dominate the information in the tSZ field. However, a more realistic forward model of the tSZ field will need to consistently connect the foregrounds to the density field and halo catalog in a given simulation realization; see Omori (2022), for an example of a robust, self-consistent model.  \nInstrumental and Atmospheric noise. The primary noise contributor to the tSZ maps are the thermal noise of the detector and atmospheric noise (for ground-based surveys). This is modelled through a two-term noise component,  \n𝑁 ℓ = 𝑁 white + 𝑁 atm GLYPH<18> ℓ ℓ knee GLYPH<19> -𝛼 knee . (A16)  \nThe first term is the frequency-independent, 'white noise' component arising from detector noise. The second term is a 1 / 𝑓 -noise that increases towards large-scales. The values we use for the different experiments are given in Table A2. For Planck , we set 𝑁 atm = 0. /one.sup/nine.sup Equation (A16) is the beamconvolved noise estimate, and to get the true noise model we deconvolve Equation (A16) as 𝑁 ℓ → 𝑁 ℓ / 𝐵 2 ℓ , where the beam 𝐵 ℓ is defined below.  \nInstrument Beam. The point-spread function (or beam) varies by experiment and detector. The values for the experiments we consider are given in Table A2. We use Gaussian beams,  \n𝐵 ℓ = exp GLYPH<18> 𝜎 2 FWHM 2 ℓ ( ℓ + 1 ) GLYPH<19> (A17)  \nwith 𝜎 FWHM = 𝜃 FWHM / √︁ 8 ln ( 2 ) . Here 𝜃 FWHM is the beam full-width half-max in radians.  \nLinear Combination (LC): Finally, the (beamdeconvolved) frequency maps are summed together through an LC, similar to that of Bleem et al. (2022). In this formalism, the final map is, /two.sup/zero.sup  \n𝑦 LC ℓ,𝑚 = ∑︁ 𝑖 𝑤 ℓ ( 𝜈 𝑖 ) 𝑇 ℓ,𝑚 ( 𝜈 𝑖 ) (A18)  \nwhere 𝑇 ℓ ( 𝜈 𝑖 ) are (harmonic space) temperature maps at different frequencies and 𝑤 ℓ ( 𝜈 𝑖 ) are the LC weights. These weights are the optimal (minimum-variance, or MV) weights that maximimize the tSZ signal, defined as  \n𝑤 ℓ ( 𝜈 𝑖 ) = 𝑓 𝜈 𝑖 SZ 𝐶 -1 ℓ 𝒇 𝑇 SZ 𝐶 -1 ℓ 𝒇 SZ . (A19)  \nwhere f SZ = { 𝑓 𝜈 0 , 𝑓 𝜈 1 , . . . } specifies that the tSZ is our target signal, and 𝐶 -1 ℓ is the inverse of the covariance matrix between the different frequency maps at a given ℓ . Following Bleem et al. (2022); Raghunathan et al. (2022a) we model the covariance matrix theoretically as,  \n𝐶 𝜈 𝑖 ,𝜈 𝑗 ℓ = 𝑁 ℓ,𝜈 𝑖 ,𝜈 𝑗 𝛿 𝑖 𝑗 + ∑︁ 𝑎 ∈ 𝑋 𝐶 𝑎 ℓ,𝜈 𝑖 ,𝜈 𝑗 (A20)  \nwhere 𝑁 ℓ is already beam-deconvolved, and 𝛿 𝑖 𝑗 ensures that instrument and atmospheric noise only affects the autocorrelation of frequency maps. /two.sup/one.sup The frequency and ℓ -dependence of the foregrounds, 𝑋 ∈ { CMB , Radio , CIB } are defined above.  \n/one.sup/nine.sup Note that Planck still exhibits some 1 / 𝑓 -noise (particularly on ℓ < 100) arising from other, non-atmospheric components. For the purposes of this work, we ignore this term.  \n/two.sup/zero.sup Under the assumption of uncorrelated Gaussian foregrounds and instrumental/atmospheric noise, the residual noise in the tSZ field is given by the spectra 𝑁 𝑦𝑦 ℓ = 𝒇 𝑇 SZ 𝐶 -1 ℓ 𝒇 SZ , which is the denominator in Equation (A19). One can obtain a realistic, noisy tSZ field by simply generating a Gaussian field with 𝐶 ℓ = 𝑁 𝑦𝑦 ℓ and adding it to the true tSZ field. However, our forward model still adopts the approach of decomposing the signal into different frequency bins and performing the LC, as this approach is more generalized.  \n/two.sup/one.sup Nominally, atmospheric noise should be correlated across frequency bands, however we follow Raghunathan et al. (2022a,b) in ignoring this.  \nFigure A1. The harmonic-space weights used in the linear combination procedure - see Equation (A19) - to combine the different frequency maps into a minimum-variance tSZ map, shown for both SPT-3G (top) and Simons Observatory (bottom). The weights for the five Planck maps are shown as dashed lines. A black dotted line shows 𝑤 = 0 for reference.  \n<!-- image -->  \nThe emission from these foregrounds is completely correlated across bands. Once this final map, 𝑦 LC ℓ,𝑚 , is made, we convolve it with a Gaussian beam of 1 . 7 ' , which is the final resolution of the map. Note that in practice, we make maps at NSIDE = 1024, which already has a larger pixel scale (3 . 4 ' ) than the beam.  \nIn summary, we generate true tSZ maps from our baryonification pipeline and decompose the signal into different frequency channels. We then add multiple other foregrounds to these frequency channels, including the CMB, Radio emission, and CIB. We also add the thermal noise from detectors and then the noise due to the atmospheric. The maps are linearly combined using weights that result in the minimum variance of the final tSZ map. The weights we use in this work are shown in Figure A1. We have verified our weights estimation procedure reproduces the results of Raghunathan et al. (2022b, their Figures 3 and 4), and we also find our dominant frequency band at a given ℓ broadly follows that of Bleem et al. (2022); /two.sup/two.sup for example the 100 GHz Planck data dominates at low ℓ while the 90 GHz SPT data dominates at high ℓ .  \nOmori (2022) provide a detailed, and precisely validated, prescription for consistently producing a realistic tSZ field  \n/two.sup/two.sup We do not expect precise agreement given our approach is still a significantly simpler version of the analysis done in Bleem et al. (2022), and given that analysis is for a previous release of the SPT data.  \nFigure B1. The impact of the chosen projection scale on the projected displacement function for halos of two different masses; see Equation (25) and Section (2.2). Increasing the limits of the lineof-sight projection integral washes out the variations in the density profile arising from baryonic imprints and lowers the amplitude of the displacements.  \n<!-- image -->  \n- including all relevant contaminant foregrounds, lensing effects, cross-correlations etc. - through the use of highresolution simulation products and data-driven models. The U/l.pc/a.pc/g.pc/a.pc/m.pc suite does not have the required resolution or fidelity of simulation products (for example, we have no merger trees or particle shells) to employ their approach. Instead, our forward model focuses solely on capturing the noise level of the tSZmaps, which is adequate for understanding the different degeneracy directions (for the different moments measurements) of the different baryonification and cosmology parameters.", 'B IMPACT OF PROJECTION SCALE': 'We extend on the discussions of methodology validation in Section 4 by checking the sensitivity of our predictions to choices in the projection operation. As noted in Section 2.2, we set the limits of the projection integral according to the simulation box (or more generally, according to the thickness of the volume that we are applying baryonification on). This is in contrast to some previous works that set the the projection scale based on the projected separation 𝑟 p (see Equation (25) for definitions of the quantity). Figure B1 presents the impact of the projection choices on the estimated projected displacement function. In the limit of an infinite projection length, the impact of baryons will be completely washed out and there is no baryonification required. Figure B1 shows that increasing the projection length does indeed reduce the predicted dis-  \nFigure C1. The dependence of the displacement function on the including/excluding the DM adiabatic relaxation. The effect changes the displacement by only 10%, and thus, this displacement is fairly insensitive to the exact modelling choices of the DM relaxation and is instead highly sensitive to the gas distribution.  \n<!-- image -->  \ncements from the model. We use a minimum projection scale of 60 Mpc, which is the width of density shells used in full-sky mocks of the lensing field (Kacprzak et al. 2023; Anbajagane et al. 2024b), and a maximum of 1000 Mpc, which is twice the projection scale we use for the Q/u.pc/i.pc/j.pc/o.pc/t.pc/e.pc.  \nThe displacements on the largest scales are affected most by varying the projection scale ( e.g. , see bottom panel of Figure B1) as at these radii the one-halo density profile becomes comparable/subdominant to the large-scale two-halo term. On smaller scales the impact is less prominent. Though, for adequately small halos ( e.g. , , 𝑀 = 10 12 M ⊙ , see top panel of Figure B1), where the two-halo term can dominate over the one-halo term even at smaller radii, the suppression in amplitude is still seen as we increase the projection scale. This result highlights the importance in choosing a consistent, physically meaningful projection scale.', 'C ALTERNATIVE CHOICES IN HALO PROFILE MODEL': 'The halo model described in Section 2 is not a unique parameterization of the different matter components in the halo; there are several modifications to the model that either extend it to enhance accuracy or simplify it to enhance computational cost. We test a few of these modifications below. In all cases, we show the differences at the profile level for brevity and do not present the summary statistics-level consistency checks.', 'C1 Adiabatic contraction of dark matter': "The collisionless matter component of the baryonification model includes galaxies as well as DM. The distribution of  \nFigure C2. The profile of total gas pressure compared with only the thermal pressure, assuming some model for the non-thermal component. We show predictions from Green et al. (2020), and for the fiducial baryonification parameters (Table 1) following the model of Shaw et al. (2010) for the non-thermal component as described in Section 2.3.  \n<!-- image -->  \n200c  \nthe latter does not follow the same NFW profile as it did in the DMO case; this has been highlighted by many works on a multitude of simulations ( e.g. , Gnedin et al. 2004; Duffy et al. 2008; Ragagnin et al. 2019; Beltz-Mohrmann & Berlind 2021; Forouhar Moreno et al. 2021; Anbajagane et al. 2022a; Shao et al. 2023; Shao & Anbajagane 2024; Sorini et al. 2024). The gravitational response of the DM to the presence of baryons has been studied by many works, and approximated through various (simple) relations. Our baryonification model adopts one such relation from Gnedin et al. (2004); Abadi et al. (2010); Teyssier et al. (2011). In Figure C1 we show the change in the displacement function if we ignore the adiabatic relaxation of the DM. This provides an estimate of the model's sensitivity to inaccuracies in the modelling of this relaxation process. The effect is at most 15% of the maximum displacement value. Thus, the displacement function is dominated by the form of the gas distribution rather than by the redistribution of the DM due to adiabatic relaxation. This is consistent with findings in hydrodynamical simulations ( e.g. , Springel et al. 2018, see their Figure 9).", 'C2 Sensitivity to non-thermal pressure models': 'The tSZ modelling described in Section 2.3 starts by assuming halos are in hydrostatic equilibrium, which means the gas distribution is supported against gravitational collapse solely due to its thermal pressure. However, simulations have shown that order ∼ 10% of the pressure support comes from non-thermal pressure, primarily constituted by turbulent motions of the gas, i.e. a velocity dispersion (Nelson et al. 2014; Green et al. 2020). We have corrected for this in our predictions by explicitly modelling the contribution of nonthermal pressure and varying/marginalizing the parameters of  \nFigure C3. The dependence of the gas density profile (top) and displacement function (middle, bottom) when including features from a large-scale cosmological shock; see Equation (C1). We vary the shock location, 𝑅 sh , and the shock width Δ ln 𝑅 sh . In the top panel, the profiles varying the shock width are offset vertically for visual purposes. The bottom two panels show the shock feature changes the displacement quite significantly. However, these effects can be captured through the existing, flexible parameterization of the gas profile. See Section C3 for more details. The results from profiles with no shock features is shown in black.  \n<!-- image -->  \nthis model. Our parameterization follows Shaw et al. (2010), which is the same choice as P24.  \nFigure C2 compares this model to other alternatives, particularly the predictions of Green et al. (2020). The latter model is an analytic prediction (see their Figure 1 for a summary of the model) that depends only on the peak height of the halo, and is validated using a set of hydrodynamic simulations from Nelson et al. (2014). For comparison, we plot thermal gas pressure profiles, using either Green et al. (2020) or our prescription in Section 2.3 to estimate the non-thermal pressure. In the latter approach, predictions are made for fiducial parameter values listed in Table 1. The phenomenological model, based on Shaw et al. (2010), qualitatively matches the form predicted by Green et al. (2020), which in turn is a advanced version of previous non-thermal models, such as that of Nelson et al. (2014).', 'C3 Sensitivity to shock features': "The gas distribution around massive halos is expected to undergo shock heating due to the generic process of gravitational infall; cold matter from nearby large-scale structure (eg. filaments) infalls into the massive halo. Cold gas has a  \nrelatively low sound speed, but the infall velocity can exceed the sound speed by factors of O( 100 ) for gas around massive halos. This causes the gas to shock heat and thereby imprint features into the gas distribution on large scales. Such features have been predicted by simulations of gas ( e.g. , Quilis et al. 1998; Miniati et al. 2000; Molnar et al. 2009; Baxter et al. 2021; Aung et al. 2020), and there are indications of these features from measurements at different wavelengths ( e.g. , Hurier et al. 2019; Pratt et al. 2021; Zhu et al. 2021; Anbajagane et al. 2022c, 2024a; Hou et al. 2023).  \nThe change in density due to a cosmological shock can be found following the Rankine-Hugonoit jump conditions (Rankine 1870; Hugoniot 1887). For these shocks, which have a very high mach number 𝑀 ≡ 𝑣 / 𝑐 𝑠 ≫ 1, where 𝑣 is the infall velocity and 𝑐 𝑠 is the sound speed of the gas - the density profile post-shock is transformed as 𝜌 → 𝜌 / 4. This estimate corresponds to a gas adiabatic index of 𝛾 = 5 / 3, which is valid under the assumption that the gas is monoatomic. We then introduce this shock feature through a modified version of the gas profile from Section 2.1,  \n𝜌 gas , sh / 𝜌 gas = GLYPH<18> 1 -0 . 25 1 + exp GLYPH<2> 1 Δ ln 𝑅 shock ln ( 𝑟 / 𝑅 shock ) GLYPH<3> GLYPH<19> + 0 . 25 . (C1)  \nThe ratio is a simple sigmoid function that asymptotes to 1 / 4 for 𝑟 ≫ 𝑅 shock and 1 for 𝑟 ≪ 𝑅 shock. Note that in Equation (10), which defines the gas profile, the normalization of the profile is set by integrating the profile to 𝑟 → ∞ . Given the shock feature changes the shape of the profile, it would nominally affect this normalization. However, we keep the normalization fixed and only change the large-scale behavior of the profile. This is because the shock feature only affects the physics of the large scales and is largely decoupled from the gas physics internal to the halo. Thus, it is unrealistic to couple large-scale features like shocks to changes in the density throughout the halo (the normalization is a scaling factor for the entire profile and thus changes the density at all scales).  \nFigure C3 shows the profiles with shocks, and the displacement functions corresponding to such profiles. There is a significant effect both at the profile level and thereby at the displacement level. As discussed before, the feature is prominent at large-scales and so the displacement is affected only on these scales. The small-scale displacement is left largely unaltered. However, the baryonification model is flexible and therefore the large-scale changes to the gas profile can be captured through the gas profile slope parameters, 𝛾 and/or 𝛿 . The impact of such shocks on the pressure profile can also be captured by the non-thermal pressure model, as that model fundamentally captures a 'lack of thermal pressure' in the halo, which is the signature that will be generated by the shocks (Anbajagane et al. 2022c, see their introduction for a  \ndescription of this process). We leave a dedicated study of shocks in the baryonification model to a future study, and note that Figure C3 highlights the potential impact it can have on the total density distribution within halos.", 'D VALIDATION AT HIGHER REDSHIFTS': 'WhileSection 4.2 shows that the baryonification model can fit the simulation measurements at 𝑧 = 0, we must also verify that the model can jointly fit measurements across redshifts. Figure D1 shows such a joint fit across all moments, up to 4th order, from 𝑧 = 0 and 𝑧 = 0 . 5.  \nThe fit is done using three models: (i) the fiducial model shown in Figure 8, with 𝜈 𝑋 = 0 (no redshift dependence in the gas profile parameters) and fit to the 𝑧 = 0 data alone, then (ii) the same fiducial model now jointly fit to 𝑧 = 0and 𝑧 = 0 . 5, and finally (iii) the same model but with those redshift-dependence parameters varied, and fit to both 𝑧 = 0 and 𝑧 = 0 . 5 data. The first model is still a decent match to the simulations, though a few predictions at 𝑧 = 0 . 5 are discrepant with the simulation beyond the 3 𝜎 bound. The second model, which is a joint fit, is a better match across both redshifts. The predictions are almost always within 1 -2 𝜎 of the datapoints. Finally, the third model, which now varies two redshift-dependence parameters, 𝜈 𝜃 ej and 𝜈 𝑀 c , provides a similar/better match depending on the subset of measurements we focus on. There are some potential discrepancies in the higher order moments of the tSZ field, but the discrepancies are still within the 2 𝜎 uncertainties. The fractional uncertainties of the higher redshift measurements -particularly those involving the tSZ field, or higher-order combinations of any fields - can be larger due to the weaker signal at higher redshift. The statistical limitations of small simulation volumes, which we had discussed earlier in Section 4.2, are only exacerbated at higher redshift as there are now even fewer massive halos generating the signals.  \nIt is interesting to note that an adequate fit is achieved even when fixing the redshift-dependence parameters, 𝜈 𝑋 . This complements the results of Section 5, and Table 2 in particular, which showed that all versions of our analysis were poorly constraining these parameters. We have thus seen two different analyses show that the 𝜈 𝑋 parameters may not be a necessary feature of the model when applying baryonification to survey data. P24 found these parameters necessary for fitting the profiles of halos across multiple redshifts. However, it may be that for summaries of the whole lensing/tSZ field measured in surveys, the mass-dependence is adequate enough without the need for any additional redshift-dependence as well. However, one limitation of our analysis here is it spans a relatively narrow range in redshift, 0 < 𝑧 < 0 . 5 (though our Fisher forecast in Section 5 uses 0 < 𝑧 < 3 . 5) and fitting simulation measurements across a broader range could once again necessitate  \nTable E1. Constraints on two cosmology parameters, and the two key baryonification parameters (for lensing) as a function of the number of baryonifciation parameters varied in the analysis (also denoted with the brackets). We forecast for a DES Y6 x SPT 3G analysis, using the 2nd and 3rd moments of the fields.  \nthe redshift-dependence parameters in this particular analyses. Further analyses, using other simulations, is needed to better confirm the necessity of these parameters for analyses of survey data.  \nIn conclusion, we find the model specified in Section 2 adequately captures the moments of the lensing and tSZ fields, from 2nd to 4th order, and across multiple redshifts. As mentioned before, more precise validation will be enabled through the use of larger simulations.', 'E IMPACT OF PARAMETER SET OF FISHER INFORMATION': "Anumber of previous works have constrained the baryonification parameters using data ( e.g. , Chen et al. 2023; Aric'o et al. 2023; Grandis et al. 2023; Bigwood et al. 2024), and they span a variety of choices in the parameters being varied. For example Chen et al. (2023) vary a single parameter, 𝑀 𝑐 , whereas Aric'o et al. (2023) varies six parameters. Additionally, Giri & Schneider (2021) explored the minimal set of parameters needed for baryonification to still predict the matter power spectrum from different simulations.  \nHere, we show our Fisher constraints on cosmology using different sets of baryonification parameters. Ordered by increasing complexity and number of parameters, we consider: varying no baryonification parameters, varying only 𝑀 𝑐 , 𝜃 ej (base), varying the two mass-dependence parameters 𝜇 𝑋 , the two redshift dependence parameters 𝜈 𝑋 , and varying the four remaining parameters. See Table 1 for a complete list of all parameters we vary. Figure E1 shows the results of this lensing plus tSZ analysis. As expected, varying all parameters degrades cosmology constraints by roughly around 60%. Table E1 also lists the constraints for cosmology, as well as 𝑀 𝑐 , 𝜃 ej, alongside the number of baryonification parameters varied per analysis.", 'F COMPUTATIONAL DETAILS': "There are a number of salient computational choices and details in the baryonification model of Section 2 that we list  \nOnly z = 0 fit  \nz = 0, 0.5 fit  \nz = 0, 0.5 fit (z-dep. param.)  \nFigure D1. Similar to Figure 8, but now showing a variety of a different fitting cases and their predictions for 𝑧 = 0 . 5 measurements (top three rows) and for the 𝑧 = 0 measurements previously shown in Figure 8 (bottom). The different curves show fits made to different combinations of redshifts. The purple line ('fixed') is from a joint fit of 𝑧 = { 0 , 0 . 5 } with all redshift dependence of the baryonification model set to 𝜈 𝑋 = 0; see Equation (16). The yellow line ('varied') is the same joint fit, but after allowing these values to be non-zero. Both the fixed and varied predictions are a decent fit to the 𝑧 = 0 . 5 measurements, with the varied one providing the best fit of the two. The fit from using only the 𝑧 = 0 measurements (blue) is shown for comparison and is the worst performing. In the bottom three rows, the blue line overlaps exactly with the yellow and so we show it as a thicker/translucent line to aid visibility.  \n<!-- image -->", 'here for completeness.': "Projection: It is computationally efficient to perform the projection integral in Fourier space, which is the choice made in the CCL code base. However, this increases the likelihood of aliasing, or 'ringing' effects in the resulting profile, particularly for profiles computed in the extreme regions of parameter space. Such profiles can deviate very significantly from simple  \npower-laws, and in these cases the F/f.pc/t.pcL/o.pc/g.pc algorithm can fail. To make our pipeline more robust, we perform all projection integrals in real-space,  \n𝜌 𝑝 ( 𝑟 𝑝 ) = ∫ 𝐿 / 2 0 2 𝑑𝑙𝜌 GLYPH<18> √︃ 𝑙 2 + 𝑟 2 𝑝 GLYPH<19> . (F1)  \nWe find that using the real-space calculation, instead of the Fourier-space one, does not slow the pipeline a noticeable  \nFigure E1. Constraints on two cosmology parameters as a function of the number of baryonifciation parameters varied in the analysis. We forecast for a DES Y6 x SPT 3G analysis, using the 2nd and 3rd moments of the fields. The constraints, and the number of baryonification parameters varied in each, are listed in Table E1. The inclusion of all 11 baryonification parameters leads to a ≈ 60% degradation in cosmology constraints.  \n<!-- image -->  \namount. Note that in practice we only compute the profiles once, at the start of the baryonification pipeline, to create a table. This lookup table is then used for generating profiles of all the halos in the simulations (see the 'Tabulation' discussion below).  \nProfile interpolation: We use tabulation and interpolation for all integrals/derivatives of the profiles appearing in this work; for example, the volume integral to get the normalization in Equation (2) and (12). For almost all interpolation steps (the one exception is described below), we use the PchipInterpolator (Fritsch & Butland 1984; Moler 2004) in S/c.pc/i.pc/p.pc/y.pc as it preserves the monotonicity of the function under the presence of outliers. The CubicSpline interpolator, which is a popular choice and also implemented in S/c.pc/i.pc/p.pc/y.pc, can break monotonicity if the input function has any rapid (but still monotonic) variation and this will then lead to ringing in the interpolated result. While our profiles are all well-behaved and do not contain sharp changes over radius, their pixel-convolved versions may present such rapid deviations if the halo radius is far smaller than the pixel-scale. The PchipInterpolator enables the pipeline to adequately handle such extreme edgecases as well. We use this interpolator to accurately compute integrals of functions, for example when converting density profiles to enclosed mass profiles.  \nWe also require interpolation for computing the derivatives - primarily the operation in Equation (20) to compute the collisionless matter profile - and here we do use the  \nCubicSpline . This is because this spline guarantees continuous first and second derivatives in the interpolated function (whereas PchipInterpolator only guarantees continuous first derivatives). In Equation (20), we take the derivative of a enclosed mass function in order to obtain a density profile, and if the second derivative of the former function is discontinuous that will cause (unphysical) jumps in the resulting density profile. The CubicSpline prevents this by enforcing smooth second derivatives in the interpolator of the enclosed mass.  \nTabulation: For all of our work here, we use a tabulated profile to speed up our calculations. These profiles (or the displacement function) are computed over some specified range in log 10 𝑟 , log 10 𝑀 200c and log 10 ( 1 + 𝑧 ) . When computing the dependence of the profiles on external parameters, such as concentration 𝑐 200c or the projection distance 𝐿 proj (see Figures 4 or B1), we extend the tabulation to include this parameter. These tables are defined on regular grids, and we use simple linear interpolation of the table through the RegularGridInterpolator routine in S/c.pc/i.pc/p.pc/y.pc. This approach greatly increases the computational speed of the model predictions, and leads to negligible drops in accuracy given all quantities (profiles, displacement functions etc.) are smoothly varying across radius, mass, and redshift. Grids of O( 10 5 ) points (1000 points in log 10 𝑟 , 10 in log 10 𝑀 200c, and 10 in log 10 ( 1 + 𝑧 ) ) take less than five minutes to build, on a singlethreaded process of an Intel Broadwell CPU, and result in interpolators that are accurate to within < 0 . 1%.  \nComputing runtimes: The baryonification process is 'embarassingly parallel' given the contributions from halos are accumulated through a loop, and each iteration in the loop can be evaluated completely disjoint from the rest. Our pipeline is parallelized across all cores in a node using the J/o.pc/b.pc/l.pc/i.pc/b.pc routines, i.e. we have implemented shared-memory parallelization. For the standard full-sky simulation used in this work - 100 density shells of NSIDE = 1024 and halo catalogs with 𝑀 > 10 14 M ⊙ ; see Section 3.3 - baryonification of the density field is done in under 2 minutes on an Intel Broadwell chip with 40 cores. The profile pasting step, which uses approximately two million halos across 0 < 𝑧 < 3 . 5, generates the tSZ field in under 2 minutes as well. The runtime scales linearly with the number of halos in the simulation volume as well as the number of shells and number of pixels. If we used halo catalogs of 𝑀 > 10 13 M ⊙ , the run time would increase by a factor of two to three; a halo catalog with 𝑀 > 10 13 M ⊙ has 10 × more halos than one with 𝑀 > 10 14 M ⊙ but these additional, less-massive halos will span a smaller area on the sky which in turn will reduce the computational load.  \nThis paper was built using the Open Journal of Astrophysics L A T E X template. The OJA is a journal which provides fast and  \neasy peer review for new papers in the astro-ph section of the  \narXiv, making the reviewing process simpler for authors and referees alike. Learn more at http://astro.theoj.org ."}