Effects of baryons on weak lensing peak statistics

TL;DR

This study examines how baryonic feedback, especially from AGN, impacts weak lensing peak statistics in convergence maps, highlighting the importance of smoothing scale choices for future surveys like DES and Euclid.

Contribution

It provides a detailed analysis of baryonic effects on weak lensing peaks using simulations and the baryonic correction model, emphasizing the dependence on smoothing scale.

Findings

Baryonic effects significantly depend on the smoothing scale applied.

A smoothing of 8 arcmin suffices for DES-like surveys to mitigate baryonic impact.

A smoothing of 16 arcmin is needed for Euclid-like surveys to control baryonic suppression.

Abstract

Upcoming weak-lensing surveys have the potential to become leading cosmological probes provided all systematic effects are under control. Recently, the ejection of gas due to feedback energy from active galactic nuclei (AGN) has been identified as major source of uncertainty, challenging the success of future weak-lensing probes in terms of cosmology. In this paper we investigate the effects of baryons on the number of weak-lensing peaks in the convergence field. Our analysis is based on full-sky convergence maps constructed via light-cones from $N$-body simulations, and we rely on the baryonic correction model of Schneider et al. (2019) to model the baryonic effects on the density field. As a result we find that the baryonic effects strongly depend on the Gaussian smoothing applied to the convergence map. For a DES-like survey setup, a smoothing of $\theta_k\gtrsim8$ arcmin is sufficient to keep the baryon signal below the expected statistical error. Smaller smoothing scales lead to a significant suppression of high peaks (with signal-to-noise above 2), while lower peaks are not affected. The situation is more severe for a Euclid-like setup, where a smoothing of $\theta_k\gtrsim16$ arcmin is required to keep the baryonic suppression signal below the statistical error. Smaller smoothing scales require a full modelling of baryonic effects since both low and high peaks are strongly affected by baryonic feedback.

Tables

Table 1. Table 1 : Set of cosmological parameters used in all our N 𝑁 N -body simulations.

$h$	$n_s$	${\Omega }_{\Lambda }$	${\Omega }_M$	${\Omega }_b$	${\sigma }_8$
0.67	0.96	0.68	0.32	0.045	0.83

Table 2. Table 2 : Baryonic benchmark models with parameters ( M c subscript 𝑀 𝑐 M_{c} , μ 𝜇 \mu ) obtained by fitting the baryon fraction of the BC model to X-ray observations. Models A, B, and C are based on different assumptions regarding the hydrostatic-mass bias ( b hse subscript 𝑏 hse b_{\rm hse} ). The gas ejection radius has been fixed to θ ej = 4 subscript 𝜃 ej 4 \theta_{\rm ej}=4 in agreement with observed gas profiles.

Name	$1-b_{\mathrm{hse}}$	${\theta }_{\mathrm{ej}}$	$M_c$ [M⊙/h]	$\mu$	Baryonic scenario
Model A	1.000	4	$2.3\times {10}^{13}$	0.31	optimistic (weak baryonic effects)
Model B	0.833	4	$6.6\times {10}^{13}$	0.21	best-guess (medium baryonic effects)
Model C	0.714	4	$1.9\times {10}^{14}$	0.17	pessimistic (strong baryonic effects)

Equations

$$\kappa (\hat{\theta })=\frac{3}{2}\frac{H_0^2{\Omega }_m}{c^2}\sum _b{W}_b\left(\frac{12\cdot {\text{NSIDE}}^2}{4\pi }\frac{L^3}{N}\frac{n_b(\hat{\theta })}{\chi (z_b{)}^2}-\left[\chi (z_b+\Delta z_b)-\chi (z_b-\Delta z_b)\right]\right),$$

$$W_{b}=\left(\int_{\Delta z_{b}}\frac{\text{d}z}{E(z)}\frac{\chi(z)}{a(z)}g(z)\right)\Big{/}\left(\int_{\Delta z_{b}}\frac{\text{d}z}{E(z)}\int_{z}^{z_{H}}\text{d}z^{\prime}n_{s}(z^{\prime})\right),$$

$$g(z)=\frac{3}{2}{\Omega }_m{\left(\frac{H_0}{c}\right)}^2(1+z){\int }_z^{z_{\text{max}}}{n}_{\mathrm{s}}(z)\frac{\chi (z^{\prime })-\chi (z)}{{\chi }^{\prime }(z)}d{z}^{\prime }$$

$$n_s(z)\propto z^2\text{exp}\left(-\frac{z}{0.24}\right)$$

$$⟨{\sigma }_{\kappa ,\text{pix}}^2⟩=\frac{{\sigma }_e^2}{A_{\text{pix}}{n}_{\text{gal}}},$$

$$⟨{\sigma }_{\kappa }^2⟩=\frac{{\sigma }_e^2}{4\pi {\theta }_{\kappa }^2{n}_{\text{gal}}},$$

$${\theta }_{\kappa }\in \{2,4,8,16\}\text{arcmin}.$$

$$P(k,z)\propto ⟨\tilde{\delta }(k,z){\tilde{\delta }}^{\ast }(k,z)⟩,$$

$$C(l)=\frac{1}{2l+1}\sum _{m=-l}^l{q}_{lm}{q}_{lm}^{\ast },$$

$$C(l)={\int }_0^{\chi (z_{\text{max}})}{g}^2(\chi (z))P\left(\frac{l}{\chi },z(\chi )\right)d\chi .$$

$$\mathcal{S}/\mathcal{N}=\frac{\kappa }{\sqrt{⟨{\sigma }^2⟩}},$$

$$\Delta n_{\text{peaks}}=n_{\text{peaks}}-n_{\text{peaks,noise}},$$

$${\rho }_{\mathrm{halo}}(r)={\rho }_{\mathrm{gas}}(r)+{\rho }_{\mathrm{cga}}(r)+{\rho }_{\mathrm{clm}}(r),$$

$${\rho }_{\mathrm{gas}}(r)\propto {\left[1+\frac{r}{r_{\mathrm{co}}}\right]}^{-\beta }{\left[1+{\left(\frac{r}{r_{\mathrm{ej}}}\right)}^2\right]}^{(7-\beta )/2},$$

$$\beta =3-{\left(\frac{M_{\mathrm{c}}}{M_{200}}\right)}^{\mu },$$

$$1-b_{\mathrm{hse}}\equiv \frac{M_{500,\mathrm{hse}}}{M_{500}},$$

$$s_{\mathrm{bcm}}\equiv \left(n_{\mathrm{bcm}}-n_{\mathrm{dmo}}\right)/{\sigma }_{\mathrm{dmo}}.$$

Full text

1 Introduction

In the currently favoured cosmological model, $\Lambda$CDM, the energy density of the universe is dominated by a cosmological constant ($\Lambda$) and a cold dark matter (CDM) component. The most accurate constraints for this model are obtained from the cosmic microwave background signal which are measured using surveys such as Planck [2, 3].

Gravitational lensing – the deflection of light around clustered matter – offers an additional avenue to test the standard model of cosmology. A particularly promising probe is weak lensing, where slight distortions of galaxy shapes are used to determine the underlying matter density distribution including the invisible dark matter component. Due to the weak nature of this effect, large area surveys observing many millions to billions of galaxies are required, such as KiDS111http://kids.strw.leidenuniv.nl/, DES222https://www.darkenergysurvey.org/, HSC333https://hsc.mtk.nao.ac.jp/ssp/ or the planned Euclid444https://www.euclid-ec.org/, LSST555https://www.lsst.org/, and WFIRST666https://wfirst.gsfc.nasa.gov/ surveys [4, 5, 6, 7, 8].

In order to use weak-lensing surveys as cosmological probes, accurate predictions of the density distribution in the universe at both linear and nonlinear scales are required. While cosmological perturbation theories only describe the linear regime, full cosmological $N$-body simulations [9, 10, 11, 12] or semi-analytical methods including the halo model [13, 14, 15, 16] are able to describe the nonlinear behaviour of structure formation. However, these approaches only consider gravitational effects describing structure formation as a fully collisionless process.

Recent work based on hydrodynamical simulations has produced mounting evidence that a gravity-only approach is no longer sufficiently accurate to predict the weak-lensing signal. For example, Refs. [17, 18, 19, 20, 21, 22, 23] show that baryonic effects, specifically the energy input from active galactic nuclei (AGN), have an impact on the matter distribution at cosmologically relevant scales. Note, however, that baryonic feedback effects are not predicted from first principles in hydrodynamical simulations, but they are included as sub-grid effects using semi-analytical recipes. As a consequence, there is no quantitative agreement between different simulations in terms of clustering statistics. Differences in the implementation of sub-grid feedback effects lead to a spread of up to 30 percent in the predictions of the power spectrum for wave-modes between $k\sim 0.1$ and 10 h/Mpc [e.g. Ref. 24, 25].

An alternative approach to full hydrodynamical simulations are models that parametrise the effects of baryons on the gravitational clustering. These parametrisations are usually performed at the level of dark matter haloes and the different model parameters are constrained using simulations [26, 27, 28, 29, 16] or direct observations [30, 1]. The advantage of a model parametrisation is that the baryonic effects can be self-consistently included in methods of cosmological parameter inference.

The present paper relies on the baryonic correction (BC) model proposed in Refs. [30, 1]. Gravity-only $N$-body simulation outputs are directly perturbed by slightly displacing particles around halo centres according to a parametrised halo profile which includes a stellar, gas, and dark matter component. Compared to the halo model approach, the BC method is based on a simulated density map, accurately reproducing all the complexity of large-scale structures including triaxial haloes, voids, filaments, etc. As a result, the model is not restricted to 2-point statistics but allows to study all summary statistics describing the matter distribution in the universe.

During the last decade, weak lensing peak statistics have emerged as a simple but powerful complementary tool to extract cosmological information that is not accessible with the angular power spectrum alone [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. The complementary nature of peak counts is highlighted by the fact that a combined analysis including the angular power spectrum results in considerably tighter cosmological constraints than any of these statistics used on their own [34]. The origin of weak lensing peaks has been closely investigated in Refs. [42, 43]. From there we know that the highest peaks are typically caused by one single massive halo. Medium and lower peaks, on the other hand, are generated either by shape noise or by line-of-sight projections of multiple smaller haloes.

In this paper we provide the first detailed study of the AGN driven baryonic effects on the peaks of weak-lensing convergence maps. While former work by Refs. [44, 45] also looked at the baryonic effects on weak-lensing peaks, they either did not include AGN feedback effects at all [44] or relied on a very mild AGN feedback recipe [45], potentially underestimating the true baryonic depletion effects. The latter statement is based on the angular power spectra published in Refs. [44, 45] showing considerably weaker baryonic effects compared to recent hydrodynamical simulations [e.g. 17, 19].

In the present study, we construct the convergence maps from a light-cone based on outputs of $N$-body simulations (which have previously been baryonified using the baryonic correction model) assuming an extended redshift distribution for the source galaxies. Together with the baryonic effects, we investigate the role of the galaxy number density and smoothing scale on the weak-lensing peak statistics. We thereby consider two survey configurations based on DES and Euclid. This will give us information on how baryons influence both current and next-generation weak lensing surveys.

The paper is structured as follows: In Sec. 2 we summarise the method to generate convergence maps from N-body simulations and we discuss the weak-lensing peak counts based on gravity-only simulations without baryonic effects. Sec. 3 provides a brief introduction to the BC model including three benchmark models that cover realistic scenarios in terms of the expected baryonic suppression. In Sec. 4 we present the main results of the paper, i.e., the effects of baryons on the weak-lensing peak statistic. Finally, we conclude our work in Sec. 5. Further details about the box replication method and resolution effects are found in the appendices A and B.

2 Numerical Methods

In this section we describe the process by which we generate weak-lensing convergence maps using $N$-body simulations and we discuss the angular power spectrum and the peak count statistic for the gravity-only case without baryonic effects. The convergence maps are based on particle light-cones which we construct by replicating standard $N$-body simulation outputs at many redshifts. In total, we generate fifty full-sky convergence maps for each survey configuration based on a set of ten statistically independent $N$-body simulations. This is possible because we randomise the particle positions in the replication process, generating five different convergence maps for every $N$-body simulation. In order to emulate different types of surveys with various numer densities, we furthermore post-process the generated convergence maps using different noise levels and smoothing scales.

2.1 N-body Simulations

The $N$-body simulations are generated using the publicly available $N$-body code PKDGRAV3 [46, 47]. We assume a flat $\Lambda$CDM cosmology with best-fitting parameters from the Planck 2015 [3] release (see Table 1). The initial conditions are generated with the MUSIC code [48] assuming a transfer function based on the fitting function from Eisenstein and Hu [49].

We simulate sub-volumes of box-length $L=512\,\,\text{Mpc}/h$ each with $N=512^{3}$ particles which are then replicated to fill the full survey volume. The particle positions within a simulation volume are offset by a random vector and the volumes themselves rotated during the replication process in an attempt to reduce any repetitive patterns that may appear due to the large number of replications. The justification behind this replication and randomization process and its effects are discussed in more detail in Appendix A.

Given this configuration and set of cosmological parameters, our simulations have a particle mass resolution of $M_{\text{particle}}=9\times 10^{10}\,M_{\odot}$/h. For the purpose of detecting the dark matter haloes in our N-body simulations we use AMIGA Halo finder [50]. We only select haloes with over 100 particles which corresponds to a halo mass resolution of $M_{\text{halo}}=9\times 10^{12}\,M_{\odot}$/h. In Appendix B we show that our resolution is sufficient to obtain converged results for the weak-lensing peak count.

2.2 Convergence Maps

We use a similar procedure as in Refs. [51, 52] to generate weak-lensing convergence maps for a single non-tomographic redshift bin of the size $z=0.1-1.5$. The lightcone is constructed using 78 concentric spherical shells, with each shell being filled by the particles of the replicated simulation volume using the respective snapshot of the N-body simulation. The simulated sub-volumes are replicated to fill the full survey volume of $\sim 6\,\,\text{Gpc}/h$ using up to 1728 replications. The particles within each shell are projected onto a HEALPix map (of resolution $N_{\rm side}=4096$), weighted according to a realistic galaxy source distribution and summed up along the line-of-sight assuming the Born approximation. The equation for the convergence of a single pixel on the HEALPix map can be formulated as:

[TABLE]

where $n_{b}(\hat{\theta})$ is the particle count within the $b$-th shell at position $\hat{\theta}$ and $\Delta z_{b}$ is the shell width at redshift $z_{b}$. The shell weights $W_{b}$ are defined as

[TABLE]

where $\chi$ is the comoving distance. The lens efficiency $g(z)$ and the source distribution $n_{s}(z)$ are given by

[TABLE]

[TABLE]

assuming a maximum redshift of $z_{\text{max}}=1.5$. The galaxy source distribution in Eq. (2.4) is a reasonable approximation for modern galaxy surveys.

After generating the raw convergence maps, we have to further account for the noise that is introduced in the observation process. For a detailed discussion of the noise and smoothing in convergence maps we refer to Ref. [53]. The pixelization of the convergence field can be approximated as smoothing using a top-hat filter. The Gaussian noise distribution is then described by

[TABLE]

where $A_{\text{pix}}$ is the area of a single pixel, $\sigma_{e}$ is the root mean square of the intrinsic ellipticity dispersion of galaxies and $n_{\text{gal}}$ the galaxy number density. We add Gaussian smoothing to reduce the noise and to better mimic realistic convergence maps generated from galaxy surveys. The noise after Gaussian smoothing is given by

[TABLE]

where $\theta_{k}$ represents the smoothing scale. The choice of the smoothing scale determines the level of noise, with a trade-off that larger smoothing scales lead to a loss of small-scale information.

In order to quantify the expected difference beteween current and future galaxy surveys, we assume two separate survey configurations for all our convergence maps. We generate a measurement by first adding noise corresponding a given galaxy number density, before smoothing with an arbitrary smoothing scale. From each full-sky map, we then take multiple measurements by cutting out multiple (non-overlapping) circular patches, where the size of each patch corresponds to the survey area. Any additional systematic effects due to the geometry of the survey footprint are ignored. Since the smoothing and noise are added before any patches are created, this should not affect the measurements of peaks.

The first survey configuration, which we will refer to as DES-like, assumes an area of $5000\ \textrm{deg}^{2}$ and a galaxy number density of $n_{\text{gal}}=6\ \text{arcmin}^{-2}$, allowing us to cut six independent patches from each full-sky map. The second survey configuration, which we will refer to as Euclid-like, has an area of $20000\ \textrm{deg}^{2}$ and a galaxy number density of $n_{\text{gal}}=30\ \text{arcmin}^{-2}$. This setup only allows the use of 2 independent patches for each full-sky map.

For both configurations we use an ellipticity dispersion of $\sigma_{e}=0.3$. We then smooth the resulting noisy maps assuming four different smoothing scales

[TABLE]

2.3 Power Spectrum

Two-point statistics correspond to the standard way of measuring the cosmological clustering process. Regarding the particle outputs from $N$-body simulations, we rely on the matter power spectrum

[TABLE]

where $\tilde{\delta}(k,z)$ are the matter perturbations in Fourier space. For the case of two-dimensional convergence maps, we use the angular power spectrum

[TABLE]

with the coefficients $q_{lm}$ representing the spherical harmonic basis functions of the decomposed convergence map. These two measures are not independent and can be directly related with each other using the Limber approximation [54]:

[TABLE]

The Limber approximation allows us to evaluate the accuracy of our generated maps by comparing the angular power spectrum directly measured on the convergence maps with the one obtained via Eq. (2.10) using the matter power spectrum from the N-body simulations.

In the left panel of Fig. 1 we show the average power spectrum measured from all our simulations at $z=0.5$ and compare it to the Takahashi et al. [15] and Mead et al. [16] semi-analytical results. The agreement is within a few percent except at very large scales (small $k$-modes), where our simulations are missing power due to the limited box-size effect (see Appendix A for more information).

The panel at right-hand-side of Fig. 1 illustrates the angular power spectrum directly measured at the map level and compares it to the results from the Limber approximation based on the same underlying simulations. The two approaches show a very good agreement, validating our pipeline for generating weak-lensing convergence maps. The differences at lower multipoles can be attributed to the replication and randomization procedure which suppresses modes at scales comparable to the simulation box-length and above. We again refer to the Appendix for more details about resolution and replication effects, including a short discussion on why the systematical effects visible in Fig. 1 do not affect the main results of the paper.

2.4 Peaks

In the previous section, we used the angular power spectrum to validate our pipeline for generating weak-lensing maps from $N$-body simulations. As a next step we now turn our focus towards investigating peak statistics of the convergence field.

We use the definition of a peak as a single pixel on our convergence map that has a higher convergence ($\kappa$) value than any of its eight neighbouring pixels. It is convenient to quantify the peak height with respect to the signal-to-noise ratio

[TABLE]

where we assume a uniform Gaussian noise ($\sigma$) across the entire convergence map. The peaks are measured in a combination of linear and logarithmic bins. For the DES-like survey we use 12 linear (up to $\mathcal{S}/\mathcal{N}\approx 3.5$) and 6 logarithmic bins. The Euclid-like survey uses 20 linear (up to $\mathcal{S}/\mathcal{N}\approx 5.2$) and 8 logarithmic bins. The logarithmic binning at large $\mathcal{S}/\mathcal{N}$ values helps to reduce the statistical error at higher peaks, where the number of peaks drops off significantly. The outermost bin with the largest value of $\mathcal{S}/\mathcal{N}$ contains approximately 10 – 20 peaks for both the DES-like and Euclid-like setup.

Fig. 2 illustrates the peak density as a function of $\mathcal{S}/\mathcal{N}$ for both DES-like and Euclid-like survey configurations and all four smoothing scales. Also shown are the peak densities from pure noise maps. The latter dominates the lower peaks.

Instead of counting the total number of peaks, it can also be convenient to define the noise-subtracted peak density given by

[TABLE]

where $n_{\rm peaks,noise}$ is the peak-count from the noise-only map. This modified statistics allows us to immediately see where the peak counts are dominated by noise and how the cosmological signal compares to the size of the statistical error.

In Fig. 3 we show the noise subtracted peak density function for various smoothing scales and both survey configurations. The results based on DES-like survey are in good qualitative agreement with previous work [see Refs. 37, 38]. A comparison between the two survey configurations highlights the expected improvement in accuracy that can be expected from future surveys. The four-times larger survey area and the five-times better galaxy number density of the Euclid-like configuration compared to the DES-like setup results in a threefold increase of the signal in the Gaussian noise subtracted peaks. This means that future weak-lensing surveys will be capable of detecting small deviations of the peak distribution and therefore require a careful treatment of all systematics. This is especially true regarding the effects of baryons as we will see in the following sections.

3 Baryonic Correction Model

The baryonic correction (BC) model consists of a method to modify outputs of $N$-body simulations in order to mimic the effects of baryons on the cosmological density field. Here we provide a summary of the method and the parametrisation. More details can be found in ST15 [30] and S19 [1].

The BC model is based on an algorithm to radially displace particles around halo centres in $N$-body simulations. The goal is to slightly perturb individual halo profiles in order to account for baryonic effects without losing the information regarding the full cosmological density field. The displacement of particles relies on a parametrisation of the halo profile which includes both dark matter and baryonic components, and which is motivated by observations. The halo profile is given by

[TABLE]

where the different terms refer to the gas (gas), the central galaxy (cga), and the collisionless matter (clm) components. The latter consists of both dark matter (dm) and the satellite galaxies (sga) of the halo. In the following, we provide a short summary of the parametrisation:

• •

The gas profile ($\rho_{\rm gas}$) is described by a power law with inner core and outer truncation, motivated by stacked X-ray observations  [see for example Refs. 55, 56]. It is given by

[TABLE]

where the core and ejection radii are fixed to $r_{\rm co}=0.1\times r_{200}$ and $r_{\rm ej}=4\times r_{200}$, for simplicity. We refer to S19 for a more general parametrisation and a discussion about the effects of varying $r_{\rm co}$ and $r_{\rm ej}$ on the cosmological density field. The slope of the gas profile is given by the equation

[TABLE]

imposing a decreasing value of $\beta$ towards smaller halo masses, again motivated by X-ray observations [56]. The free model parameters $M_{c}$ and $\mu$ control the pivot scale and the slope of this mass dependence.

• •

The central galactic profile ($\rho_{\rm cga}$) is described in the baryonic correction model as a simple power law that is exponentially truncated beyond the half-light radius. Note, however, that the stellar profile only affects small scales and has therefore little importance for the weak-lensing signal.

• •

The collisionless matter profile ($\rho_{\rm clm}$) is given by a truncated NFW profile which is assumed to undergo adiabatic relaxation. The latter leads to a steepening of the profile in the centre and a flattening in the outskirts induced by the stellar and the gas components.

Note that the BC model only slightly perturbs the cosmological density field without destroying its complex structure, which includes voids, sheets, filaments, and haloes of triaxial shape. However, it does induce a flattening of profiles, a change in virial mass of haloes, and a slight displacement of halo positions in agreement with expectations from baryonic physics. The fact that the BC model only slightly perturbs the nonlinear matter distribution makes it much more realistic than more analytical approaches such as the halo model, where basic assumptions of randomly distributed halo positions or perfectly symmetric halo profiles have to be made.

The version of the BC model sketched out above has two free model parameters ($M_{c}$, $\mu$) regulating the amount of gas inside and outside of haloes. We now follow S19 and use observed gas fractions from X-ray data to constrain these parameters.

An important problem regarding gas fractions from X-ray observations is that they are derived assuming hydrostatic equilibrium of the gas. We account for this potential systematic by including a hydrostatic mass bias

[TABLE]

where $M_{500}$ is the actual halo mass at $r_{500}$, while $M_{\rm 500,hse}$ is the mass obtained from X-ray observations under the assumption of hydrostatic equilibrium. Results from hydrodynamical simulations suggest the hydrostatic bias ($b_{\rm hse}$) to be between 0 and 40 percent. Based on this, we define three benchmark models that correspond to an optimistic, best-guess, and pessimistic scenario with $1-b_{\rm hse}=1.0$, $0.83$, $0.71$, respectively.

The top panels of Fig. 4 show the gas fractions from a collection of X-ray observations [57, 58, 59] assuming the three different values for the hydrostatic mass bias mentioned above. This means that the first panel shows the uncorrected gas fractions from the literature, while the second and third panel is corrected to account for a hydrostatic mass bias of 20 and 40 percent, respectively. A large hydrostatic mass bias leads to an under-estimation of the total halo mass (see Eq. 3.4). Correcting for this causes the data points in the upper-central and upper-right panel to shift further downwards and to the right compared to the data points in the upper left panel. A more detailed explanation of the connection between X-ray gas fractions and the hydrostatic mass bias can be found in S19 [1]. Next to the X-ray data points, the top panels show the gas fractions obtained with the BC model where, the model parameters ($M_{c}$, $\mu$) are fitted to the data. This fitting procedure defines the BC models A, B, and C that we will used as benchmark scenarios throughout this paper. The characteristics of the benchmark models are summarised in Table 2. The values of the benchmark parameters are listed in Table 2.

The bottom panel of Fig. 4 shows the resulting power spectra of the three benchmark models. They feature a maximum suppression of about 15, 20, and 25 percent, repsectively. The curves agree with the power spectra of models A-avrg, B-avrg, and C-avrg from S19 to better than one percent for $k<5$ h/Mpc. The reason why they are not identical is the underlying $N$-body simulations which have a larger volume but lower resolution compared to the ones used in S19. We refer to Appendix B for a detailed discussion about resolution effects.

Throughout this paper, we will use the benchmark models defined above to quantify the range of realistic baryonic correction effects. They stand for an optimistic (model A), a best-guess (model B), and a pessimistic (model C) scenario regarding the baryonic effects on the cosmological density field. This characterisation is reasonable because the hydrostatic mass bias has been identified in S19 to be the largest systematic uncertainty of the analysis pipeline.

Before moving on to the result section, it is important to emphasise that, although consisting of an approximative method, the BC model has been shown to be in good agreement with full hydrodynamical simulations. As a test, S19 calibrated the BC model parameters using the measured gas fractions from hydrodynamical simulations at redshift zero. They then predicted the matter power spectrum based on the BC model (using the fitted parameters) and compared it with the power spectrum of the same hydrodynamical simulation. At redshift zero the agreement between the BC model and the simulations was better than two percent for all modes below $k=5$ h/Mpc. At redshift one and two the agreement degraded slightly but was still better than 3 and 5 percent, respectively. The decrease in accuracy towards higher redshift is not surprising since the model parameters have been calibrated against gas fractions at redshift zero. The comparison of S19 included the hydrodynamical simulations OWLS [60, 17], cosmo-OWLS with AGN8.0 and AGN8.5 [61, 19], Horizon-AGN [62, 22], and Illustris-TNG [21].

4 Results

With the BC model introduced in the previous section, it is straight-forward to include baryonic effects into the pipeline for weak lensing convergence maps. We simply apply the particle displacement algorithm to all $N$-body outputs at different redshifts, before building the particle light-cone and producing the angular maps. For each of the ten $N$-body runs, we perform the baryonic correction three times with model parameters corresponding to the three benchmark models summarised in Table 2.

As a first step we investigate the effects of baryons on the angular power spectrum measured from the convergence maps. Fig. 5 illustrates the ratio between the angular power spectrum with and without baryons ($C_{\rm bcm}$ and $C_{\rm dmo}$, respectively). Depending on the benchmark model, the baryonic effects lead to a suppression of up to 20 percent starting beyond $l\sim 200-400$ and dominating the statistical error at $l\sim 500-800$. These scales are relevant for current weak-lensing surveys such as KiDS, DES, or HSC. Recent results from the last of these have shown measurements of the angular power spectrum up to $l\sim 3000$ with error bars comparable to the amplitude of the baryon effects shown here [see Ref. 63].

The results of Fig. 5 are generally consistent with the findings of S19 [1]. Note however, that compared to S19, the amplitude of the baryon suppression is slightly smaller by about 1-2 percent below $l\sim 3000$ and $\sim 5$ percent below $l\sim 6000$. These differences can be attributed to the lower mass resolution of the simulations used here. See Appendix B for a more detailed analysis of resolution effects.

We now turn our attention towards the main goal of the paper which is the measurement of baryonic effects on the weak-lensing peaks. A convenient statistical measure to determine the significance of the baryon effects ($s_{\rm bcm}$) is the relative difference of the peak distributions divided by the noise, i.e.,

[TABLE]

A significance of $|s_{\rm bcm}|\ll 1$ means that the baryonic effects are subdominant compared to the statistical error and can be ignored. On the other hand, for the case of $|s_{\rm bcm}|>1$, the baryonic effects are expected to dominate the error budget and have therefore to be properly included in the analysis pipeline.

Fig. 6 illustrates the significance of the baryonic effects for the peak distributions of a DES-like survey configuration, where each panel refers to a different smoothing scale. Unsurprisingly, the smaller the smoothing scale the larger the significance of the baryonic feedback effects. A smoothing of $\theta_{k}=8$ arcmin or higher is sufficient to guarantee that baryonic effects stay subdominant over all scales of the DES configuration. This is not the case for smaller values of $\theta_{k}$, where the high peaks become significantly affected by baryons. For our most extreme case of $\theta_{k}=2$ arcmin, the baryons lead to a 1-2 $\sigma$ bias for scales above $\mathcal{S}/\mathcal{N}=3$. Smaller peaks below this threshold remain less affected by baryonic effects.

The results shown in Fig. 6 validate the approach of Kacprzak et al. [37], a cosmological parameter study based on peak counts of the DES science verification data. While this paper did not include baryonic effects, it used a sufficiently large smoothing of the map to fully mitigate all potential systematics from baryonic feedback. The situation could be different regarding the more recent peak count analysis from KiDS [40, 39] which applied a smoothing scale of $\sim 2$ arcmin and used both low and high signal-to-noise peaks. Our findings suggest that these results could indeed be affected by baryons at a significant level. Note, however, that more quantitative conclusions can only be obtained with a dedicated study that assumes a KiDS-like survey configuration in terms of mean galaxy number density, redshift distribution, and survey area.

In the future, baryonic feedback will become an important systematical effect for weak lensing surveys. In Fig. 7 we show the same information as in Fig. 6 but this time for a Euclid-like survey configuration. As expected, the strongly increased signal-to-noise leads to larger systematics from baryonic feedback. A smoothing of $\theta_{k}=16$ arcmin or above is required if baryonic effects are ignored in the peak modelling. However, such a large smoothing suppresses cosmological information and is therefore no viable option. At smaller smoothing scales below $\theta_{k}=8$ arcmin, baryonic feedback is significantly altering the signal. Regarding high peaks with $\mathcal{S}/\mathcal{N}>2$, we expect a bias of the order 2-4 $\sigma$ for $\theta_{k}=4$ arcmin and 5-8 $\sigma$ for $\theta_{k}=2$ arcmin. Furthermore, note that at these small smoothing scales the low peaks with $\mathcal{S}/\mathcal{N}<2$ also become affected at the 1-$\sigma$ level for $\theta_{k}=4$ arcmin and the 4-$\sigma$ level for $\theta_{k}=2$ arcmin.

While low $\mathcal{S}/\mathcal{N}$ peaks will also be influenced by baryons for future surveys like Euclid, the effect on high peaks remains clearly stronger. We speculate that the reason for this behaviour is that high peaks are caused by individual haloes whose profiles are directly affected by baryonic feedback. Low peaks, on the other hand, are caused by a combination of multiple line-of-sight haloes and shape noise [42, 43], where only the haloes are modified by baryonic feedback while the noise term stays the same.

Towards small smoothing scales, the relative peak counts of both DES and Euclid type surveys show a characteristic wave shape (see top panels of Fig. 6 and 7) with $n_{\rm bcm}>n_{\rm dmo}$ at high and $n_{\rm bcm}<n_{\rm dmo}$ at low values of $\mathcal{S}/\mathcal{N}$. Since the total number of peaks changes by less than one percent between the baryonic corrected and the dark matter-only sample, we conclude that baryonic effects cause peaks to simply shift towards lower $\mathcal{S}/\mathcal{N}$ values without destroying them.

Based on the analysis shown above we conclude that for current, stage-III weak-lensing surveys, the baryonic suppression effects are at most comparable to the statistical error bars. We therefore expect that ignoring baryonic feedback does not strongly bias cosmological parameter estimates. For future, stage-IV weak-lensing surveys, on the other hand, it will be essential to properly include baryon effects in the prediction pipeline. Most notably, this is true for small smoothing scales below $\theta_{k}\sim 8$ arcmin which contain important information about cosmology.

5 Conclusion

In this paper we have investigated the effects of baryonic feedback on peak number statistics of the weak-lensing convergence field. The convergence maps were constructed via light-cones of $N$-body simulations which have been integrated along the line-of-sight using appropriate lensing kernels and simple noise configurations. The effects of baryons were implemented using the baryon correction model of S19 which can be directly applied to N-body simulations and has previously been shown to provide an accurate prescription of the baryon suppression effects at scales relevant for weak lensing [1]. Following S19 we defined three benchmark baryonic models motivated by X-ray observations of gas fractions which provide a realistic range of scales for the baryonic suppression effects.

In a first step we validated our pipeline by measuring the angular power spectrum using fifty full-sky convergence maps based on light-cones of replicated $N$-body outputs. A comparison to results from the literature based on the Limber approximation shows good agreement regarding both the baryon suppression effect and the absolute power spectrum of the dark-matter-only case.

As a further step we investigated the convergence peak count, assuming two distinct configurations for the survey footprint, the mean galaxy density, and the noise implementation. The first configuration is motivated by the characteristics of DES, while the second is based on the expected setup of Euclid. This allowed us to quantify the baryonic feedback effects of both current stage-III and future stage-IV weak-lensing experiments.

Not surprisingly, we found that the strength of the baryon suppression on the peak statistics strongly depends on the scale of the Gaussian smoothing applied to the convergences maps. For the DES-like configuration, a smoothing scale of $\theta_{k}\sim 8$ arcmin or larger is sufficient to wash-out the baryonic effects so that the suppression signal becomes smaller than the statistical error from limited number counts. For the Euclid-like configuration, on the other hand, a smoothing of at least $\theta_{k}\sim 16$ arcmin is required to reduce the baryon suppression signal enough so that it falls below the statistical error.

Smaller smoothing scales than the ones mentioned above require a proper modelling of baryon suppression effects for the weak-lensing convergence map. Assuming a DES-like survey configuration and a Gaussian smoothing of $\theta_{k}=2$ ($\theta_{k}=4$) arcmin, the baryon effects become larger than the statistical error for high peaks with $\mathcal{S}/\mathcal{N}>3$ ($\mathcal{S}/\mathcal{N}>4$). Low peaks are also modified by baryon feedback, but this effect stays much smaller than the sample variance error due to the limited survey footprint.

The consequences of baryon feedback are more serious for the case of a Euclid-like survey setup. Assuming a smoothing scale of $\theta_{k}=4$ arcmin or below, the baryon effects exceed the expected statistical errors over the entire range of $\mathcal{S}/\mathcal{N}$. This means that for Euclid, not only high but also low peaks are significantly affected by baryons and a proper modelling of baryon feedback effects becomes indispensable.

Acknowledgments

We acknowledge support from the Swiss National Science Foundation. This includes the grants PZ00P2_161363 and 200021_169130.

Appendix A Replication and Randomization Effects

The goal of this appendix is to show that the results or the paper are unaffected by the box-replication procedure that we apply to the $N$-body outputs in order to construct full-sky light-cones and convergence maps. We test the replication procedure by running simulations with four different box-sizes leading to different replication counts. The angular power spectra of the simulations with box-lengths $L=128,\,256,\,512,\,1024$ Mpc/h with and without randomization are illustrated in Fig. 8.

The top panel of Fig. 8 shows the angular power spectra measured from the convergence maps that were obtained via the light-cone construction discussed in Sec. 2.2. Different lines refer to different box-sizes ($L$) and therefore different numbers of box replica (since all light-cones cover the same range in redshift of $z=0.1-1.5$). The plot shows that the setup used in the main analysis (i.e. the blue line with box-size $L=512$ Mpc/h and $12^{3}$ replica) is well converged to the level of a few percent.

In the bottom panel of Fig. 8 we compare the randomised to the non-randomised version of the light-cone production, again for the standard case of $L=512$ Mpc/h and $12^{3}$ replica. The angular power spectra of both versions agree very well except for a small ($\sim 2$ percent) difference at the larges physical scales. Furthermore, both power spectra are in good agreement with the one obtained via the Limber approximation and the three dimensional matter power spectra from the $N$-body outputs.

In a second step we perform the same type of test for the peaks of the convergence maps, using the Euclid-like survey configuration and a $\theta_{k}=$ 2 arcmin smoothing scale. Fig. 9 shows the resulting peak differences for the different box-sizes (left) and the randomised/non-randomised case (right). Comparing this result to Fig. 3 we see see the effect of the replication or randomization on the peak counts is negligible compared to the cosmological signal.

From testing the angular power spectrum and the peak statistics we conclude that errors due to box-size, numbers of replications, or the randomization method remain subdominant over all relevant ranges of scales. Furthermore, since we mainly focus on the relative effects of the baryonic versus dark-matter-only cases, all potential sources of errors investigated above are expected to cancel out.

Appendix B Resolution Effects

As mentioned in Section 2.1, we only modify haloes with more than 100 particles, i.e. with a halo mass above $M_{\text{halo}}=9\times 10^{12}$ M⊙/h. This is not enough to get fully converged power spectra up to the scales relevant for future weak-lensing surveys [see S15 30]. In Fig. 10 we show that the angular power spectrum from our simulations (solid lines) are not identical to the converged results from ST19. There is a difference of a few percent above $l=3000$ which can be attributed to the fact that baryon effects in haloes below our resolution limit contribute to the matter power spectrum at small scales. Note that the percent level difference at lower at $l<3000$ is not a resolution effect but is a result of the fact that the dash-dotted lines were derived using the Limber approximation (see S19). For comparison, we also plot the case of a lower mass resolution (of $N=256$ with the same box-size $L=512$ Mpc/h), which shows a significantly weaker baryon suppression (dashed lines). Although our results are not converged at the percent level for $l>3000$, we estimate the results to be good enough for the purpose of this paper.

Finally we investigate potential resolution effects on the peak statistics. Fig. 11 shows the absolute and relative effect of all three benchmark models with both simulations using the Euclid-like survey and a $\theta_{k}=2$ arcmin smoothing scale. The solid lines correspond to the standard resolution ($N=512^{3}$) used in the main analysis, while the dashed lines show results of a reduced resolution ($N=256^{3}$). The figure shows that while the high peaks are fully converged this is not necessarily the case for the lower peaks. Unfortunately, we do not have an even higher resolution to establish if there is a remaining discrepancy.