Map-level baryonification: Efficient modelling of higher-order correlations in the weak lensing and thermal Sunyaev-Zeldovich fields

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- $N$th-order 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 simulation-based 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 map-level baryonification is publicly available at https://github.com/DhayaaAnbajagane/Baryonification.