Impact of baryons in cosmic shear analyses with tomographic aperture mass statistics

Non-Gaussian cosmic shear statistics based on weak-lensing aperture mass ($M_{\rm ap}$) maps can outperform the classical shear two-point correlation function ($\gamma$-2PCF) in terms of cosmological constraining power. Reaching the full potential of these new estimators however requires an accurate modeling of the physics of baryons as the extra non-Gaussian information mostly reside at small scales. We present one such modeling based on the Magneticum hydrodynamical simulation for the KiDS-450, DES-Y1, and a Euclid-like surveys. We compute the bias due to baryons on the lensing PDF and the distribution of peaks and voids in $M_{\rm ap}$ maps and propagate it to the cosmological forecasts on the structure growth parameter $S_8$, the matter density parameter $\Omega_{\rm m}$ and the dark energy equation of state $w_0$, using the SLICS and cosmo-SLICS sets of dark matter only simulations. We report a negative bias of a few percents on $S_8$ and $\Omega_{\rm m}$ and also measure a positive bias of the same level on $w_0$ when including a tomographic decomposition. These biases are increased to the order of $\sim 5$% when combining $M_{\rm ap}$ statistics with the $\gamma$-2PCF as these estimators show similar dependency on the AGN feedback. We verify that these biases constitute less than a $1\sigma$ shift on the probed cosmological parameters for current cosmic shear surveys. Baryons, however, need to be accounted for at the percent level for future Stage IV surveys and we propose to include uncertainty on the AGN feedback amplitude by marginalizing over this parameter using multiple simulations such as those presented in this paper. Finally, we explore the possibility of mitigating the impact of baryons by filtering the $M_{\rm ap}$ map but find that this process would require to suppress the small-scale information to a point where the constraints would no longer be competitive.


Introduction
Weak lensing cosmic shear is one of the most powerful cosmological probe of the late-time Universe. So far, most analyses have focused on studying the correlation between the shape distortions of pairs of galaxies as a function of their separation: the so-called shear two-point correlation function (γ-2PCF; e.g., Troxel et al. 2018;Hikage et al. 2019;Asgari et al. 2020). However, it is becoming clear that other estimators, and in particular those based on weak-lensing mass maps, outperform the standard γ-2PCF (e.g. Dietrich & Hartlap 2010;Ajani et al. 2020;Zürcher et al. 2020;Coulton et al. 2020;Martinet et al. 2020). Indeed, mass maps are highly sensitive to the non-Gaussian part of the matter distribution which arises from the non-linear growth of structures, which contains information that is overlooked by two-point estimators. Consequently the combination of both probes yields even tighter constraints as seen in recent applications to observational data (e.g., Martinet et al. 2018;Harnois-Déraps et al. 2020). With future Stage IV surveys, this combination is also expected to improve not only our measurement of the growth of structure parameter S 8 by a factor of two, but also that of the dark energy equation of state w 0 ) and of the sum of neutrino masses Σm ν (Li et al. 2019), by factors of three and two, respectively.
These non-Gaussian estimators are however difficult to predict theoretically due to limits in our understanding of the nonlinear growth of structures (see e.g., Fan et al. 2010;Lin & Kilbinger 2015;Shan et al. 2018;Giocoli et al. 2018b; Barthelemy et al. 2020, for some attempts) and are instead modeled with N-body simulations. This can significantly increase the computational cost of such analyses, but resorting to N-body simulations is also necessary to accurately model the γ-2PCF at scales affected by non-linearities (e.g., Euclid Collaboration: Knabenhans et al. 2020). Moreover, many public simulation suites can be exploited for this purpose, including the Scinet LIght-Cone Simulations (Harnois-Déraps et al. 2015, SLICS), the cosmo-SLICS (Harnois-Déraps et al. 2019) or the MassiveNuS .  Fabjan et al. 2010; and star formation (Springel & Hernquist 2003). Both effects redistribute the matter in and around dark matter haloes, however the exact amplitude of the feedback varies between simulations (see Chisari et al. 2019, for a recent review on feedback in cosmological hydrodynamical simulations). Extending the simulation suites of Martinet et al. (2020) based on the SLICS and the cosmo-SLICS, we construct Euclid-like mocks from the Magneticum hydrodynamical simulation. We measure the impact of baryons on the γ-2PCF, on the lensing peaks, minima, and on the lensing PDF, in the form of a multiplicative baryon bias correction factor. We next propagate this correction into the cosmological inference pipeline described in Martinet et al. (2020), and investigate the impact on the parameter forecasts. Finally, we explore various mitigation schemes to decrease the baryondominated small-scale contribution to the M ap computation.
We introduce the Magneticum simulation and compare it to other state-of-the-art hydrodynamical simulations in Sect. 2. We present the methodology that we employ to measure the bias due to baryons in Sect. 3. We then measure their impact on different data vector (DV) in Sect. 4.1 and propagate the effect to the cosmological forecasts in Sect. 4.2. We explicit different mitigation setups in Sect. 5 and conclude in Sect. 6. Finally, we adapt our mocks to the KiDS and DES surveys in Appendix A and measure the impact of baryons on the cosmological constraints by Martinet et al. (2018); Harnois-Déraps et al. (2020) in these two surveys.

The Magneticum hydrodynamical simulation
The Magneticum suite (Biffi et al. 2013;Saro et al. 2014;Steinborn et al. 2015Steinborn et al. , 2016Dolag et al. 2015Dolag et al. , 2016Teklu et al. 2015;Remus et al. 2017a;Castro et al. 2020) is a compilation of N-body and hydrodynamical simulations describing the cosmic evolution of the Universe. In total, the suite follows up to 2×10 11 particles divided in dark matter, gas, stars, and black holes. The simulations were performed with the TreePM+SPH code P-Gadget3 -a more performing and efficient version of the publicly available Gadget-2 code (Springel 2005) developed con- Table 1. Subset of the Magneticum simulation suite used in this work. From left to right: box size, gravitational softening and the particle masses for the different components (DM, gas, and stars), the number of lens planes built, redshift range of the past-light cone, its field of view, and the map angular resolution. A "−" indicates that the parameter value of Box 2b/hr is identical to that of Box 2/hr. comitantly to its successor Gadget-4 (Springel et al. 2020). Notably, the smoothed-particle-hydrodynamics (SPH) solver implements the improved model of Beck et al. (2016). The particle dynamics is coupled to different astrophysical effects such as radiative cooling, heating by a uniform evolving UV background, star formation (Springel & Hernquist 2003), stellar evolution and chemical enrichment processes (Tornatore et al. 2007). Cooling is implemented following the metallicity dependent formulation presented in Wiersma et al. (2009), using cooling tables produced by the publicly available CLOUDY photo-ionization code (Ferland et al. 1998). Lastly, AGN feedback and black hole growth are modeled as described in . The sub-grid physics model of our simulation reproduces a long list of observations, from galaxy properties (Teklu et al. 2015;Remus et al. 2017b) and AGN population Steinborn et al. 2016), to the inter-galactic and intercluster medium Bocquet et al. 2016;Gupta et al. 2017). Of special attention for this work is the robustness of the AGN feedback of our model, which has a clear footprint on the matter power-spectrum, suppressing its amplitude with respect to the DM-only case on k ∼ 10 h Mpc −1 by ∼ 15%. The specific suppression range and amplitude depend strongly on the AGN model. In Fig. 1 we compare the ratio of the matter power-spectrum computed from Box 2 (see Table 2 (Pillepich et al. 2018), and OWLS (Schaye et al. 2010). As can be observed, Magneticum provides an AGN feedback suppression consistent with other simulations, in particular with BA-HAMAS and in a lesser extend with TNG-300 that are used in other recent mass map analyses. We note that the AGN feedback also controls the gas fraction of halos (and the conversion efficiency into stars) such that the simulations presented in Fig. 1 also differ in these quantities although the matter power spectrum is a good indicator to assess the current theoretical uncertainties in modeling baryonic physics.

Past light-cone reconstruction
The Magneticum past light-cone reconstructions are performed using the Simulation LIght conE buildeR code -SLICER 2 -closely following the pipeline presented in Castro et al. (2018) and forked from a MapSim branch (Giocoli et al. 2015(Giocoli et al. , 2018aHilbert et al. 2020). Shortly, light-cones are built in post-processing, assigning particles to predetermined 2D mass maps according to the triangular-shaped cloud mass assignment scheme. The geometry of the past-light-cones is a square based pyramid (in angular coordinates) where the observer is located at the z = 0 vertex, and which extends to z max . The opening angle is chosen to be 10 degrees and the angular resolution of the 2 https://github.com/TiagoBsCastro/SLICER light-cone mass planes is 3.6 arcsec. SLICER allows mass maps thicknesses that are a rational fraction of the box size, which we have chosen to be half of the box size. Every two mass maps, particles were randomly shifted and reflected with respect to one of the box axis (accounting for periodic boundary conditions) in order to avoid the repetition of the same cosmic structure along the line of sight.
Next, mass maps are converted into surface density maps Σ(x, y) as where n indicates the number of particles, m j the interpolated contribution of the j th particle to the pixel at position x, y and L p the physical size of the pixel in units of Mpc h −1 . Given a source plane, convergence maps κ(x, y) are created by weighting these maps by the corresponding critical surface mass density and integrating over the past light-cone. Here, c indicates the speed of light, G the Newton's constant and D l , D s and D ls the angular diameter distances between observer-lens, observersource, and lens-source, respectively. The shear components (γ 1 , γ 2 ) are obtained from the standard inversion technique of Kaiser & Squires (1993). Since Box 2b Hydro has not been run down to z = 0, our past light-cones are built from grafting Box 2 in the range z = [0, 0.248] with Box 2b for z = [0.248, 3.44]. We have used the DM-only runs to validate our approach, and verified that for the source redshift distribution of interest in this paper, grafting affects only marginally the lensing statistics for z s 1.0 with respect to the same statistics inferred from a contiguous lightcone 3 . In total we have produced 20 pseudo-independent lightcones -10 Hydro and 10 DM-only -with properties summarized in

Methodology
In this section, we first give a brief description of our cosmological forecasting pipeline in Sect. 3.1, then detail in Sect. 3.2 the calculation of baryonic bias.

From shear to cosmology
Cosmological forecasts are computed with the same pipeline as in Martinet et al. (2020). While we recapitulate the salient points of this analysis here, we refer the reader to this publication for more details. • Measurements: From each mock, we compute a 1024 × 1024 pixels M ap map (Schneider 1996) as a convolution between the tangential ellipticity t around a pixel θ 0 , and the Schirmer et al. (2007) compensated Q filter adapted to the detection of matter halos: with This filter is a simpler form of the NFW (Navarro et al. 1997) profile with an additional exponential attenuation in the center and edge of the aperture. x c controls the tilt between the core and the edge of the profile and we introduce an optional inner radius parameter θ in to govern the exponential cut off in the inner part. The sum in Eq. (3) is carried out over galaxies at position θ i within an aperture of radius θ ap = 10 by default (corresponding to an effective smoothing scale of θ ap x c = 1.5 ) and centered on θ 0 , and it is normalized by the galaxy density n gal within the aperture. We estimate the noise due to intrinsic ellipticities in the mocks to define several M ap map statistics based on the signal-to-noise ratio (S/N) distribution of pixels of these maps: peaks and voids (pixels with values greater or smaller than their 8 neighbors), as well as the full distribution of pixels (1D M ap , often referred to as the lensing PDF). Following Martinet et al. (2020), these distributions are organised in 8 bins between −2.5 < S/N < 5.5 for peaks, 8 bins between −5 < S/N < 3 for voids, and in 9 bins between −4 < S/N < 5 for the 1D M ap . For each mock, we also compute the γ-2PCF ξ ± (ϑ) using the Athena software (Kilbinger et al. 2014), with the estimators defined as (Schneider et al. 2002): , where the sum is over galaxy pairs i j separated by an angle ϑ and potentially in different tomographic bins. The results are binned in 8 angular bins logarithmically spaced between 0.1 and 60.5 for ξ + and 0.5 and 300 for ξ − . The quantities t,× are the tangential and "cross" components of the ellipticities.
• Tomography: We separate the mock data into five redshift bins: 0 < z 1 < 0.47, 0.47 < z 2 < 0.72, 0.72 < z 3 < 0.96, 0.96 < z 4 < 1.33, and 1.33 < z 5 < 3. We reconstruct M ap maps from the galaxy samples in individual redshift slices (auto-M ap ) and from the combination of multiple redshift bins from two and up to five (the cross-M ap terms introduced in Martinet et al. 2020). The latter showed that including the cross-maps yields a significant improvement in the forecast precision. For the γ-2PCF, we include both the auto-and the cross-tomographic correlations.
• Likelihood: Finally, we compute the likelihood on a fourdimensional grid, taking the "fiducial" cosmo-SLICS simulations as our observations. We use a Student-t likelihood (Sellentin & Heavens 2016), which generalizes the multivariate Gaussian, and we account for the noise through the covariance matrix. We also improve the accuracy of the emulator used in Martinet et al. (2020), initially limited by the number of points in the interpolation grid, by performing a second interpolation in the parameter space restricted to the hyper-volume where the likelihood is non-zero. We measure the 1σ uncertainty on each cosmological parameter by finding the range of values enclosing 68% of the likelihood previously marginalized over all other parameters. Since the posterior distribution on h extends up to the prior limit imposed by the simulation space, we only discuss the forecasts for the other three probed parameters 4 : S 8 , w 0 , and Ω m . All predictions are computed for a 100 deg 2 area to ensure that the few percent uncertainties in the model (primarily due to the accuracy of the N-body simulations) are lower than the precision of the forecasts.

Measuring biases
In this article, we use the pipeline described above to study the impact of the baryonic bias on a Stage-IV cosmic shear data analysis. We note however that we could propagate instead biases due to shear measurement uncertainty, mean photometric redshift inaccuracy, galaxy intrinsic alignments or source-lens coupling (see e.g., Kacprzak et al. 2016;Harnois-Déraps et al. 2020), which we leave to future work. The strategy is to bias the observation DV and compare the positions of the maximum of the likelihood to the no-bias case. We use the DM-only and equivalent hydrodynamical runs of the Magneticum simulations to measure the cosmological bias due to baryonic effects. We create galaxy mocks that reproduce the same Stage-IV survey properties as the cosmo-SLICS. In particular, galaxy redshifts, positions, and intrinsic ellipticities are identical to that of the model. We generate 50 realizations of the hydro and DM light-cones: 10 different line-of-sights to lower the sample variance, each populated with 5 different realizations of intrinsic ellipticities to converge on an average shape noise contribution. We compute the M ap map in every mock, extract the DV, and measure the ratio between the average over the 50 DVs measured in the hydrodynamical and in the DM-only mocks.
This multiplicative correction factor is computed for each S/N bin and serves to infuse our Magneticum baryon model in our (DM-only) observation data. Infusing baryons into the model avoids being affected by small residual differences between the

Impact of baryons
We examine in this section the impact of baryons on the different DVs (Sect. 4.1) and on the cosmology inferred by a Stage-IV survey (Sect. 4.2).

Impact on computed statistics
We show in Fig. 2 the impact of baryons on the γ-2PCF, presented as the fractional difference between ξ ± measured in the hydrodynamical and in the equivalent DM-only runs. This figure shows the baryonic bias in absence of tomography to better resolve the effect, however similar curves are observed in the tomographic case. We measure a decrease in the amplitude of ξ ± of up to 10 − 15% at scales below a few arcmin, which is fully consistent with the suppression in the matter power spectrum seen in Fig. 1.
In Fig. 3 we show the equivalent measurement on the M ap statistics: voids, peaks, and 1D M ap . The impact of baryons on these DVs is more complex than in the case of the γ-2PCF given the different physical origins of each part of the M ap distributions. The most plausible scenario is described in Osato et al. (2020): Overdense regions are diluted due to AGN feedback, leading to smaller amplitude of high-S/N structures. This expelled material can be deposited in low density regions, which likely explains the reduction of pixels with highly negative S/N. We note that the cause of the latter effect is not fully understood yet and could involve other baryonic processes such as interaction of internal and accretion shocks (e.g., Zhang et al. 2020). The increase of the distribution at S/N close to zero accounts for the density smoothing due to this redistribution of matter which leads to a higher number of small S/N structures. This reasoning is deducted from the lensing PDF behavior but holds for peaks and voids. We do not probe the impact of radiative cooling seen at very high S/N in Osato et al. (2020) as we focus here on a lower, more conservative S/N range. This would also necessitate to use a pixel scale finer than our fiducial 0.59 . Indeed these effects appear to be significant only at scales lower than 0.5 in the M ap map according to the location of the ξ + upturn in Fig. 2 and consistent with the propagation of the radiative cooling scale of ∼ 15 h Mpc −1 into the lensing PDF using Eq. (11) of Castro et al. (2018).
Quantitatively, we find a decrease of 5% to 13% in the extremal values of our DVs, depending on the considered DV and S/N bin. In particular, for peaks of S /N = 4 we measure a depletion of ∼ 5% due to baryons, in perfect agreement with measurements from Fong et al. (2019) and Coulton et al. (2020) on the BAHAMAS simulations, which present a similar amplitude of the baryonic feedback (see Fig. 1). This is also of the same order of magnitude as the results from the TNG simulations ) and the baryonification method (Weiss et al. 2019).
When applying tomography, we also note a decrease of the impact of baryons in each slice compared to the combined case. Although Osato et al. (2020) also found, by using thin source slices, that the baryon bias is lower at higher redshift, the lower effect in our case is more likely due to an increase of the noise due to lower galaxy densities as noted by Harnois-Déraps et al.
(2020) in their tomographic peak counts analysis of DES data. This is supported by the fact that we find similar baryonic effects in all five redshift slices, which we designed to have identical galaxy densities.
In Figs. 2 & 3 the error bars are computed for the expected 15 000 deg 2 of the Euclid survey by area-rescaling the SLICS covariance matrix as Cov → Cov (100/15 000). These figures show that the impact of baryons on all tested estimators will be significant with respect to the statistical precision of Stage IV cosmic shear surveys. For the present 100 deg 2 analysis, the changes are below the statistical noise, whence the necessity to

Propagation to cosmological constraints
We propagate the impact of baryons on the cosmological parameter forecasts following the method described in Sect. 3. We focus on the peak statistics as it is the most widely used mass map estimator in the literature, but we also presents results for voids and lensing PDF. All values reported in this section are summarized in Table 2. Fig. 4 presents our results on peaks, first without including tomography. We show the 1 and 2σ contours of the marginalized 2D and 1D likelihoods for S 8 , Ω m , and w 0 . The blue contours and curves correspond to the forecasts for the DM-only observations, while the data yielding the green ones are infused with the baryon bias. We see a shift towards smaller values of S 8 (∆S 8 = S DM+Baryons 8,best − S DM 8,best = −0.028 (−3.4%)) and Ω m (∆Ω m = −0.018 (−6.1%)) when including baryons. As expected, including baryonic feedback mimics the effect of having lower matter density/structure growth, the two being highly degenerate for cosmic shear. These results are fully consistent with the Ω m shift recently found by Coulton et al. (2020) for peaks measured in the BAHAMAS simulations, but larger than the earlier work of Osato et al. (2015), who reported a −1.5% shift in Ω m from simulations with weaker AGN feedback. Finally, we see negligible changes in w 0 , in part because the sensitivity to that parameter is quite low in the non-tomographic case.
When including tomography (see Fig. 5), the constraining power increases significantly. We measure a very similar effect as in the no-tomography case for S 8 and Ω m but with slightly smaller shifts (∆S 8 = −0.024 (−2.9%) and ∆Ω m = −0.005 (−1.7%)). This behavior is also found in the DES-Y1 Magneticum mocks (see Harnois-Déraps et al. 2020, and Appendix A): with tomography, the noise in each mass map in-  creases and tends to wash out the effect of baryons. We also find a small positive shift in the maximum of the 1D likelihood for w 0 (∆w 0 = 0.035 (3.5%)). We note, however, that this result is only supported by a small distortion of the likelihood, which otherwise agrees fairly well with the DM-only case. Although this effect could be physically motivated it is likely due to degeneracies in the parameter space notably between S 8 and w 0 , and it will be interesting to see if it persists in future analyses.
If we now combine peaks with the γ-2PCF in the tomographic setup (Fig. 6), the results are similar to those from peaks alone, but accentuated in amplitude. We find ∆S 8 = −0.037 (−4.4%), ∆Ω m = −0.019 (−6.5%), and ∆w 0 = 0.047 (4.7%). Both peaks and γ-2PCF are affected in a similar manner by baryons and thus their combination show a larger effect. Although one could have hoped that the impact of baryons would diminish when adding the information of the γ-2PCF that partly come from larger scales, this study demonstrates on contrary that baryonic physics do not vanish in this combination and need to be accounted for in future surveys.   Table 2 shows the forecast biases for all the different tested DVs and their combination with the γ-2PCF. Overall baryons impact the different M ap estimators in a similar manner as expected from the similarities in how they affect each DV in Fig. 3. We report smaller biases for voids than for peaks, as already noted in Coulton et al. (2020) in the no-tomography case, and confirm this trend when including tomography, with possibly a very low bias on w 0 . When combined with the γ-2PCF we also note a difference of sign in the w 0 bias which could highlight a different sensitivity of voids to this parameter, but is more likely due to a residual small peak in the likelihood from the interpolation of the DV in this case. 1D M ap presents a comparable bias to peaks with slightly lower shifts as well. Although voids are the less affected by baryons, this analysis shows that all estimators present biases of a few percents on at least one of the probed cosmological parameters. Considering the degeneracies between cosmological parameters, this highlights the necessity of accounting for baryons when modeling the dependence of non-Gaussian M ap estimators on cosmology. Because these results depend on the particular implementation of baryonic feedback processes in the Magneticum simulation, we run an additional test where we only infuse half of the DV baryonic bias to compute the cosmological forecasts. Although we cannot accurately model the response of the M ap statistics to AGN feedback with only one simulation, this case likely corresponds to a much lower feedback amplitude. We find a bias on cosmological parameters which approaches the percent value in the tomographic case. This strong dependency of the cosmological parameters on the amplitude of the infused baryonic effect suggests that the impact of baryons could be mitigated by integrating a modeling of AGN feedback in the likelihood, a possibility which is not studied here.
Although these few percent biases are worrisome for future Stage IV cosmic shear surveys, we note that they remain fairly small compared to the statistical precision of current surveys. Our 100 deg 2 mocks at Euclid depth include about ten million galaxies, a similar number to Stage III surveys. In this case, the biases are always below 1σ for every parameter and configuration, except for the S 8 parameter when combining M ap estimators and γ-2PCF where the bias can reach up to 2σ. In Appendix A, we tailor our mocks to the KiDS-450 and DES-Y1 surveys to verify the impact of baryons on the peak statistics analyses conducted in Martinet et al. (2018) and Harnois-Déraps et al. (2020), respectively. Our findings validate the choice of neglecting baryons in current Stage III peak counts analyses.

Mitigating baryons with small scale cut
In the case of γ-2PCF, the impact of baryons can be mitigated by discarding small scales which are the most affected, as seen from Fig. 2. Although this process decreases the statistical precision as it removes part of the signal, it allows one to gain on accuracy without needing to run computationally expensive hydrodynamical simulations. This trade off was notably chosen by the DES collaboration in the first year data release (Troxel et al. 2018). Alternatively, baryons can be modeled with halo-based codes (HMCode, Mead et al. 2020) or from libraries of power spectra (van Daalen et al. 2020), allowing the inclusion of smaller angular scales and increasing the statistical power as in Hildebrandt et al. (2017) and Huang et al. (2020).
Since we do not yet have access to a library of hydrodynamical simulations to estimate the baryonic bias, we follow the work of Weiss et al. (2019) and explore different possibilities to mitigate the effect of baryons on M ap statistics. As suggested by their analysis we vary the size of the aperture θ ap but we also investigate possible variations to the Schirmer et al. (2007) Q filter shape entering Eq. (3). In particular, we vary the aperture filter size θ ap , the tilt parameter x c and the inner filter radius parameter θ in .
The top part of Fig. 7 shows the impact of varying θ ap , θ in , and x c on the shape of the Q filter. The fiducial values that we use in the rest of the paper are θ ap = 10 , θ in = 0.4 , and x c = 0.15. In this figure, we only present results when varying one parameter at a time for θ ap = 20 , θ in = 2.4 , and x c = 0.45 to better highlight the effect of each parameter. We however investigate multiple values in the range 3.3 ≤ θ ap ≤ 106.7 , 0.4 ≤ θ in ≤ 3 , and 0.05 ≤ x c ≤ 1 and explore variations of all three parameters together. Increasing any of these three parameter values reduces the relative importance of the small scales in different manners. Increasing θ ap adds galaxies further away from the aperture center; increasing θ in removes galaxies close to the aperture center without distorting the general weighting; finally, increasing x c up-weights distant galaxies without modifying the number of galaxies captured by the aperture. We re-compute the baryonic bias from M ap constructed with the new Q filters and observe, in the bottom part of Fig. 7, that none of these methods is completely efficient at removing the impact of baryons. This is because with M ap statistics we cannot discard close galaxies as is done for γ-2PCF: we only reduce their weight when centered on them which does not prevent them from entering other apertures on a θ ap scale. We note a mild improvement when increasing θ in and x c . Increasing θ ap to 20.0 still decreases the impact of baryons by ∼ 50% at the extreme S/N values, and only the A&A proofs: manuscript no. spe2 largest θ ap > 100 can bring it to zero, but the cost in precision is high, as we show next.
We apply these filter modifications to the measurements from the cosmology and covariance mocks, and carry out for each of these a full cosmological forecasts. We present the variations of the bias on the inferred cosmology for these different configurations in Table 3, relative to the fiducial Q filter. We also report the change in the forecast precision due to the reduction in smallscale information. We show these results for peaks in a configuration without tomography but we found similar behavior for other M ap estimators and including tomography: the bias due to baryons is reduced (e.g. ∆S 8 /∆S Q fid. 8 < 1.0), but at the cost of a reduction in the statistical precision (δS 8 /δS Q fid. 8 > 1.0). Quantitatively, we find with θ ap = 20 that the bias is reduced by a factor of almost two and three on S 8 and Ω m respectively, but with a loss of respectively 12% and 17% on the forecast precision, and a loss of 15% on w 0 . The two other variations are less efficient but also retain more of the cosmological information. The cut at θ in = 2.4 decreases the bias by ∼ 30% at the cost of a 5%, ∼ 10%, and ∼ 20% wider statistical precision on S 8 , w 0 , and Ω m , respectively. We find similar results for all the Q filter con-figurations using higher parameter values than the fiducial: the gain in accuracy is always balanced by a significant loss in precision. When using smaller values of θ ap , θ in , or x c , the impact of baryons is however increased because of the larger contribution of the small-scale baryonic features but the constraining power is also degraded as we chose the fiducial Q filter to maximize the forecast precision in Martinet et al. (2020).
Overall, none of the small-scale cuts we applied to mitigate baryonic effects are able to decrease the bias below a few percent accuracy while preserving a strong statistical precision. Using the effective scale of θ ap x c = 16 recommended in the analysis of Weiss et al. (2019) for Euclid-like mocks, we confirm that the bias becomes consistent with zero. However, such large smoothing scale results in a decrease of the constraining power by a factor of more than 2, 1.5, and 3 on S 8 , w 0 , and Ω m , respectively, motivating a full forward-model approach of the baryonic bias.

Conclusion
In this paper we investigate the impact of baryons on various M ap statistics: peaks, voids, and the lensing PDF. Baryonic physics is modeled with the state-of-the-art Magneticum hydrodynamical simulations, and its impact on the data vector is propagated into a full cosmological forecasts on S 8 , w 0 , and Ω m , for a Stage-IV lensing survey. The likelihood sampling exploits the cosmological pipeline of Martinet et al. (2020), which is based on the SLICS and cosmo-SLICS DM-only simulations. Our results are summarized below: • Baryons are biasing the measured M ap estimators by about 5 − 10% on most S/N bins, notably decreasing the number counts at extreme S/N values, while increasing the number of intermediate S/N features. This is a direct consequence of strong baryonic feedback, which dilutes the density profile of massive halos and decreases their S/N in the M ap map.
• In our Stage-IV survey setup without tomography, the baryonic feedback propagates into a negative bias of about −3% on S 8 for every estimator. The bias on Ω m depends on the estimator and ranges from zero for voids to −6% for peaks. When including a tomographic decomposition with five redshift slices including cross-tomographic bins, these biases are slightly lowered, likely due to the increased shape noise in each tomographic slice, but remain of the order of a few percents. We observe in our tomographic setup positive bias of the order of 3 − 5% on w 0 , however it is not clear at the moment whether this has a physical origin or is caused by parameter degeneracies that are not fully captured by our likelihood.
• Biases on all parameters are increased to ∼ 5% when combining any M ap statistics with γ-2PCF in the tomographic analysis. This combined analysis is maximally affected as the contributions of baryons are in the same direction and of similar amplitude for the individual probes.
• After investigating a range of scale cuts on the M ap statistics, we find that it is difficult to efficiently lower the impact of baryons without significantly degrading the statistical power. In line with Weiss et al. (2019), we find that only an overly large aperture size could lower the bias to sub-percent level. A large portion of the signal is lost with such filtering, leading to less competitive cosmological constraints. Table 3. Changes in the measured bias due to baryons and in the associated forecast precision for various mitigation schemes in the case of M ap peaks without tomography. The comparison is performed with respect to the fiducial Q filter used in the rest of the paper with θ ap = 10 , θ in = 0.4 , and x c = 0.15. We vary one parameter at a time, the two others being held to their fiducial values. "∆" refers to the bias due to baryons, and "δ" to the 1σ precision forecast. We built Magneticum mocks to measure the impact of baryons on peak statistics analyses of current stage III surveys, namely in KiDS-450 (Martinet et al. 2018) and DES-Y1 , and show that it remains below the statistical uncertainties associated to these surveys. This will, however, not be the case for future Stage IV surveys for which baryons need to be accounted for in order to reach a percent-level precision while remaining unbiased.
In this article, we present a correction scheme based on the DV: we measure a corrective factor from hydrodynamical simulations and apply it to mock observations in order to estimate how biased would become cosmological constraints in case this step was omitted. We note, however, that these results fully depend on the amplitude of the baryonic feedback modeled by the Magneticum simulations, which is still uncertain as seen from the scatter between the different state-of-the-art simulations. A more accurate correction would consist in modeling the impact of baryons on M ap statistics from a set of hydrodynamical simulations with various feedback amplitude values, and to marginalize over the extra free parameters when computing the cosmological constraints. Coulton et al. (2020) recently showed the feasibility of this approach using the BAHAMAS simulations ran with three different amplitudes of the AGN feedback. The baryonification method described in Schneider & Teyssier (2015) particularly suits the needs of such analysis and is thus a promising tool to design future sets of N-body simulations that explore both the cosmology dependence and the response to baryons of non-Gaussian statistics. tions. A total of 2600 mock surveys are used for the wCDM peak count model, and 1240 for the covariance matrix. The impact of baryons on the tomographic constraints from M ap peak statistics in DES-Y1 are shown in Fig. A.2. The blue contours correspond to the 1 and 2σ constraints on the S 8 and Ω m parameters after marginalisation over the photometric redshift and the shear calibration uncertainties, while the black contours further include the correction due to baryonic physics. As expected, we find again a positive shift towards larger values of S 8 and Ω m . Quantitatively, the shift is of ∆S 8 = 0.013 (1.8%), lower than the 2.9% found with the Euclid-like mocks, again due to the lower galaxy density in the DES-Y1 data, which is ∼ 5 arcmin −2 . In terms of statistical precision, this bias corresponds to a 0.32σ shift and can be safely ignored in the DES-Y1 analysis. Harnois-Déraps et al. (2020) explored other sources of systematics in their analysis, and show in their Figure 17 that baryons and intrinsic alignments are the two most important effects, while source-lens clustering and uncertainty in the nonlinear growth of structure could be safely ignored.