# Comparing weak lensing peak counts in baryonic correction models to
hydrodynamical simulations

Max E. Lee,1 Tianhuan Lu,1 Zoltán Haiman,1 Jia Liu,2 Ken Osato3,4
1Department of Astronomy, Columbia University, New York, NY 10027, USA
2Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583,
Japan
3Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University,
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan
4LPENS, Département de Physique, École Normale Supérieure, Université PSL,
CNRS, Sorbonne Université, Université de Paris,
24 rue Lhomond, 75005 Paris, France
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ
Report number:YITP-21-162)

###### Abstract

Next-generation weak lensing (WL) surveys, such as by the Vera Rubin
Observatory’s LSST, the Roman Space Telescope, and the Euclid space mission,
will supply vast amounts of data probing small, highly nonlinear scales.
Extracting information from these scales requires higher-order statistics and
the controlling of related systematics such as baryonic effects. To account
for baryonic effects in cosmological analyses at reduced computational cost,
semi-analytic baryonic correction models (BCMs) have been proposed. Here, we
study the accuracy of BCMs for WL peak counts, a well studied, simple, and
effective higher-order statistic. We compare WL peak counts generated from the
full hydrodynamical simulation IllustrisTNG and a baryon-corrected version of
the corresponding dark matter-only simulation IllustrisTNG-Dark. We apply
galaxy shape noise expected at the depths reached by DES, KiDS, HSC, LSST,
Roman, and Euclid. We find that peak counts in BCMs are (i) accurate at the
percent level for peaks with ${\rm S/N}<4$, (ii) statistically
indistinguishable from IllustrisTNG in most current and ongoing surveys, but
(iii) insufficient for deep future surveys covering the largest solid angles,
such as LSST and Euclid. We find that BCMs match individual peaks accurately,
but underpredict the amplitude of the highest peaks. We conclude that existing
BCMs are a viable substitute for full hydrodynamical simulations in
cosmological parameter estimation from beyond-Gaussian statistics for ongoing
and future surveys with modest solid angles. For the largest surveys, BCMs
need to be refined to provide a more accurate match, especially to the highest
peaks.

††pubyear: 2021††pagerange: Comparing weak lensing peak counts in baryonic
correction models to hydrodynamical simulations–Comparing weak lensing peak
counts in baryonic correction models to hydrodynamical simulations

##  1 Introduction

Weak gravitational lensing (WL) measures the correlated distortions in the
shapes of distant galaxies due to light deflections induced by gravitational
potentials along the line of sight direction (see, e.g. Bartelmann & Schneider
2001; Hoekstra & Jain 2008; Kilbinger 2015 for reviews). Because of its
sensitivity to the clustering of matter, WL is a promising cosmological probe.
WL has already been shown to be a powerful tool for constraining parameters of
the standard $\Lambda$CDM model, in particular the total matter density
parameter, $\Omega_{\mathrm{m}}$, the amplitude of matter fluctuations
$\sigma_{8}$, and their commonly used combination
$S_{8}=\sigma_{8}(\Omega_{\mathrm{m}}/0.3)^{0.5}$ (Hikage et al., 2019; Hamana
et al., 2020; DES Collaboration et al., 2021).

The most common approaches to WL analysis rely on quantities related to the
two-point correlation function of the shear field, or its Fourier transform,
the power spectrum (Heymans et al., 2013; Kitching et al., 2014; Hildebrandt
et al., 2017; Hikage et al., 2019). However, it has been found that these
statistics do not capture all the information from the non-linear and non-
Gaussian regime of the WL field. This has led to an array of work attempting
to extract this additional information, including non-Gaussian $N$-point
statistics such as the 3-point function and its Fourier transform the angular
bispectrum (Fu et al., 2014; Coulton et al., 2019; Halder et al., 2021; Halder
& Barreira, 2022), Minkowski functionals (Maturi et al., 2010; Munshi et al.,
2011; Kratochvil et al., 2012; Petri et al., 2013; Marques et al., 2019; Gatti
et al., 2021), convergence peak and minima statistics (Jain & Van Waerbeke,
2000; Yang et al., 2011; Liu et al., 2015; Kacprzak et al., 2016; Li et al.,
2019; Martinet et al., 2021; Zürcher et al., 2021; Harnois-Déraps et al.,
2021), as well as machine learning methods (Gupta et al., 2018; Fluri et al.,
2019; Ribli et al., 2019; Matilla et al., 2020; Lu et al., 2021) and related
techniques (Cheng et al., 2020).

Peak statistics have been shown to be a simple and effective non-Gaussian
statistic of the WL field for constraining parameters and testing cosmological
models. They represent pixels in the WL field that are larger in value when
compared to their eight nearest neighbors (Jain & Van Waerbeke, 2000; Dietrich
& Hartlap, 2010; Kratochvil et al., 2010; Liu et al., 2015). These peaks are
correlated with massive individual dark matter haloes, or with numerous
intermediate-mass haloes, making them sensitive to the structural evolution
and history of the universe (Yang et al., 2011; Kilbinger, 2015; Osato et al.,
2015; Liu & Haiman, 2016; Osato et al., 2021). It has recently been found in
Sabyr et al. (2021) that beyond haloes, peaks are also sensitive to inter-halo
matter.

With future surveys covering $18,000\,\mathrm{deg}^{2}$ and galaxy number
densities of $50\,\mathrm{arcmin}^{-2}$ (LSST Science Collaboration et al.,
2009; Laureijs et al., 2011), systematic effects, such as intrinsic alignments
of galaxy shapes (Schneider & Bridle, 2010; Sifón et al., 2015), and baryonic
physics (Jing et al., 2006; Chisari et al., 2019) will dominate the errors on
small scales. As beyond-Gaussian statistics, such as peak counts, are
particularly sensitive to small scales, it is critical to understand such
systematics for these statistics.

In this work, we focus our attention on the influence of baryonic physics on
weak lensing peak statistics. Baryonic effects can be incorporated by
hydrodynamical simulations that follow subgrid prescriptions calibrated to
observations (Schaye et al., 2010; McCarthy et al., 2017; Springel, 2010;
Nelson et al., 2019), by simpler semi-analytical baryonic correction models
(Schneider & Bridle, 2010; Schneider et al., 2019; Aricò et al., 2020), or
through machine learning methods (Villaescusa-Navarro et al., 2020; Dai &
Seljak, 2021; Lu & Haiman, 2021). While hydrodynamical simulations which have
been appropriately calibrated are ideal testing grounds for baryonic effects
and cosmological analyses, they require expensive simulation suites.

The semi-analytical baryonic correction models (BCMs) are attractive for
cosmological analyses, as they correct computationally cheaper dark matter-
only simulations without the need to explicitly simulate expensive baryonic
processes. In Schneider & Bridle (2010), Schneider et al. (2019), and Aricò et
al. (2020), BCMs with increasingly complex semi-analytical prescriptions for
baryonic effects, such as star formation, active galactic nuclei (AGN)
feedback, and gravitationally bound and unbound gas have been proposed and
compared with full hydrodynamical simulations. In particular, Aricò et al.
(2020) have recently shown that their BCM is able to reproduce the power
spectra of the 3D matter field to within $1\%$ accuracy in a broad range of
state-of-the-art hydrodynamical simulations.

To date, BCMs have been calibrated and tested almost exclusively on two-point
statistics. The exception is the recent study by Aricò et al. (2021a), who
examined a joint fit of their BCM to the 3D matter power spectrum and
bispectrum, and found a $<3\%$ agreement with hydrodynamical simulations. In
this work, we perform a similar comparison, focusing on WL peak counts—a
statistic that contains information from all orders. We seek to answer the
following questions:

  1. 1.

How well do peak counts from BCM match hydrodynamical simulations?

  2. 2.

What is the physical cause of any apparent deviations?

  3. 3.

Can peak counts from BCMs be used as a replacement for hydrodynamical
simulations in current and future survey analyses?

  4. 4.

Can the agreement between BCM and hydrodynamical peak counts be improved?

This paper is organised as follows. In § 2 we describe the simulations used,
our weak lensing field computations, and the application of the Aricò et al.
(2020) baryonic correction model. We then present our comparisons between peak
counts derived from hydrodynamical and BCM weak lensing maps in § 3. In § 4 we
consider causes for the deviations between BCM and hydrodynamical peak counts,
and explore the influence of galaxy shape noise on our results. This allows a
quantification of how much the BCM will deviate from a hydrodynamical
simulation if used in current and future analyses. We further propose and test
a simple solution—excluding the highest peaks—to help mitigate these
discrepancies. We summarise our main results and the implications of this work
in § 5.

##  2 Methods

In this section we discuss the hydrodynamical and dark matter-only simulations
used to generate WL maps. We then review the BCM of Aricò et al. (2020) and
its application to dark matter-only simulations. We end by describing the
process of ray-tracing and the extraction of peak statistics.

###  2.1 $N$-body and Hydro Simulations

In this work we use the highest-resolution run of the magneto-hydrodynamical
and dark matter-only cosmological simulations from IllustrisTNG-300, and
IllustrisTNG-300-Dark (Pillepich et al., 2017; Springel et al., 2018; Nelson
et al., 2019; Naiman et al., 2018; Marinacci et al., 2018, hereafter TNG and
TNG-Dark, respectively). TNG uses the adaptive moving-mesh code AREPO
(Springel, 2010) on a $(205\,h^{-1}\mathrm{Mpc})^{3}$ simulation box with
$2500^{3}$ dark matter particles and $2500^{3}$ initial gas cells. While both
TNG and TNG-Dark contain the same initial random seed, the subgrid model of
TNG incorporates baryonic physics, such as AGN feedback, star-formation and
evolution, as well as radiative cooling (Pillepich et al., 2017).

For TNG, the mass of a dark matter particle is $2.98\times
10^{7}\,h^{-1}\,M_{\odot}$ and an initial gas cell mass is $1.1\times
10^{7}\,h^{-1}\,M_{\odot}$, while TNG-Dark has a dark matter particle mass of
$4.72\times 10^{7}\,h^{-1}\,M_{\odot}$. Both simulations start at a redshift
of $z=127$ and use $100$ snapshots to arrive at $z=0$, though in our work we
focus our attention on snapshots between $0\leq z\leq 2.5$, covering the range
of redshifts in WL surveys.

###  2.2 Baryonic Correction Model

Figure 1: The density profiles of the four baryonic correction model (BCM)
components and the total correction ($\times 2$ for visualization purposes)
applied to a halo with mass $M_{200}\simeq 10^{14}\,h^{-1}M_{\odot}$. The
lines here represent theoretical predictions for the BCM components, while the
points show the measured values sampled from the halo after applying the BCM.
This validates our implementation of the BCM outlined in A20 (through the
close match with Fig. 1 of A20).

We use the four-parameter semi-analytical baryonic correction model outlined
in Aricò et al. (2020, hereafter A20). A20’s corrections account for cooling
of gas in halo centres, the ejection of gas from the haloes, hydrostatic
equilibrium of hot bound gas, star formation, and adiabatic relaxation of dark
matter due to baryons. Here we briefly summarise the key components and free
parameters outlined in A20 as well as in Lu & Haiman (2021); interested
readers are referred to these papers for more details.

In the BCM, density profiles of haloes are modified such that,

$\rho_{\mathrm{DMO}}(r)\to\rho_{\mathrm{BCM}}(r)=\rho_{\mathrm{CG}}(r)+\rho_{\mathrm{HBG}}(r)+\rho_{\mathrm{EG}}(r)+\rho_{\mathrm{RDM}}(r).$
(1)

We begin with a dark matter halo density profile $\rho_{\mathrm{DMO}}(r)$,
which is then modified to be composed of a central galaxy component,
$\rho_{\mathrm{CG}}(r)$, hot bound gas $\rho_{\mathrm{HBG}}(r)$, ejected gas,
$\rho_{\mathrm{EG}}(r)$, and relaxed dark matter $\rho_{\mathrm{RDM}}(r)$,
each with its own density profile. Once these density profiles are known, the
enclosed-mass profile of the halo can be found by integrating Eq. (1) to find
$M_{\mathrm{BCM}}(r)$. We then invert this function to find
$r_{\mathrm{BCM}}(M)$, and the corresponding correction of each simulated
particle’s radial position from the halo’s centre of mass,

$\Delta r(M)=r_{\mathrm{BCM}}(M)-r(M).$ (2)

For each halo with a virial mass greater than $10^{12}\,h^{-1}\,M_{\odot}$ we
extract all of the particles from its main sub-halo. This cut is reasonable
because, as shown in A20, haloes with $M\leq 10^{13}\,h^{-1}\,M_{\odot}$ have
an impact of $\leq 2\%$ on the total baryonic suppression at scales of
$k\approx 5\,h\,\mathrm{Mpc}^{-1}$. For each sub-halo, we then fit the dark
matter particles to a truncated NFW profile defined as (Navarro et al., 1996)

$\begin{split}\rho_{\mathrm{NFW}}=\begin{cases}\frac{\rho_{0}}{(r/r_{\mathrm{s}})(1+r/r_{\mathrm{s}})^{2}}\quad&{r\leq
r_{200}}\\\ 0&{r>r_{200}}\end{cases},\end{split}$ (3)

where $\rho_{0}$ is the characteristic density of the halo,
$r_{\mathrm{s}}=r_{200}/c$, $c$ is the concentration parameter, and $r_{200}$
is the radius at which the average density is $200$ times the critical density
of the universe. We define $r_{200}$ as the halo radius and $M_{200}$ as the
halo mass. The initial density profile $\rho_{\mathrm{DMO}}$ is then set to
the NFW profile in Eq. (3).

The central galaxy component follows a power-law with an exponential cutoff as
in Kravtsov et al. (2018). The slope and cutoff of the density profile is
purely determined by the halo mass and radius, but the mass fraction
$M_{\mathrm{galaxy}}/M_{\mathrm{halo}}$ is found from sub-halo abundance
matching and uses the best-fit parameters found by Behroozi et al. (2013). To
provide further flexibility, $M_{1}$, which is the mass of a halo with a
galaxy mass fraction of $M_{\mathrm{galaxy}}/M_{\mathrm{halo}}=0.023$ at
$z=0$, is left as a free parameter.

We model hot bound gas in the halo as a fluid in hydrostatic equilibrium up to
$\sqrt{5}r_{200}$ at which point it follows the slope of the truncated NFW
profile of Eq. (3). The fraction of bound gas is determined by the power law
$(M_{\mathrm{h}}/M_{\mathrm{c}})^{\beta}$ where $M_{\mathrm{h}}$ is the halo
mass, $M_{\mathrm{c}}$ is the mass at which half of the gas is retained inside
of the halo, and $\beta$ is the slope of the gas fraction. Both
$M_{\mathrm{c}}$ and $\beta$ are free parameters of the BCM.

Ejected gas assumes a Maxwell–Boltzmann velocity distribution with isotropic
ejection. Particles are allowed to reach a maximum radius of
$r_{\mathrm{ej}}=\eta\cdot 0.75r_{\mathrm{esc}}$ where $r_{\mathrm{esc}}$ is
determined by considering the distance a gas particle with constant halo
escape velocity, $v_{\mathrm{esc}}$, travels in half of the Hubble time
(Schneider & Bridle, 2010),

$r_{\mathrm{esc}}=0.5t_{\mathrm{Hubble}}\sqrt{\dfrac{8}{3}G\Delta_{200}\rho_{c}}\,r_{200}\sim
0.5\sqrt{\Delta_{200}}\,r_{200}.$ (4)

In the BCM of A20, $\eta$ is the fourth and final free parameter.

After considering the addition of baryonic effects, we introduce an adiabatic
relaxation to the dark matter particles. This allows for positional adjustment
given the gravitational interactions with the baryons. Thus regions of deep
potentials, near the centre of haloes where we place the central galaxies, the
dark matter particles contract. Regions near the edges of haloes, where
ejected gas carries particles away, cause an expansion of the dark matter
particles.

We follow the original implementation of the BCM as in A20, in which we move
each simulated DM particle in a given halo according to the BCM’s semi-
analytical prescription. Note that more efficient methods which rely on
surface density projections (Lu & Haiman, 2021), and emulators (Aricò et al.,
2021b) exist. For the purpose of this work, we find that the computational
cost of individual particle displacement is not prohibitive, but for future
work, we plan to implement these inexpensive simplifications.

We used $8$ nodes111Intel Xeon Phi 7250 KNL nodes on TACC’s Stampede2, see
https://portal.tacc.utexas.edu/user-guides/stampede2 for further node details
and $\sim 1000$ core hours to correct each simulation snapshot from TNG-Dark
on TACC’s Stampede2. Using twenty TNG-Dark snapshots (see § 2.3), this
corresponds to a total of $\sim 20,000$ core hours. We have fully parallelised
our BCM implementation using mpi4py222https://github.com/mpi4py/mpi4py(Dalcin
& Fang, 2021) which generates corrected maps far more efficiently. For
simplicity, in this study we adopt the parameters already identified to best
fit the 3D matter power spectrum in A20, as listed in Table 1. In principle,
optimal parameters could instead be found by fitting the WL peak counts
directly. This could potentially improve the BCM’s accuracy, but would require
substantially more computational resources and work. We save this
investigation for a future study, and simply note here that our results below
are conservative.

As a test of our implementation, we compare the density profiles of each BCM
component in Fig. 1. This plot gives a good intuition as to what the BCM is
doing in general. The central galaxy component causes a sharp density increase
at small radius which exponentially falls off, the warm gas then dominates the
baryonic effects up to the halo edge at $r_{200}$. Near the outskirts of the
halo, the ejected gas becomes dominant, causing particles originally inside
the $r_{200}$ boundary to become expelled. Our Fig. 1 is in excellent
agreement with Figure 1 in A20, validating our BCM implementation in TNG-Dark.

Parameter | Value | Meaning
---|---|---
$M_{1}\,[10^{11}\mathrm{M}_{\odot}h^{-1}]$ | 0.22 |  $M_{200}$ of halo with galaxy-mass fraction
|  | of $0.023$ at $z=0$
$M_{\mathrm{c}}\,[10^{14}\mathrm{M}_{\odot}h^{-1}]$ | 0.23 |  $M_{200}$ of halo with half gas bound
$\beta$ | 4.09 | slope of gas fraction vs. halo mass
$\eta$ | 0.14 | Maximum distance of ejected gas in units
|  | of $r_{\mathrm{esc}}$
Table 1: The parameter values used for the baryonic correction model. These
are adopted from the fit in A20 to the 3D matter power spectrum, which have
been shown to yield $\leq 1\%$ errors.

###  2.3 Ray-Tracing

Figure 2: We show a typical $\kappa$TNG map at $z=1.5$ (_left panel_), as well
as the difference between a $\kappa$TNG map and $\kappa$Dark (_middle panel_)
and $\kappa$BCM (_right panel)_ that were generated with the same initial
seeds. In the right two panels, we see a blue-red dipole pattern frequently
occurring. We interpret this as being due to baryonic effects outside of
haloes in TNG causing a shift of the centres of large haloes not captured by
$\kappa$Dark and $\kappa$BCM (see § 4 for discussion).

Here we briefly summarise the process outlined and used in Osato et al.
(2021), which generated weak lensing convergence maps ($\kappa$ maps) from the
TNG and TNG-Dark simulations. We use twenty TNG snapshots in intervals of
approximately equally spaced comoving distance, each with a comoving volume of
$(205\,\text{Mpc}\,h^{-1})^{3}$. This allows the generation of a light cone
with a $5\times 5\,\mathrm{deg}^{2}$ coverage333There is a tiling of snapshots
used for farther redshifts to cover the full extent of the light cone. See
Osato et al. (2021) for details.. We place source galaxies at four redshifts
corresponding to $z=0.4,0.7,1.0,1.5$ and model lensing convergence using the
multiple lens plane algorithm of Jain et al. (2000).

In this algorithm, each snapshot is split into two projections spanning a
comoving distance of $102.5\,\mathrm{Mpc}\,h^{-1}$. The projection is then
pixelated onto a $4096\times 4096$ mesh, on which the 2D lensing potential and
deflection angles are computed. The change in lensed position ($\beta$), due
to the deflection angle at each plane, with respect to its initial position
($\theta$) then represents the standard distortion matrix in weak lensing,

$\begin{split}A_{ij}(\vec{\theta},\chi)&=\dfrac{\partial\beta_{i}(\vec{\theta},\chi)}{\partial\theta}\\\
&=\begin{pmatrix}1-\kappa-\gamma_{1}&-\gamma_{2}+\omega\\\
-\gamma_{2}-\omega&1-\kappa+\gamma_{1}\end{pmatrix}\end{split}$ (5)

where $\kappa$ is the weak lensing convergence, $\gamma_{1}$ and $\gamma_{2}$
are the shear components, and $\omega$ represents rotation. We employ the ray-
tracing scheme of Hilbert et al. (2009) to trace the light rays back from
$z=0$ to the source redshift through each lensing potential along the light
path. This method has been shown in many past studies to be memory efficient
and accurate for WL simulations.

The resultant suite, dubbed $\kappa$TNG, contains $20,000$ total pseudo-
independent convergence maps, each with $1024\times 1024$ pixels, half of
which are from ray-tracing the TNG and the other half from TNG-Dark
simulations. We emphasise that these maps are pseudo-independent because each
map was generated with the same twenty snapshots from TNG (TNG-Dark, BCM),
though the particles in each snapshot were randomly translated, rotated by
$0$, $90$, $180$, or $270$ degrees, and flipped along each of the three axes
prior to the generation of lensing potential planes. In $\kappa$TNG this
process was applied $100$ times, followed by a repetition of this process to
each lens-plane another $100$ times to generate $10000$ pseudo-independent
realizations. This process of translation and rotation to generate independent
maps has been studied in Petri et al. (2016) and applied in Osato et al.
(2021) providing statistically independent power spectra and peak counts.

Because of the computational cost of performing ray-tracing, we do not compute
$10,000$ BCM $\kappa$ maps ($\kappa$BCM maps), and instead generate $2,000$
pseudo-independent convergence maps for the purpose of suppressing sample
variance. For each of $\kappa$BCM map, we find the concordant map (same random
seeds) from $\kappa$TNG and $\kappa$Dark giving us a total of $6,000$ maps.
Prior to analysis, all maps are smoothed with a Gaussian filter defined as

$W(\theta)=\frac{1}{\pi\theta_{\mathrm{G}}^{2}}\exp\left(-\frac{\theta^{2}}{\theta_{\mathrm{G}}^{2}}\right),$
(6)

where $\theta_{\mathrm{G}}=1\,\text{arcmin}$. This filter removes spurious
small-scale noise, and was used in Osato et al. (2021). In the cases described
in § 4, we apply noise prior to smoothing. We illustrate in Fig. 2 a
$\kappa$TNG map, and the differences $\kappa\mathrm{TNG}-\kappa\mathrm{Dark}$,
and $\kappa\mathrm{TNG}-\kappa\mathrm{BCM}$ but save a discussion of this
figure for § 4 below.

###  2.4 Peak Statistics

We are interested in exploring the non-Gaussian statistic of peak counts from
$\kappa$ maps. Peak counts are defined as the number of pixels with a given
$\kappa$ value, that are greater than their surrounding eight neighbors. We
explore peak counts from noiseless maps, and maps with Gaussian shape noise
added to each pixel with a mean $\mu=0$ and variance

$\sigma^{2}=\frac{\sigma_{e}^{2}}{2n_{\mathrm{gal}}A_{\mathrm{pix}}}.$ (7)

Here $\sigma_{e}=0.4$ is the assumed mean intrinsic ellipticity of galaxies,
$n_{\mathrm{gal}}$ is the surface number density of lensed galaxies, and
$A_{\mathrm{pix}}$ is the solid angle of a pixel. We consider noise levels
corresponding to current surveys with observed $n_{\mathrm{gal}}=10$ or
$20\,\mathrm{arcmin}^{-2}$, and future surveys with $n_{\mathrm{gal}}=30$ and
$50\,\mathrm{arcmin}^{-2}$.

We denote the average peak counts over all realizations ($N=2000$) of a method
as $\langle\mathbf{n}\rangle$. Throughout this work, we represent our peak
count histograms as a function of their significance in convergence maps, or
${\rm S/N}$. Here the “noise” $N$ is taken as the root-mean-square $\kappa$
value ($\kappa_{\mathrm{rms}}$) of convergence maps, both noisy and noiseless,
averaged over all $2000$ realizations. We show in Table 3 the
$\kappa_{\mathrm{rms}}$ values from $\kappa$BCM and $\kappa$Dark in noiseless
convergence maps for a range of redshifts. Choosing any of the columns in
Table 3 is sufficient and only effects the $x$-axis scale by $<$ a percent
level. This justifies the use of only $\kappa$TNG’s $\kappa_{\rm rms}$ values
throughout this work. $\kappa_{\rm rms}$ changes with noise level and
redshift, and for each we compute the associated average $\kappa_{\rm rms}$
from $\kappa$TNG to find the S/N. We create 18 equally spaced bins from
$-1.5\leq{\rm S/N}\leq 7.5$ for each peak count histogram.

The significance of difference in peak counts between two methods, such as
between the $\kappa$BCM and $\kappa$TNG, depends on two quantities: (i) the
difference in the expectation value of the peak counts histogram
$\Delta\langle\mathbf{n}\rangle$, and (ii) the peak counts covariances
$\mathbf{C}$. The latter can be defined in either group of maps by

$C_{ij}=\frac{1}{N-1}\sum_{k=1}^{N}\left(n_{i}^{(k)}-\langle
n_{i}\rangle\right)\left(n_{j}^{(k)}-\langle n_{j}\rangle\right),$ (8)

where $N=2000$ denotes the number of realizations. For simplicity, we take the
peak counts in the $\kappa$TNG maps to calculate the covariances here. We have
found that evaluating the covariance matrix in the $\kappa$TNG or the
$\kappa$BCM suites yields nearly identical $\chi^{2}$ values, and makes no
difference to our results. Since the peak counts histogram generally follows a
multivariate normal distribution (Gupta et al., 2018), the quantity

$\chi^{2}\equiv{\left(\Delta\langle\mathbf{n}\rangle\right)}^{T}\,\hat{\mathbf{C}}^{-1}\,\Delta\langle\mathbf{n}\rangle,$
(9)

is a $\chi^{2}$ test (degrees of freedom $d=18$) for the difference between
two peak counts histograms drawn from both methods, where

$\hat{\mathbf{C}}^{-1}=\frac{N-d-2}{N-1}\mathbf{C}^{-1},$ (10)

and prefactor $(N-d-2)/(N-1)$ is used to debias the precision matrix
estimation (Hartlap et al., 2007). In Fig. 3, we validate that our definition
of $\chi^{2}$ in Equation (9) indeed follows a $\chi^{2}$ distribution with
$d=18$ when peak counts are sampled from $\kappa$TNG maps, i.e.
$\Delta\langle\mathbf{n}\rangle$ is replaced by
$\left(\mathbf{n}-\langle\mathbf{n}\rangle\right)$. This allows us to directly
interpret the $\chi^{2}$ as the confidence at which one can reject the null
hypothesis that the average value of the $\kappa$BCM peak histogram was drawn
from the $\kappa$TNG distribution. In particular, the standard 1, 2, and
3$\sigma$, or 68.27, 95.45 and 99.73% confidence levels correspond to
probabilities of finding $\chi^{2}\leq 20.28$, 29.24 and 39.17, respectively,
by random chance, for a $\chi^{2}$ distribution with 18 degrees of freedom.

The above calculation only applies to comparing individual peak counts
histograms of $5\times 5\,\mathrm{deg}^{2}$ maps from different methods, which
corresponds to a $25\,\mathrm{deg}^{2}$ survey. To derive the significance of
difference for larger surveys, we shall generalize the calculation. We note
that the peak counts histogram of a survey with a area $A$ is equal to the
average of $A/(25\,\mathrm{deg}^{2})$ draws of $25\,\mathrm{deg}^{2}$
histograms, which does not change the expectation value
$\langle\mathbf{n}\rangle$ but scales the covariance matrix by

$\mathbf{C}(A)=\left(\frac{A}{25\,\text{deg}^{2}}\right)^{-1}\mathbf{C}(25\,\mathrm{deg}^{2}),$
(11)

where $\mathbf{C}(25\,\mathrm{deg}^{2})$ follows Equation (8). $\chi^{2}$ of
this larger survey is modified accordingly:

$\chi^{2}(A)=\frac{A}{25\,\text{deg}^{2}}\,\chi^{2}(25\,\mathrm{deg}^{2}),$
(12)

where $\chi^{2}(25\,\mathrm{deg}^{2})$ follows Equation (9). Intuitively, when
we switch from a small survey to a larger one, the peak counts will be of
similar values but the uncertainties of the counts will become smaller;
therefore, with the same absolute difference between $\kappa$BCM and
$\kappa$TNG, our ability of distinguishing one method from the other is
improved ($\chi^{2}$ is proportional to $A$). Hereafter, we refer to
$\chi^{2}(A)$ as simply $\chi^{2}$ when it is calculated for each specific
survey.

Figure 3: The distribution of the $\chi^{2}$ statistic defined in Eq. (9),
evaluated within the $\kappa$TNG simulation suite. This statistic follows a
$\chi^{2}$ distribution for 18 degrees of freedom (using 18 peak-count bins),
as expected if the peak counts followed a multivariate Gaussian distribution.
This allows us use $\chi^{2}$ to quantify the significance of the difference
between the $\kappa$BCM and $\kappa$TNG peak-count histograms.

###  2.5 Summary

We summarise our process below.

  1. 1.

Extract primary sub-halo from all haloes with $M_{200}\geq
10^{12}\,h^{-1}\mathrm{M}_{\odot}$ from TNG-Dark, for $20$ snapshots between
$z=0$ and $z\approx 2.5$, using halo catalogs released with TNG-Dark

  2. 2.

Fit NFW profiles (Eq. (3)) to each of the haloes to find
$\rho_{\mathrm{DMO}}(r)$

  3. 3.

Compute BCM density profiles
$\rho_{\mathrm{BG}},\rho_{\mathrm{EG}},\rho_{\mathrm{CG}},\rho_{\mathrm{RDM}}$

  4. 4.

Find the new density profile (Eq. (1)) and apply radial shifts for each
particle (Eq. (2))

  5. 5.

Compute 2000 $\kappa$ maps following the prescription of Osato et al. (2021)
using BCM corrected 3D matter fields

  6. 6.

Add galaxy shape noise following Eq. (7)

  7. 7.

Smooth all maps with a $\theta_{\mathrm{G}}=1\,\mathrm{arcmin}$ Gaussian
filter (Eq. (6))

  8. 8.

Compute the peak counts as described in § 2.4 for $\kappa$TNG, $\kappa$Dark,
and $\kappa$BCM

  9. 9.

Calculate covariance matrices (Eqs. (8) and (11)) and finally the $\chi^{2}$
(Eqs. (9) and (12)) between the $\kappa$BCM and $\kappa$TNG peak count
histograms.

##  3 Results

Figure 4: The suppression in the weak lensing convergence power spectrum defined as $S(\ell)=\langle P_{\mathrm{Baryons}}(\ell)\rangle/\langle P_{\mathrm{Dark}}(\ell)\rangle$. _Top panel:_ suppression in both the full hydrodynamical simulation ($\kappa$TNG, in blue) and in our BCM implementation (in orange), finding good agreement between the two. _Bottom panel:_ the fractional difference between the two remains mostly below $\approx 1\%$. Both lines show averages over 2000 realizations, without noise added to the convergence maps. Figure 5: Peak histograms from noiseless $\kappa$TNG (_upper left panel_) and $\kappa$BCM (_upper right panel_). In both cases, we show the deviations in peak counts compared to $\kappa$Dark (_lower panels_). At ${\rm S/N}\gtrsim 4$ we find that the BCM underpredicts the number of peaks when compared to $\kappa$Dark maps, whereas TNG leads to a higher number of peaks.  Survey | Approx. $n_{\mathrm{gal}}\,[\mathrm{arcmin}^{-2}]$ | Approx. Area [$\mathrm{deg}^{2}$] | Approx. Median z
---|---|---|---
Roman | $50$ | $2200$ | $1.5$
LSST | $30$ | $18000$ | $1.0$
Euclid | $30$ | $15000$ | $1.0$
HSC | $20$ | $1500$ | $0.7$
DES | $10$ | $5000$ | $0.4$
KiDS | $10$ | $1350$ | $0.4$
Table 2: A summary of the approximate survey depth, area, and median redshift
we adopted to correspond to the six weak lensing surveys Roman, LSST, Euclid,
HSC, DES and KiDS.

In this section, we present results from weak lensing statistics compared
between $\kappa$TNG, $\kappa$Dark, and $\kappa$BCM. We first validate our BCM,
ray-tracing, and peak count implementations through comparisons to previous
work. We then discuss the peak count histograms in $\kappa$TNG and in
$\kappa$BCM.

As a validation step, we compute power spectra of each of the WL convergence
maps, and find the average power spectrum suppression, defined as

$S(\ell)\equiv\frac{P_{\mathrm{Baryons}}(\ell)}{P_{\mathrm{Dark}}(\ell)},$
(13)

for both Baryons $=$ TNG and Baryons $=$ BCM. Here $P(\ell)$ indicates the
average power spectrum over the 2000 realizations, measured using the package
nbodykit444https://github.com/bccp/nbodykit(Hand et al., 2018). The results in
Fig. 4 show that the average suppressions in the $\kappa$BCM and $\kappa$TNG
suites are very similar, with the absolute power spectra agreeing to within
$\sim 1\%$ at all $\ell$. On the other hand, the $\kappa$BCM power spectrum is
systematically slightly lower than $\kappa$TNG, which is likely caused by a
somewhat excessive AGN feedback component in the BCM. This effect has been
shown by A20 and in Lu et al. (2021) to dominate the power spectrum
suppression on small scales when using BCM’s. While in this work we use the
BCM parameters fit to the 3D matter power spectrum from A20, in the future we
will investigate directly fitting BCM parameters to other WL statistics. Here,
we find that the $\kappa$BCM’s difference from $\kappa$TNG is small and
consistent with previous studies, validating our implementation (Aricò et al.,
2021b; Lu & Haiman, 2021).

Next, for each smoothed convergence map, we identify the peaks and measure the
$\kappa_{\mathrm{rms}}$ values following § 2.4. In Fig. 5, we show the peak
count histogram obtained for $\kappa$TNG and $\kappa$BCM in the top panels, as
well as their comparison to $\kappa$Dark in the bottom panels. We plot these
histograms as functions of their significance (${\rm S/N}$) in $\kappa$ maps,
where ${\rm S/N}\equiv\kappa_{\rm peak}/\kappa_{\mathrm{rms}}$. Throughout
this work, we use approximate median redshifts, galaxy densities, and areas
found by inspecting details of LSST Science Collaboration et al. (2009); DES
Collaboration et al. (2021); Aihara et al. (2018); Laureijs et al. (2011); de
Jong et al. (2013); Spergel et al. (2015) (see Tab. 2). Fig. 5 serves two
purposes: (i) we are able to match our results to Osato et al. (2021),
validating our convergence map and peak count pipeline, and (ii) we can
compare the deviations between peaks (_lower panels_). We see qualitatively
that the differences between the $\kappa$TNG vs. $\kappa$Dark histograms are
well matched by $\kappa$BCM vs. $\kappa$Dark up to ${\rm S/N}\sim 4$. Above
this value, $\kappa$TNG has more peaks than $\kappa$Dark, yielding a positive
slope, whereas the BCM systematically underpredicts the number of peaks,
causing a negative slope. We explore the implications of this more in § 4.2

Next, we quantify how significant the deviations are between $\kappa$TNG and
$\kappa$BCM peaks. We consider surveys with depths ranging from
$10\,\mathrm{arcmin}^{-2}\leq n_{\mathrm{gal}}\leq 50\,\mathrm{arcmin}^{-2}$
and with median redshifts ranging from $0.4\leq z\leq 1.5$. We simplify our
analysis by placing all galaxies at a single redshift and add shape noise
following Equation (7). We smooth and compute the average $\kappa_{\rm rms}$
and peak count histogram over all $2000$ map realizations. To visually
illustrate the significance of the difference between the two histograms, we
consider the Poisson uncertainty in each bin, scaled to the appropriate survey
size. In Fig. 6, we then show the relative errors along with the Poisson
uncertainties, for depths, redshifts and solid angles roughly corresponding to
those in six different WL surveys. The parameters we adopted for these surveys
are listed in Table 2.

Figure 6: We compare the peak counts in $\kappa$BCM and $\kappa$TNG as a
function of peak height significance (S/N; bottom labels) and $\kappa$ (top
labels). The four panels show different combinations of median source galaxy
redshift and noise level, roughly corresponding to the redshifts and galaxy
number densities in different WL surveys, as labeled. We find that the
agreement between peaks from $\kappa$BCM and $\kappa$TNG is $<2\%$ up to ${\rm
S/N}\approx 4$, beyond which the discrepancies increase sharply. The shaded
regions correspond to $1\sigma$ Poisson uncertainties in each bin for DES,
HSC, Euclid, LSST, and Roman-like surveys (with the assumed areas listed in
Table 2), illustrating the significance of the deviations in peak counts. The
bottom panels show the contribution to the $\chi^{2}$ from each bin (see text
for a numerical calculation, taking into account correlations between bins).
For reference, we also show the (much larger) differences between $\kappa$TNG
and $\kappa$Dark (dashed curves).

Fig. 6 shows a $\sim 2\%$ agreement between $\kappa$TNG and $\kappa$BCM peak
statistics over the range $-1.5<{\rm S/N}<4$ for most galaxy densities
considered. The shaded regions in this figure mark the Poisson uncertainties
scaled to the rough areas of the DES, KiDs, HSC, Euclid, LSST, and Roman
surveys. We further plot the relative difference between $\kappa$TNG and
$\kappa$Dark, to show that the overall effect of the baryons is much larger
than the errors introduced by using the BCM. It is clear from this
visualization that $\kappa$BCM at lower redshifts and higher noise levels
achieves greater relative accuracy than at high redshift and with low noise.
Comparing the curves to the shaded areas also gives a rough estimate of the
significance of the deviations in each S/N bin. We further include the
corresponding $\kappa$ values which confirms a weak dependence of $\kappa_{\rm
rms}$ with respect to $n_{\rm gal}$ when compared to the the dependence of
$n_{\rm peaks}$ on $n_{\rm gal}$.

In the next section, we account quantitatively for correlations between bins
by computing the covariance matrix, the $\chi^{2}$ values, and the
corresponding probabilities. We also study how the significance of the
$\kappa$BCM vs. $\kappa$TNG difference is impacted by the presence of shape
noise.

$z$ | $\kappa_{\mathrm{rms}}^{\mathrm{TNG}}$ | $\kappa_{\mathrm{rms}}^{\mathrm{BCM}}$ | $\kappa_{\mathrm{rms}}^{\mathrm{Dark}}$
---|---|---|---
0.4 | 0.0040 | 0.0039 | 0.0041
0.7 | 0.0069 | 0.0068 | 0.0070
1.0 | 0.0099 | 0.0099 | 0.0100
1.5 | 0.0136 | 0.0136 | 0.0137
Table 3: The $\kappa_{\mathrm{rms}}$ values at the four different redshifts
with no noise used in this work. We show $\kappa_{\mathrm{rms}}$ for
$\kappa$TNG, $\kappa$BCM, and $\kappa$Dark to show how close they are, and to
justify the use of $\kappa_{\mathrm{rms}}^{\mathrm{TNG}}$ as the “noise” $N$
in determining the significance (S/N) of peaks in $\kappa$ maps. Each value
was computed by averaging the root-mean-square over all $2000$ convergence
maps.

##  4 Discussion

In the previous sections we described the process of modifying the TNG
simulations to account for baryons using the BCM of A20, how we ray-traced
through this simulation, and the statistics measured from the resulting maps.
The main result we found is that peak counts using the BCM agree at the
percent level with the full hydrodynamical TNG simulation up to a peak height
of ${\rm S/N}\approx 4$. In this section, we further quantify this result and
discuss its implications. We also look more closely at matches between
individual peaks at the map level, and assess whether the BCM’s accuracy is
sufficient for current and upcoming surveys.

###  4.1 Accuracy of BCM as a function of survey depth and area

The goal of using BCMs for WL peak count analysis is to avoid the large
computational costs associated with high-resolution hydrodynamic simulations
with baryonic physics. A natural question then is to what extent BCMs will be
a viable alternative in the analyses of current and future surveys? One could
either attempt to simultaneously fit cosmological and baryonic parameters, or
simply marginalize over the baryonic physics through the BCM parameters.

Fig. 6 partly addresses this question. It shows that the systematic deviation
in peak counts introduced by the BCM stays within $2\%$ up to peak heights of
${\rm S/N}\approx 4$. Furthermore, this is within the Poisson uncertainty in
each bin, except for the largest (nearly all-sky) surveys. The shaded region
in this plot represents the square root of the covariance matrix’s diagonal
elements. The ratio of the $\mathbf{\Delta n}$ curves to the width of the
shaded area is therefore the contribution to the $\chi^{2}$ from each bin
(corresponding to the diagonal elements in the full matrix product Eq. (9)).
These are shown in the bottom panels, and reveal that the largest
contributions to $\chi^{2}$ are from peaks with ${\rm S/N}\gtrsim 4$ for all
redshifts and noise levels. However, we note that for LSST and Euclid, at
$z=1$ and $n_{\mathrm{gal}}=30$, significant contributions to $\chi^{2}$ come
from smaller peaks ($0\lesssim{\rm S/N}\lesssim 4$) as well.

We next include correlations between bins, i.e. compute the full $\chi^{2}$
including the off-diagonal terms of the covariance matrix in Eq. (9). We
compute peak histograms at each of the four redshifts $z=0.4,0.7,1$ and $1.5$
used in this work, roughly corresponding to the median redshift of current and
future surveys. At each redshift, we add shape noise corresponding to
$n_{\mathrm{gal}}=10,20,30,50\,\mathrm{arcmin}^{-2}$, respectively, which
roughly represent the depths of these surveys (see Table 2). All galaxies are
placed at a single source plane at the median redshift for simplicity. We
adopt the same binning as in § 2.4 and follow the procedure outlined in that
section to compute $\chi^{2}$ for each combination of noise level, source
redshift, and survey area. Note that $\chi^{2}$ is simply proportional to
survey size.

In Fig. 7 we show a map of the resulting $\chi^{2}$ values as a function of
depth ($n_{\rm gal}$) and survey size. The lines mark the 1, 2, and 3$\sigma$
contours, and the four panels correspond to increasing redshifts, with the
relevant surveys marked in each panel areas of KiDS, DES, HSC, Euclid, and
Roman. This plot shows that the BCM is sufficient to use for most surveys
considered in this work. In particular, at the depth/area combinations of the
existing/ongoing surveys DES, KiDs and HSC, the BCM predictions are well
within the 1$\sigma$ uncertainties. The deepest survey we consider, Roman,
still remains within the 2$\sigma$ uncertainty. On the other hand, in the two
largest surveys, LSST and Euclid, the $\kappa$BCM peak counts become
discrepant from $\kappa$TNG well beyond the $3\sigma$ level. This is clearly
due to the large number of peaks in these surveys, which correspondingly
strongly reduces the statistical errors on the peak counts. In the next
section, we attempt to alleviate these discrepancies between $\kappa$TNG and
$\kappa$BCM by removing bins from consideration.

Figure 7: The significance of the difference between the peak counts in the
$\kappa$BCM and the full hydrodynamical $\kappa$TNG maps, as a function of
shape noise ($n_{\rm ngal}$ and survey area. The four panels roughly
correspond to the median redshifts in existing/ongoing and future WL surveys,
as labeled in each panel. The symbol marking each survey also indicates the
approximate depth and area of that survey. The three curves in each panel mark
the 1, 2, and 3$\sigma$ uncertainties of the $\kappa$TNG model. The
$\kappa$BCM predictions are within these uncertainties of most surveys, with
the exception of LSST and Euclid. The small statistical errors in these two
large future surveys render the $\kappa$BCM predictions discrepant well beyond
the 3$\sigma$ level.

###  4.2 Differences between $\kappa$BCM and $\kappa$TNG peak counts

We next take a closer look at how the population of peaks in $\kappa$BCM and
$\kappa$TNG differ, on the map level.

Reexamining Fig. 2, we find regions in the leftmost panel containing multiple
high density peaks (for example, the region bounded by $0.5\leq X\leq
1.5\,\mathrm{deg}$ and $3\leq Y\leq 4\,\mathrm{deg}$). These peaks appear to
coincide with alternating _blue-red_ “dipole” patterns in the right two
panels. A natural explanation for this pattern is that baryonic effects shift
the location of the peaks on the $\kappa$TNG maps, compared to $\kappa$Dark.
The locations of haloes in dark-matter only vs. hydrodynamical simulations are
generally different. Since the BCMs change only the halo profiles, but not
their locations, they are unable to capture this baryonic effect.

Figure 8: A subset of peaks, in a zoom-in of a small region, of the difference
$\kappa$ map ($\kappa$TNG-$\kappa$BCM, where red is positive and blue is
negative) in the rightmost panel of Fig. 2. The size of each circle (square)
is proportional to the height of the peak in $\kappa$TNG ($\kappa$BCM). The
figure illustrates the small offset between the peak centres in $\kappa$TNG
and $\kappa$BCM evident in $>25\%$ of the peaks, causing an alternating colour
“dipole” pattern. We interpret this as being caused by different halo centre
locations in the TNG and BCM density fields (see $X\approx 1.15$, $Y\approx
3.5$). The BCM cannot reproduce baryonic effects well outside haloes, or the
impact of baryonic effects on the halo locations.  Figure 9: The difference
between peak heights for $5$ $\kappa$BCM and $\kappa$TNG maps (corresponding
to $3337$ peaks) at $z=1.5$ plotted against the peak heights from $\kappa$TNG.
We see a systematic under-prediction by $\kappa$BCM of peak heights at ${\rm
S/N}\gtrsim 4$.

To test this hypothesis further, in Fig. 8 we show the population of peaks in
a small sub-region of the convergence maps from $\kappa$TNG and $\kappa$BCM.
The size of each circle (square) is proportional to the $\kappa$ value of the
peak in $\kappa$TNG ($\kappa$BCM). The background colours reproduce the
rightmost panel of Fig. 2, representing $\kappa$TNG-$\kappa$BCM. Fig. 8 shows
that most of the peaks exist in both $\kappa$TNG and $\kappa$BCM, and have
similar $\kappa$ values. On the other hand, the heights do not match
perfectly, and there are two $\kappa$TNG peaks that are missed in $\kappa$BCM.
Furthermore, the centres of many peaks are visibly offset. A clear example of
this is the peak at $X\approx 1.3\,\mathrm{deg}$, $Y\approx
3.55\,\mathrm{deg}$. The square ($\kappa$BCM peak) is shifted to the left,
while the circle ($\kappa$TNG peak) is shifted to the right. This matches the
dipole colour pattern of the background, which shows a $\kappa$BCM-dominated
region to the left of a $\kappa$TNG-dominated region. When comparing the
locations of peaks in general, we find that $\gtrsim 25\%$ of the peaks in
$\kappa$TNG and $\kappa$BCM have offsets by at least one pixel.

A plausible origin of this effect is that the BCM is unable to account for the
formation histories and locations of dark matter haloes. While baryonic
effects, such as AGN feedback, are implemented at a single snapshot in space
and time by the BCM, these effects are not propagated through the simulation,
resulting in different halo centres. Additionally, baryonic processes can
modify the density distribution well outside haloes, as well, beyond the
regions captured by the BCMs. Since non-halo matter is known to contribute to
peak counts (Sabyr et al., 2021), this effect can also contribute to changing
peak locations and heights, when compared to dark-matter only simulations.

Beyond peak placement, the BCM assigns systematically lower $\kappa$ values to
peaks. As shown in the right panels of Fig. 5 and Fig. 6, the number of peaks
in the highest S/N regions are underestimated by the BCM when compared to
either $\kappa$Dark or $\kappa$TNG maps. Further evidence for this can be seen
in Fig. 8. Regions with high peaks that share the same peak centre (e.g.
$X\approx 1.18\,\mathrm{deg}$, $Y=3.3\,\mathrm{deg}$) tend to be redder,
corresponding to a larger $\kappa$TNG peak height. This is likely due to the
BCM over-accounting for AGN feedback, an effect which we commented on in § 3.
This effect would disperse more matter out of haloes, reducing their density
and lower the peak $\kappa$ values. As the peak-count histograms show, their
slope becomes steeper at higher S/N. Because of this, we expect that the
number of these high peaks (${\rm S/N}\gtrsim 4$), are the most affected by
moving peaks into lower-$\kappa$ bins.

This motivates us to investigate the peak heights more carefully. We find all
the peaks that are present in both $\kappa$TNG and $\kappa$BCM for 10
individual maps, where we identify the matching between the peaks in both maps
by their distance being smaller than a threshold of $2.9\,{\rm arcmin}$. In
Fig. 9 we compare peak heights from $\kappa$TNG and $\kappa$BCM. From this
representation it is clear that the BCM systematically assigns lower peak
heights on average at ${\rm S/N}\gtrsim 4$. Intuitively, the increased AGN
feedback present in BCM will manifest as smaller peak heights for the peaks
that are made up of individual massive haloes and which have expelled an
overestimated amount of matter. For smaller peaks which rely on multiple low-
mass haloes and are therefore less sensitive to to this effect, we see that
the systematic underestimation in the peak heights of the BCM are less severe
or non-existent.

We also note that there is a small number of peaks in the $\kappa$TNG maps
that do not have a clear match in $\kappa$BCM and vice versa. More
specifically, the number of unmatched peaks that we had to exclude in Fig. 2
is roughly 2% of the total $\kappa$TNG peaks. Likewise, there are peaks in
$\kappa$BCM without a match in $\kappa$TNG, which we find also account for
about 2% of the total. The asymmetry in number among this population of
unmatched peaks is small, and, furthermore, the missed peaks have a maximum
S/N of 3.5. We therefore conclude that the reason for BCM’s underprediction of
the number of high S/N peaks is the systematic underestimation of the peak
heights, rather than missing peaks entirely.

###  4.3 Mitigation of high S/N peak discrepancies

We have seen that the BCM is unable to reproduce hydrodynamical simulations at
the pixel level, and that peak counts between $\kappa$TNG and $\kappa$BCM
deviate at ${\rm S/N}\gtrsim 4$. Further we find that future surveys with
large areas such as LSST and Euclid produce $\chi^{2}$ differences well above
the $1\sigma$ random coincidence level.

With this motivation, we consider a remedy to lower the value of $\chi^{2}$
for LSST and Euclid from $\chi^{2}_{\mathrm{LSST}}=77.98$ and
$\chi^{2}_{Euclid}=64.97$ to the $1\sigma$ threshold $\chi^{2*}_{\rm D.O.F}$
which depends on the number of bins. To do this, we remove the most discrepant
peaks and use only the remaining subset of histogram bins. Since both the
highest and lowest peaks are highly discrepant, we consider various
combinations of lower-S/N and upper-S/N cuts in the $\chi^{2}$ computation. In
Fig. 10, we show the maximum allowable survey area corresponding to $1\sigma$
$\chi^{2*}_{\rm D.O.F}$ as a function of these S/N thresholds. Each pixel
corresponds to a S/N cut, with the bottom left representing utilizing the full
$18$ bins, and the diagonal elements containing only a single bin. _Green_ and
_blue_ pixels show cuts that generate a $1\sigma$ maximum allowable survey
area greater than or equal to Euclid and LSST respectively. As the figure
shows, additional bins can either decrease or increase the allowed area. This
is because bins in which BCM is accurate decrease $\chi^{2}$ with respect to
$\chi^{2*}_{\rm D.O.F}$ (accounting for the extra D.O.F) , while discrepant
bins cause an increase. This figure shows that reducing $\chi^{2}$ to the
$1\sigma$ level for LSST and Euclid is not possible by just excluding the
highest and/or lowest peaks. The $1\sigma$ criterion can be met only by
removing the vast majority of peaks from the analysis, which would mean a
significant loss of cosmological sensitivity. This result could also be
anticipated from the bottom panel in the $n_{\rm gal}=30$, $z=1.0$ case in
Fig. 6. This shows that the discrepancy between $\kappa$BCM and $\kappa$TNG is
largest at ${\rm S/N}\gtrsim 4$, but the lower-S/N bins still contribute
significantly.

Figure 10: We compute the maximum survey area that would produce $1\sigma$
deviations between $\kappa$TNG and $\kappa$BCM, when the lowest and/or highest
S/N peaks are excluded from the $\chi^{2}$ calculation. Each pixel represents
a different combination of these cuts given by the x-axis lower bound and
y-axis upper bound. We highlight the regions above Euclid’s area (_green_),
LSST’s area (_blue_) and all smaller areas (_red_).

As we saw in Figs. 4 and 9, there is a systematic under-prediction of peak
heights and power spectra. We speculated that this was due to an
overestimation of the AGN feedback component, thus removing too much matter
from the most massive haloes. This suggests a different mitigation approach:
to change the BCM model parameters or possibly modify the BCM model itself.
Aricò et al. (2021a) shows that fitting the BCM parameters jointly to the
power spectrum and the bispectrum provides a vast improvement to the small-
scale bispectrum in hydrodynamical simulations, compared to using BCMs fit
only to the power spectrum. Furthermore, they find that this joint fitting
does not significantly degrade the fit to the power spectrum. In future work,
we plan to fit the BCM parameters to the peaks themselves, which we expect
would increase the agreement between peak heights at the largest peaks.
Another possible alternative is to change the BCM model to provide a more
flexible AGN feedback component that is reduced for the more massive halos,
removing the systematic under-prediction for the largest peaks.

##  5 Summary and conclusions

In this work, we investigated the level of agreement in weak lensing peak
counts between semi-analytical baryonic correction models (BCMs), and
hydrodynamic simulations. We applied a BCM developed in Aricò et al. (2020),
which has been shown to match 3D matter power spectra and lensing convergence
power spectra in IllustrisTNG, and other hydrodynamic simulations. After
verifying that our implementation and parameterisation is comparable to the
results of Aricò et al. (2020), we performed ray-tracing as developed in Osato
et al. (2021), and computed the peak statistics from the resulting convergence
($\kappa$) maps.

We find that the ray-traced BCM maps ($\kappa$BCM) is able to reproduce the
peak statistics from ray-traced TNG maps ($\kappa$TNG) at $\sim$ percent
accuracy for a range of redshifts, survey depths, and sky areas up to a peak
height of ${\rm S/N}\approx 4$, but significantly under-predicts the number of
peaks above this threshold. We investigated the pixel-level agreement between
$\kappa$BCM and $\kappa$TNG maps, and found that while $\kappa$BCM matches
almost all of the $\kappa$TNG peaks, it yields systematically lower peak
heights, misses a small fraction of the peaks, and yields a small offset in
the peak locations. We believe that the latter effect is due to TNG capturing
the impact of baryons well outside haloes and halo formation histories
affected by baryons that the BCM cannot account for. The most significant
effect for the peak-count statistic arises from $\kappa$BCM’s underestimation
of the heights of the largest-S/N peaks. We speculate that this is caused by
the BCMs over-accounting for AGN feedback and removal of excess mass from
haloes, and suggest that this could be corrected by further modifying the
BCMs.

We investigated the statistical significance between the $\kappa$BCM and
$\kappa$TNG peak histograms, using a simple $\chi^{2}$ statistic. We found
that deviations between $\kappa$BCM and $\kappa$TNG are statistically
insignificant at the depths and sky coverages typical of current and on-going
surveys (well within $1\sigma$ uncertainties). On the other hand, the
discrepancy rises to the $\approx 2\sigma$ level for a deep survey such as
Roman, and becomes even more significant (well above $3\sigma$) for large
survey areas, such as in the surveys by LSST and Euclid. We considered a
simple remedy in which we excluded the most discrepant lowest and/or highest
peaks from the $\chi^{2}$ computation, but found that this was insufficient in
reducing $\chi^{2}$ to the level of $1\sigma$.

In this work, we did not perform any optimization by parameter fitting for the
BCM using the peaks, and instead adopted the values presented by Aricò et al.
(2020) that best fit the 3D matter power spectrum. We will investigate the
differences in parameters when fitting power spectra and peaks in a future
work, to assess whether a single set of parameters simultaneously fits both
statistics, or if peak-count predictions could be further improved through
calibration of BCM parameters. The influence of BCMs on other non-Gaussian
statistics, such as WL minima and Minkowski functionals, as well as using more
realistic survey parameters (such as redshift distributions) are also of
interest, and will be also explored.

Most importantly, in this work, we studied only the statistical difference
between peaks counts in $\kappa$BCM and $\kappa$TNG. In future work, we will
assess the impact of BCM biases on cosmological parameter estimation. Our
results presented here are somewhat conservative, in the sense that when the
peak counts are statistically indistinguishable, they can not bias
cosmological parameter-inference. On the other hand, when the $\kappa$BCM peak
counts statistically differ from $\kappa$TNG, they may or may not cause a
significant cosmological-parameter bias, depending on how the cosmological vs.
BCM parameters impact the peak counts.

## Acknowledgements

We acknowledge support by NASA ATP grant 80NSSC18K1093 (to ZH), the use of the
NSF XSEDE facility Stampede2, and the Columbia University High-Performance
Computing cluster Ginsburg for simulations and data analysis used in this
study. KO is supported by JSPS Research Fellowships for Young Scientists. This
work was supported by Grant-in-Aid for JSPS Fellows Grant Number JP21J00011.

## Data Availability

The data underlying this article were accessed from Stampede2. The derived
data generated in this research will be shared on reasonable request to the
corresponding author.

## References

  * Aihara et al. (2018) Aihara H., et al., 2018, PASJ, 70, S4
  * Aricò et al. (2020) Aricò G., Angulo R. E., Hernández-Monteagudo C., Contreras S., Zennaro M., Pellejero-Ibañez M., Rosas-Guevara Y., 2020, MNRAS, 495, 4800
  * Aricò et al. (2021a) Aricò G., Angulo R. E., Hernández-Monteagudo C., Contreras S., Zennaro M., 2021a, MNRAS, 503, 3596
  * Aricò et al. (2021b) Aricò G., Angulo R. E., Contreras S., Ondaro-Mallea L., Pellejero-Ibañez M., Zennaro M., 2021b, MNRAS, 506, 4070
  * Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Cheng et al. (2020) Cheng S., Ting Y.-S., Ménard B., Bruna J., 2020, MNRAS, 499, 5902
  * Chisari et al. (2019) Chisari N. E., et al., 2019, The Open Journal of Astrophysics, 2, 4
  * Coulton et al. (2019) Coulton W. R., Liu J., Madhavacheril M. S., Böhm V., Spergel D. N., 2019, J. Cosmology Astropart. Phys., 2019, 043
  * DES Collaboration et al. (2021) DES Collaboration et al., 2021, arXiv e-prints, p. arXiv:2105.13549
  * Dai & Seljak (2021) Dai B., Seljak U., 2021, Proceedings of the National Academy of Science, 118, 2020324118
  * Dalcin & Fang (2021) Dalcin L., Fang Y.-L. L., 2021, Computing in Science Engineering, 23, 47
  * Dietrich & Hartlap (2010) Dietrich J. P., Hartlap J., 2010, Monthly Notices of the Royal Astronomical Society, 402, 1049
  * Fluri et al. (2019) Fluri J., Kacprzak T., Lucchi A., Refregier A., Amara A., Hofmann T., Schneider A., 2019, Phys. Rev. D, 100, 063514
  * Fu et al. (2014) Fu L., et al., 2014, MNRAS, 441, 2725
  * Gatti et al. (2021) Gatti M., et al., 2021, arXiv e-prints, p. arXiv:2110.10141
  * Gupta et al. (2018) Gupta A., Zorrilla Matilla J. M., Hsu D., Haiman Z., 2018, Phys. Rev. D, 97, 103515
  * Halder & Barreira (2022) Halder A., Barreira A., 2022, arXiv e-prints, p. arXiv:2201.05607
  * Halder et al. (2021) Halder A., Friedrich O., Seitz S., Varga T. N., 2021, MNRAS, 506, 2780
  * Hamana et al. (2020) Hamana T., et al., 2020, PASJ, 72, 16
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Harnois-Déraps et al. (2021) Harnois-Déraps J., Martinet N., Castro T., Dolag K., Giblin B., Heymans C., Hildebrandt H., Xia Q., 2021, MNRAS, 506, 1623
  * Hartlap et al. (2007) Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399
  * Heymans et al. (2013) Heymans C., et al., 2013, MNRAS, 432, 2433
  * Hikage et al. (2019) Hikage C., et al., 2019, Publications of the Astronomical Society of Japan, 71
  * Hilbert et al. (2009) Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A, 499, 31
  * Hildebrandt et al. (2017) Hildebrandt H., et al., 2017, MNRAS, 465, 1454
  * Hoekstra & Jain (2008) Hoekstra H., Jain B., 2008, Annual Review of Nuclear and Particle Science, 58, 99
  * Jain & Van Waerbeke (2000) Jain B., Van Waerbeke L., 2000, ApJ, 530, L1
  * Jain et al. (2000) Jain B., Seljak U., White S., 2000, ApJ, 530, 547
  * Jing et al. (2006) Jing Y. P., Zhang P., Lin W. P., Gao L., Springel V., 2006, ApJ, 640, L119
  * Kacprzak et al. (2016) Kacprzak T., et al., 2016, MNRAS, 463, 3653
  * Kilbinger (2015) Kilbinger M., 2015, Reports on Progress in Physics, 78, 086901
  * Kitching et al. (2014) Kitching T. D., et al., 2014, MNRAS, 442, 1326
  * Kratochvil et al. (2010) Kratochvil J. M., Haiman Z., May M., 2010, Phys. Rev. D, 81, 043519
  * Kratochvil et al. (2012) Kratochvil J. M., Lim E. A., Wang S., Haiman Z., May M., Huffenberger K., 2012, Phys. Rev. D, 85, 103513
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Li et al. (2019) Li Z., Liu J., Zorrilla Matilla J. M., Coulton W. R., 2019, Phys. Rev. D, 99, 063527
  * Liu & Haiman (2016) Liu J., Haiman Z., 2016, Phys. Rev. D, 94, 043533
  * Liu et al. (2015) Liu J., Petri A., Haiman Z., Hui L., Kratochvil J. M., May M., 2015, Phys. Rev. D, 91, 063507
  * Lu & Haiman (2021) Lu T., Haiman Z., 2021, Monthly Notices of the Royal Astronomical Society, 506, 3406–3417
  * Lu et al. (2021) Lu T., Haiman Z., Zorrilla Matilla J. M., 2021, arXiv e-prints, p. arXiv:2109.11060
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Marques et al. (2019) Marques G. A., Liu J., Zorrilla Matilla J. M., Haiman Z., Bernui A., Novaes C. P., 2019, J. Cosmology Astropart. Phys., 2019, 019
  * Martinet et al. (2021) Martinet N., Harnois-Déraps J., Jullo E., Schneider P., 2021, A&A, 646, A62
  * Matilla et al. (2020) Matilla J. M. Z., Sharma M., Hsu D., Haiman Z., 2020, Phys. Rev. D, 102, 123506
  * Maturi et al. (2010) Maturi M., Angrick C., Pace F., Bartelmann M., 2010, A&A, 519, A23
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * Munshi et al. (2011) Munshi D., van Waerbeke L., Smidt J., Coles P., 2011, Monthly Notices of the Royal Astronomical Society, 419, 536
  * Naiman et al. (2018) Naiman J. P., et al., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1206–1224
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, The Astrophysical Journal, 462, 563
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Osato et al. (2015) Osato K., Shirasaki M., Yoshida N., 2015, ApJ, 806, 186
  * Osato et al. (2021) Osato K., Liu J., Haiman Z., 2021, MNRAS, 502, 5593
  * Petri et al. (2013) Petri A., Haiman Z., Hui L., May M., Kratochvil J. M., 2013, Phys. Rev. D, 88, 123002
  * Petri et al. (2016) Petri A., Haiman Z., May M., 2016, Phys. Rev. D, 93, 063524
  * Pillepich et al. (2017) Pillepich A., et al., 2017, Monthly Notices of the Royal Astronomical Society, 475, 648
  * Ribli et al. (2019) Ribli D., Pataki B. Á., Zorrilla Matilla J. M., Hsu D., Haiman Z., Csabai I., 2019, MNRAS, 490, 1843
  * Sabyr et al. (2021) Sabyr A., Haiman Z., Zorrilla Matilla J. M., Lu T., 2021, arXiv e-prints, p. arXiv:2109.00547
  * Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  * Schneider & Bridle (2010) Schneider M. D., Bridle S., 2010, MNRAS, 402, 2127
  * Schneider et al. (2019) Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmology Astropart. Phys., 2019, 020
  * Sifón et al. (2015) Sifón C., Hoekstra H., Cacciato M., Viola M., Köhlinger F., van der Burg R. F. J., Sand D. J., Graham M. L., 2015, A&A, 575, A48
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Springel (2010) Springel V., 2010, Monthly Notices of the Royal Astronomical Society, 401, 791–851
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, arXiv e-prints, p. arXiv:2010.00619
  * Yang et al. (2011) Yang X., Kratochvil J. M., Wang S., Lim E. A., Haiman Z., May M., 2011, Phys. Rev. D, 84, 043529
  * Zürcher et al. (2021) Zürcher D., Fluri J., Sgier R., Kacprzak T., Refregier A., 2021, J. Cosmology Astropart. Phys., 2021, 028
  * de Jong et al. (2013) de Jong J. T. A., Verdoes Kleijn G. A., Kuijken K. H., Valentijn E. A., 2013, Experimental Astronomy, 35, 25