# A deconstruction of methods to derive one-point lensing statistics

Viviane Alfradique  1,2111These authors contributed equally to this work.
Tiago Castro  0000-0002-5225-1923 3,4,5,6**footnotemark: * Valerio Marra
0000-0002-6292-3228 3,4,7 Miguel Quartin  0000-0002-7773-1579 1,8,9
Carlo Giocoli10,11 and Pierluigi Monaco3,4,6,12 0000-0001-5853-6164
0000-0002-9590-7961 0000-0003-2083-7564 1Instituto de Física, Universidade
Federal do Rio de Janeiro, 21941-972, Rio de Janeiro, RJ, Brazil
2Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180
Rio de Janeiro, RJ, Brazil
3INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste,
Italy
4INFN – Sezione di Trieste, I-34100 Trieste, Italy
5IFPU – Institute for Fundamental Physics of the Universe, via Beirut 2,
34151, Trieste, Italy
6ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e
Quantum Computing, Via Magnanelli 2, Bologna, Italy
7Departamento de Física, Universidade Federal do Espírito Santo, 29075-910,
Vitória, ES, Brazil
8PPGCosmo, Universidade Federal do Espírito Santo, 29075-910, Vitória, ES,
Brazil
9Observatório do Valongo, Universidade Federal do Rio de Janeiro, 20080-090,
Rio de Janeiro, RJ, Brazil
10INAF - Osservatorio Astronomico di Bologna, via Piero Gobetti 93/3, I-40129
Bologna, Italy
11INFN, Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy
12Università di Trieste, Dipartimento di Fisica, via Tiepolo 11, 34143 Trieste

###### Abstract

Gravitational lensing is a crucial tool for exploring cosmic phenomena,
providing insights into galaxy clustering, dark matter, and dark energy. Given
the substantial computational demands of $N$-body simulations, approximate
methods like PINOCCHIO and turboGL offer viable alternatives for simulating
lensing probability density functions (PDFs). This paper evaluates these
methods in contexts where baryonic effects are negligible, focusing on dark
matter-dominated models and assessing their effectiveness across both weak and
strong lensing regimes. Our comparative analysis reveals that these methods
are particularly effective for applications involving electromagnetic and
gravitational wave point sources, where strong lensing events are infrequent.
Both PINOCCHIO and turboGL perform well in modeling the weak-lensing region
influenced by mildly nonlinear structures. However, they lose accuracy in
capturing small-scale nonlinear matter fields, owing to oversimplified
assumptions about internal halo structures and reliance on perturbation
theory. The analysis shows that PINOCCHIO achieves an 8–15% agreement with
$N$-body simulations for the second-to-fourth moments of lensing PDFs. These
findings aim to inform future studies on gravitational lensing of point
sources, which are increasingly relevant with upcoming supernova and
gravitational wave datasets.

###### keywords:

large-scale structure of Universe – cosmology:observations – cosmological
parameters

##  1 Introduction

Electromagnetic and gravitational wave point sources, such as supernovae and
compact binary coalescences, have proven to be invaluable for studying the
universe’s evolution across extensive redshift ranges. Upcoming observatories
like the Euclid mission,222https://www.euclid-ec.org/ the Nancy Grace Roman
Space Telescope,333https://roman.gsfc.nasa.gov/ and the Vera C. Rubin
Observatory444https://www.lsst.org/ are set to accurately observe tens of
thousands of supernovae up to a redshift of $z\sim 2$. Simultaneously, third-
generation gravitational wave detectors such as the Einstein
Telescope555https://www.et-gw.eu and Cosmic
Explorer666https://cosmicexplorer.org will detect hundreds of thousands of
gravitational events up to $z\sim 5$ and beyond. To harness this wealth of
data effectively, it is crucial to meticulously model all potential systematic
effects, with gravitational lensing being a significant concern, particularly
for sources at $z\gtrsim 1$, as it strongly influences distance measurements.

Gravitational lensing affects electromagnetic and gravitational waves as they
traverse our inhomogeneous Universe. This phenomenon describes the deflections
that waves undergo when interacting with the matter density field, causing
distant sources to appear magnified or demagnified. Such deflections can lead
to significant distance measurement errors, modifying the intrinsic
distribution of astronomical sources and introducing variance and non-
Gaussianities. Furthermore, this deflections also gives rise to distortions in
the shapes of the background galaxies. A gravitational lens induces a
correlation on galaxy image orientations which can be measured statically and
reveal the properties and evolution of large-scale structures. This is known
as cosmic shear (see e.g. Bartelmann and Schneider, 2001; Kilbinger, 2015, for
a review).

Gravitational lensing notably limits the accuracy of distant standard candles
and sirens as cosmological probes. Estimates suggest that lensing adds a
relative uncertainty of 0.5%–1% to distance measurements of stellar binary
black holes up to 3 Gpc (Shan et al., 2021). Recent studies, such as those by
Mpetha et al. (2024), have indicated that neglecting this effect can bias
cosmological analyses of bright sirens at $z\leq 0.5$ by more than 1$\sigma$,
resulting in $\Delta H_{0}=\pm$ 0.25 km/s/Mpc for binary neutron stars (BNS)
and $\Delta H_{0}=\pm$ 5 km/s/Mpc for supermassive black hole binaries.
Canevarolo and Chisari (2023) further demonstrated that gravitational lensing
could introduce biases up to 3.5$\sigma$ for the $\Omega_{\rm m0}$ measurement
and up to 2.5$\sigma$ for $H_{0}$ using BNS up to $z=2$, and Canevarolo et al.
(2024) that lensing may also bias the distribution of neutron star masses.
Moreover, Wu et al. (2023) showed that ‘delensing’ methods, which correct the
effect after estimating the lensing magnification independently through weak
lensing maps reconstructed from galaxy surveys, could halve these lensing
errors on average (Shapiro et al., 2010).

Accurate modeling of gravitational lensing effects on the Hubble diagram
residuals can yield significant cosmological insights. By meticulously
addressing the systematic effects of lensing, researchers can effectively
transform what appears as noise into a valuable signal. Numerous studies
(Quartin et al., 2014; Castro and Quartin, 2014; Castro et al., 2016b, a;
Macaulay et al., 2017, 2020; Scovacricchi et al., 2017) have investigated the
correlation between the Hubble diagram residuals and lensing signals,
examining how these relationships depend on the background cosmology and the
perturbations in intervening matter. These efforts have demonstrated the
potential to achieve a precision of 0.6% for the amplitude of the power
spectrum $\sigma_{8}$ and 7% for the growth rate index $\gamma$ in forthcoming
surveys (Amendola et al., 2015).

The lensing probability distribution function (PDF) is the fundamental one-
point statistic for describing the lensing of point sources. The most direct
and precise method for determining the lensing PDF involves applying the ray-
tracing technique in full hydrodynamical cosmological simulations (Castro et
al., 2018). However, these simulations demand extensive computational
resources. In regimes with negligible baryonic effects on lensing, $N$-body
simulations using collisionless particles provide a sufficient alternative and
significantly reduce computational demands. Despite the reliability of the
results from $N$-body simulations (Takahashi et al., 2011; Hilbert et al.,
2020), the computational requirements and storage capacities remain
challenging.

An even faster approach, which combines precision and speed, is to use
approximate methods to describe the nonlinear evolution of perturbations
(Manera et al., 2012; Kitaura and S., 2013; Tassev et al., 2013; Monaco,
2016). The PINpointing Orbit Crossing Collapsed HIerarchical Objects
(PINOCCHIO) algorithm (Monaco et al., 2002; Taffoni et al., 2002; Munari et
al., 2016) integrates Lagrangian perturbation theory (LPT) with a semi-
analytical model to approximate the nonlinear hierarchical formation of
structures. This method generates dark matter halo catalogs using a fraction
of the resources needed for $N$-body simulations and achieves roughly 10%
accuracy in clustering statistics (power spectrum, bispectrum, and two-point
correlation function). For the weak-lensing regime, the stochastic
gravitational lensing (sGL) approach provides an even faster alternative
(Kainulainen and Marra, 2009, 2011a, 2011b). This method calculates the 1-halo
term of the lensing PDF using stochastic configurations of inhomogeneities
based on a halo model and subsequently incorporates the 2-halo term through
convolution. A similar analytical method was recently proposed by Thiele et
al. (2020).

In this paper, we deconstruct the essential physics needed to calculate
lensing PDFs accurately. We explore the effect of the 2-halo term, the
influence of non-collapsed matter in the field, the accuracy of PINOCCHIO, and
the degree to which baryonic effects can be overlooked. We specifically
analyze four relevant lensing quantities in one-point statistics: the PDFs for
convergence ($\kappa$), shear ($\gamma$), magnification ($\mu$), and magnitude
magnification ($\Delta m$), the latter being particularly crucial for the
statistics of lensing of standard candles like supernovae and gravitational
waves. We will describe our methodology in Section 2, present our comparative
results in Section 3, and summarize our findings in Section 4. In Appendix A,
we revise the impact of the angular resolution on the lensing PDF.

##  2 Lensing Simulation Methods

Table 1: Simulation methods ordered according to their complexity and the range of physical phenomena modeled. model name |  description
---|---
tgl 1h |  Stochastic gravitational lensing simulation with the 1-halo turboGL code.
tgl 1h+2h |  Extends tgl 1h by including a theoretical 2-halo term convolution.
PINOCCHIO${}_{\text{halos}}$ |  Ray-tracing on PINOCCHIO past light cone (PLC) with field particles randomly redistributed.
PINOCCHIO |  Ray-tracing on standard PINOCCHIO PLC.
$N$-body |  Ray-tracing on $N$-body simulation PLC.
Hydro |  Ray-tracing on hydro-dynamical simulation PLC.

This study aims to assess the accuracy of various simulation methods in
deriving the lensing Probability Density Function (PDF) on the regime in which
baryonic physics is sub-dominant. The methods under investigation are
summarized in Table 1, ranked by their complexity and the breadth of physical
phenomena they encompass. For the PINOCCHIO and $N$-body simulations, we
further evaluate the angular power spectrum and the generated lensing and
density maps. The Hydro simulations are used simply to define the boundaries
of reliability of the $N$-body results. For all methods, our fiducial model
for the Universe is the $\Lambda$CDM model with the cosmological parameters
set very close to the ones obtained by the combination of DES Y1 + Planck 2015
+ JLA SNe + BAO (Abbott et al., 2018): $\Omega_{\rm m0}=0.301,\,\Omega_{\rm
b0}=0.048,\,\Omega_{\Lambda 0}=0.699,\,h=0.682,\,\sigma_{8}=0.798$, and
$n_{\rm s}=0.973$. We proceed to describe each method in detail.

###  2.1 The sGL method

The turboGL (tgl) code777https://github.com/valerio-marra/turboGL is the
numerical implementation of the semi-analytical stochastic gravitational
lensing (sGL) method, which was introduced in Kainulainen and Marra (2009,
2011a). The sGL method operates by modeling the matter density contrast
according to the halo model (HM), where the inhomogeneous Universe is
approximated as a collection of different types of halos whose positions
satisfy the linear power spectrum.

####  2.1.1 The halo model

The halo model assumes that on small scales (large wavenumbers $k$) the
statistics of matter correlations are dominated by the internal halo density
profiles. In contrast, on large scales, the halos are assumed to cluster
according to linear theory. The model does not include intermediate low-
density structures such as filaments and walls. Following Kainulainen and
Marra (2011b), in tgl the total power spectrum is obtained by the simple
addition of the two components:888More sophisticated and accurate versions of
the HM are available in the literature; see Asgari et al. (2023), and
references therein.

$P_{m}(k,z)=P_{L}(k,z)+P_{H}(k,z)\,.$ (1)

The first term on the right-hand side is also called the 2-halo component, and
the second is the 1-halo component. The usefulness of the halo model stems
from the fact that both terms in Eq. (1) can be computed without having to
resort to numerical simulations. Specifically (Peacock and Smith, 2000):

$\displaystyle P_{L}(k,z)$ $\displaystyle={2\pi^{2}\over
k^{3}}\delta_{H0}^{2}\left({ck\over
H_{0}}\right)^{3+n_{s}}T^{2}(k)D^{2}(z)\,,$ (2) $\displaystyle P_{H}(k,z)$
$\displaystyle=\int_{M_{\rm min}}^{\infty}{\rm
d}n(M,z)\left({M\,W_{k}(M,z)\over\rho_{M0}}\right)^{2}\,,$ (3)

where $T(k)$ is the transfer function, $D(z)$ the growth function, $H_{0}$ is
the present-day Hubble parameter, $n_{s}$ is the spectral index, and
$\delta_{H0}$ is the amplitude of perturbations on the horizon scale today,
which we fix by setting $\sigma_{8}$. Furthermore, in Eq. (3), $M_{\rm min}$
is the smallest halo mass that is considered, ${\rm d}n$ is the number density
of halos in the mass range ${\rm d}M$, which is defined via the halo mass
function $f(M,z)$ that gives the fraction of the total mass in halos of mass
$M$ at the redshift $z$:

${\rm d}n(M,z)\equiv n(M,z){\rm d}M={\rho_{M0}\over M}\,f(M,z){\rm d}M\,,$ (4)

where $n(M,z)$ is the number density and $\rho_{M0}$ the matter density today.
Finally, $W_{k}$ is the Fourier transform of the halo density profile:

$W_{k}(M,z)=\frac{1}{M}\int_{0}^{R}\rho(r,M,z){\sin kr\over
kr}4\pi\,r^{2}\,{\rm d}r\,,$ (5)

and $R$ is the halo radius. For the halo profile $\rho(r,M,z)$, we use the
Navarro-Frenk-White (NFW) profile (Navarro et al., 1996). Modifications to the
inner cuspy NFW profile by baryonic feedback are neglected.

####  2.1.2 Lensing convergence

The lens convergence $\kappa$ in the Born approximation is given by the
following integral evaluated along the unperturbed light path (Bartelmann and
Schneider, 2001):

$\kappa(z_{s})=\int_{0}^{r_{s}}\textrm{d}r\,\rho_{M0}\,G(r,r_{s})\,\delta_{M}\big{(}r,t(r)\big{)}\,,$
(6)

where $G$ is the lensing efficiency of inhomogeneity at the comoving radius
$r$ (basically the inverse of the critical surface density defined below):

$G(r,r_{s})=\frac{4\pi
G}{c^{2}a}\,\frac{f_{k}(r)f_{k}(r_{s}-r)}{f_{k}(r_{s})}\,.$ (7)

The functions $a(t)$ and $t(r)$ are the scale factor and geodesic time for the
background FLRW model, $r_{s}=r(z_{s})$ is the comoving position of the source
at redshift $z_{s}$ and
$f_{k}(r)=\sin(r\sqrt{k})/\sqrt{k},\;r,\;\sinh(r\sqrt{-k})/\sqrt{-k}$
depending if the curvature $k>,=,<0$, respectively.

Because in Eq. (6) the contributions to the total convergence are combined
additively, it is useful to decompose the density field into a sum of Fourier
modes:

$\delta_{M}(\bm{r},t)=\int_{0}^{\infty}{d^{3}\bm{k}\over(2\pi)^{3}}\,e^{i\bm{k}\cdot\bm{r}}\,\delta_{M}(\bm{k},t)\,.$
(8)

Eq. (8) can be used to separate the contributions to the convergence due to
small-scale inhomogeneities from the ones due to large-scale inhomogeneities
so that, similarly to Eq. (1), it is:

$\kappa(z)=\kappa_{L}(z)+\kappa_{H}(z)\,.$ (9)

####  2.1.3 2-halo convergence PDF

Given Eqs. (2), (6) and (8), it is straightforward to obtain the variance of
the convergence PDF:

$\sigma^{2}_{\kappa_{L}}=\int_{0}^{r_{s}}\textrm{d}r\rho_{M0}^{2}G^{2}(r,r_{s})\int_{0}^{k_{L}}{k\,\textrm{d}k\over
2\pi}P_{L}\big{(}k,z(r)\big{)}\,.$ (10)

The third moment (skewness) can be obtained via the bispectrum (three-point
correlation function), which should also contribute to the variance (although
the contribution is expected to be small at these linear scales). Here, we
take an alternative approach and model the full PDF relative to the 2-halo
term as a log-normal distribution with zero mean and the variance given by Eq.
(10), see Marra et al. (2013) for details. This is motivated by the well-known
fact that the lensing PDF inherits the log-normal character of the density
contrast PDF.

In Eq. (10), the cut-off scale $k_{L}$ separates the 2-halo contribution,
characteristic of larger scales, from the 1-halo term, relevant at smaller
scales. The precise determination of $k_{L}$ is theoretically challenging,
especially because, at a given redshift, the dependence on the halo masses is
non-trivial. Consequently, we employ a pragmatic approach to establish this
scale, which involves adjusting it so that the lensing PDF derived from tgl
matches that obtained from PINOCCHIOhalos. A detailed description of the
method involving PINOCCHIOhalos is provided in Section 2.3. The latter’s
matter field representation closely matches indeed the one used by tgl, making
it a suitable reference for this calibration. We determined that
$k_{L}=\exp(3.9-4.6z)\text{ Mpc}^{-1}$ achieves a good agreement, showing that
$k_{L}$ increases at lower redshifts. This aligns with the fact that the halo
model does not represent low-density filamentary structures, which become
relatively more important at lower redshifts. Addressing this discrepancy
necessitates the inclusion of higher mode wavenumbers in calculating the
2-halo variance, thereby compensating for the model’s lack of power in these
regions.

####  2.1.4 1-halo convergence PDF

The contribution to the convergence from the 1-halo term could be similarly
computed via the halo power spectrum:

$\sigma^{2}_{\kappa_{H}}=\int_{0}^{r_{s}}\textrm{d}r\rho_{M0}^{2}G^{2}(r,r_{s})\int_{k_{H}}^{k_{\rm
max}}{k\,dk\over 2\pi}P_{H}\big{(}k,z(r)\big{)}\,,$ (11)

where $k_{\rm max}$ is the cut-off to account for the finite angular
resolution. The other cut-off scale $k_{H}$ is instead relative to the size of
the halos: for $k\ll k_{H}$, halos contribute to Eq. (11) only via shot noise.

Thanks to the sGL method, we can improve upon the previous result and obtain
an estimate equivalent to including the information from all the $n$-point
correlation functions. As explained in Kainulainen and Marra (2011a), the halo
contribution to Eq. (6) can be rewritten as:

$\displaystyle\kappa_{H}=\sum_{i,u}\,\kappa_{1iu}(z_{s})\Big{(}k_{iu}-\Delta
N_{iu}\Big{)}\,,$ (12)

where $\kappa_{1iu}$ is the convergence due to one object in the bin $iu$:

$\kappa_{1iu}(z_{s})\equiv G_{i}(r_{s})\,\Sigma_{iu}(t_{i})\,,$ (13)

where $\Sigma_{iu}$ is the surface mass density of the halo in the bin $iu$.
Here, the index $i$ labels the lens planes, while the index $u$ labels the
parameters that define the surface density; for halos these are the impact
parameter and the halo mass. Now, as far as the 1-halo term is concerned, the
integers $k_{iu}$ are distributed as Poisson random variables:

$P_{k_{iu}}=\frac{(\Delta N_{iu})^{k_{iu}}}{k_{iu}!}\,e^{-\Delta N_{iu}},$
(14)

where the parameter $\Delta N_{iu}$ is the expected number of objects in the
bin volume $\Delta V_{iu}$:

$\Delta N_{iu}=\Delta n_{iu}\,\Delta V_{iu}=\Delta n_{iu}\,\Delta
r_{i}\,\Delta A_{iu}\,,$ (15)

where $\Delta n_{iu}$ is the comoving density of objects corresponding to the
parameters in the bin $iu$, $\Delta r_{i}$ is the bin in the geodesic path,
and $\Delta A_{iu}$ is the corresponding cross-sectional area of the halos in
co-moving units.

From Eq. (12), using the properties of the Poissonian distribution, one can
see that $\langle\kappa_{H}\rangle=0$ and (Kainulainen and Marra, 2011a, b):

$\displaystyle\sigma^{2}_{\kappa_{H}}$
$\displaystyle=\sum_{i,u}\,\kappa_{1iu}^{2}\,\Delta N_{iu}\,,$ (16)
$\displaystyle\mu_{3,\kappa_{H}}$
$\displaystyle=\sum_{i,u}\,\kappa_{1iu}^{3}\,\Delta N_{iu}\,,$ (17)

from which, going back to integral form, one obtains:

$\displaystyle\sigma^{2}_{\kappa_{H}}$ $\displaystyle={1\over
N_{O}}\int_{0}^{r_{s}}\textrm{d}r\,G^{2}(r,r_{s})\int_{M_{\rm
min}}^{\infty}\textrm{d}n(M,z(r))$
$\displaystyle\times\int_{0}^{R(M,z(r))}\textrm{d}A(b,M)\,\Sigma^{2}(b,M,z(r)))\,,$
(18) $\displaystyle\mu_{3,\kappa_{H}}$ $\displaystyle={1\over
N_{O}}\int_{0}^{r_{s}}\textrm{d}r\,G^{2}(r,r_{s})\int_{M_{\rm
min}}^{\infty}\textrm{d}n(M,z(r))$
$\displaystyle\times\int_{0}^{R(M,z(r))}\textrm{d}A(b,M)\,\Sigma^{3}\big{(}b,M,z(r)\big{)}\,.$
(19)

The integral limits for the last two integrals are implicitly defined, ${\rm
d}A(b)\equiv 2\pi b\,{\rm d}b$ and $\Sigma$ is the halo surface density:

$\Sigma(b,M,z)=a^{3}\int_{b}^{R}\frac{2\,r\,{\rm
d}r}{\sqrt{r^{2}-b^{2}}}\,\rho(r,M,z)\,.$ (20)

Again, we want to point out that these expressions use the halo profiles in
real space and not in Fourier space, thus including higher-order correlation
terms beyond the power spectrum.

Obtaining the moments beyond $\mu_{3}$ analytically is more challenging as
there is no closed form. It is easier to estimate them by drawing a set of
Poissonian numbers $k_{iu}$, from which the full convergence PDF and its
moments can be easily estimated.

####  2.1.5 Full convergence PDF

The full convergence PDF is finally obtained by convolving the 1- and 2-halo
PDFs. The former is obtained via a histogram of a (Poissonian) realization;
the latter is based on the log-normal template previously discussed. The final
PDF is in the image plane.

The 1-halo contribution is based on binning the comoving space around the
photon geodesic. We adopt a fine enough binning such that the numerical
results converge. Furthermore, we smooth the NFW density profile via a
Gaussian kernel in order to have a fixed angular resolution. The PDF will then
be relative to this corresponding smoothing angular resolution, as is the case
for the other methods adopted in this study. Finally, the modeling of the
halos follows the prescription adopted in PINOCCHIO; that is, we adopt the
mass function by Watson et al. (2013, Eq. 12, FOF) and the concentration
parameter model by Bhattacharya et al. (2013, $200_{c}$). Furthermore, the
smallest halo we model ($M_{\rm min}$ in the previous equations) has a mass
that corresponds to the minimum halo mass adopted in the PINOCCHIO
simulations, that is, $2.6\times 10^{9}h^{-1}M_{\odot}$ (10 dark matter
particles, see Table 2).

###  2.2 Light-cone construction and ray tracing

The past light cone (PLC) lens maps for the PINOCCHIO and $N$-body simulations
were constructed using the SLICER999https://github.com/TiagoBsCastro/SLICER
code (Castro, 2024), introduced in Castro et al. (2018).101010SLICER is an
improved and fully rewritten version of MapSim (Giocoli et al., 2015). This
code distributes particles within the light cone based on the simulation
snapshots. The process begins with SLICER configuring the field of view and
setting the maximum source redshift, as specified in the input file.
Subsequently, the code determines the required number of lens planes to ensure
a continuous construction of the light cone. Each simulation snapshot is then
processed sequentially. The code reads the files, selecting particles within
the predetermined field of view. Prior to projecting these particles onto the
lens plane, each snapshot undergoes a randomization process. This step
involves reflecting and translating the center of the particle distribution,
in accordance with the periodic boundary conditions. Additionally, a specific
face of the simulation cube is chosen to align with the line of sight. The
final stage involves constructing the lens planes. This is achieved by mapping
the positions of particles to the nearest preselected lens plane. The angular
positions of the particles are maintained during this process. The surface
density is then pixelated using the triangular cloud method, as described in
Hockney and Eastwood (1988) and Bartelmann (2003). The grid pixels are
carefully selected to maintain a consistent angular size across all planes. To
prevent repetition of cosmic structures in the field of view (FoV) and to
ensure that the map sizes are equivalent to the simulation box size (set at
$150\textrm{ Mpc}/h$), all maps are constructed with a uniform FoV value of
$1^{\circ}$, see Table 2. Consequently, using SLICER, we successfully
generated 20 light cones up to $z_{s}=5$. These cones comprise 118 lens maps
with varying resolutions: $2048^{2}$, $1024^{2}$, $512^{2}$, and $256^{2}$
pixels. These resolutions correspond to angular resolutions of $1.76$, $3.52$,
$7.03$, and $14.06\textrm{ arcsec}$, respectively.

Table 2: Basic properties from the simulations used in this work. From left to right: the size of the box, gravitational softening, number of DM particles and their masses, the number of lens planes built up to $z=5$, the field of view of the constructed past light-cone, and the number of PLC realizations (different random seeds). $L_{\textrm{box}}$ | $\epsilon_{\rm{soften.}}$ | $N_{\textrm{part.}}$ | $m_{\textrm{DM}}$ | $N_{\rm planes}$ | FoV | PLCs
---|---|---|---|---|---|---
(Mpc/h) | (kpc/h) |  | $(M_{\odot}/h)$ | ($z\leq 5$) | (deg.) |
$150$ | $1.4$ | $1024^{3}$ | $2.6\times 10^{8}$ | $118$ | $1.0$ | 20

After creating the lens planes, the lensing quantities (convergence $\kappa$,
shear $\gamma$, and magnification $\mu$) are computed using the ray-tracing
technique, which integrates the non-linear matter on the unperturbed
trajectory of the source plane to the observer. Castro et al. (2018) showed
that the latter Born approximation yields accurate results for 1-point lensing
statistics. In the thin lens approximation, i.e., disregarding the size of
each lens, the convergence for a single lens plane is

$\kappa\left(\bm{\theta}\right)=\frac{\Sigma\left(\bm{\theta}\right)}{\Sigma_{\rm{crit}}},$
(21)

where $\bm{\theta}$ is the observed angular position,
$\Sigma\left(\bm{\theta}\right)$ is the surface mass density and
$\Sigma_{\rm{crit}}$ is the critical surface density

$\Sigma_{\rm{crit}}=\frac{c^{2}D_{l}}{4\pi GD_{s}D_{ls}}.$ (22)

In the above, c is the speed of light, G is the gravitational constant,
$D_{l}$, $D_{s}$ and $D_{ls}$ are the angular diameter distances between the
observer-lens, observer-source, and lens-source, respectively. Thus, we use
equation (21) to obtain the differential convergence maps, which are summed up
to the source redshift.

The shear maps were computed using the Shear Map Reconstruction (SMR) PYTHON
script distributed with SLICER. The shear components $\gamma_{1},\gamma_{2}$
are written as a function of the deflection potential $\Psi$:

$\displaystyle\gamma_{1}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{2}\left(\Psi_{\theta_{1}\theta_{1}}-\Psi_{\theta_{2}\theta_{2}}\right),$
(23) $\displaystyle\gamma_{2}\left(\bm{\theta}\right)$
$\displaystyle=\Psi_{\theta_{1}\theta_{2}}=\Psi_{\theta_{2}\theta_{1}},$ (24)

where

$\Psi\left(\bm{\theta}\right)=\frac{1}{\pi}\int\kappa\left(\bm{\theta^{\prime}}\right)\ln\left|\bm{\theta}-\bm{\theta^{\prime}}\right|{\rm
d}^{2}\theta^{\prime}\,.$ (25)

SMR uses FFTW to reconstruct the lensing potential from the convergence maps
and fourth-order finite differences to calculate its shear components.

Finally, the magnification $\mu$ is computed using the relation

$\mu\equiv\frac{1}{\left(1-\kappa\right)^{2}-\gamma^{2}},$ (26)

where $\gamma=\sqrt{\gamma_{1}^{2}+\gamma_{2}^{2}}$ is the shear modulus. In
the weak lensing regime, the above relation is reduced to

$\mu\approx 1+2\kappa+\gamma^{2}+3\kappa^{2}.$ (27)

We computed the angular power spectra $P_{\kappa}$ and $P_{\gamma}$ for all
the convergence and shear maps, respectively. In linear theory, it is
$P_{\kappa}=P_{\gamma}$, and for a source at $z_{s}$ it is written as an
integral of the density contrast along the unperturbed path of light, using
the Limber approximation (Limber, 1953; Munshi et al., 2008; Bartelmann and
Schneider, 2001):

$P_{\kappa}\left(l,z_{s}\right)\\!\\!=\\!\\!\frac{9H_{0}^{4}\Omega_{\rm
m0}^{2}}{4c^{4}}\\!\\!\int_{0}^{\chi_{s}}\\!\\!\\!\\!\\!\textrm{d}\chi^{\prime}\\!\left(\frac{\chi_{s}-\chi^{\prime}}{\chi_{s}\,a(\chi^{\prime})}\right)^{2}\\!\\!\\!P_{m}\left(\frac{l}{\chi^{\prime}},\chi^{\prime}\right),$
(28)

where $\chi(z)$ is the comoving distance at $z$, $\chi_{s}=\chi(z_{s})$, and
$P_{\rm m}$ is the linear matter power-spectrum.

###  2.3 PINOCCHIO

PINOCCHIO generates simulated halo catalogs using the ellipsoidal collapse to
compute the collapse time, defined as the moment of orbit crossing, and
Lagrangian perturbation theory (LPT) to displace the formed halos and
uncollapsed particles. Briefly, PINOCCHIO works as follows: a linear field is
smoothed using a Gaussian filter at various scales, the collapse time is
calculated for each particle (grid point) with the ellipsoidal collapse; these
are then sorted by collapse time in chronological order and displaced
according to the LPT theory; later, the code groups the collapsed particles
into halos and filaments in order to properly reproduce the hierarchical
formation of the structures; and finally the LPT theory is used to move the
halos to their final position and thus describe their evolution. As the
particles that are members of a halo have passed through orbit crossing, their
displacements are extrapolations of the LPT displacement outside its validity;
therefore, PINOCCHIO cannot provide reliable profiles needed for accurate
lensing statistics. The particles belonging to halos are then redistributed
assuming a NFW profile and the mass-concentration relation given by
Bhattacharya et al. (2013, $200_{c}$).

We use PINOCCHIO configured in the third order of the LPT theory to build a
set of simulations defined according to Table 2. In an effort to better
compare PINOCCHIO and sGL results, we consider two PINOCCHIO configurations
(see Table 1): i) the label PINOCCHIO represents the complete results
(halos+filaments+uncollapsed particles that do not belong to any structure) of
the standard PINOCCHIO code, and ii) the label $\texttt{\small
PINOCCHIO}_{\textrm{halos}}$ considers only the halos and randomizes the rest.

As a note, the construction of a lightcone from snapshots is an inefficient
use of PINOCCHIO, because the code itself outputs halos in a lightcone, and
the position of particles outside halos can be easily recovered at any time
using LPT. However, in this study, we prefer to apply to the PINOCCHIO runs
the same analysis applied to the $N$-body simulations.

Table 3: Ranges where the PDFs from the Hydro and DM versions of Magneticum Box 3 hr differ by less than 50% (1.76 arcsec). Statistics | Redshift | Range
---|---|---
Convergence | 1 | [$-\infty$, 0.38]
Convergence | 2 | $[-0.12,0.50]$
Convergence | 3 | $[-0.17,0.50]$
Convergence | 5 | $[-0.25,0.49]$
Shear | 1 | [$1.1\times 10^{-5}$, $5.8\times 10^{-2}$]
Shear | 2 | [$2.3\times 10^{-5}$, $1.3\times 10^{-1}$]
Shear | 3 | [$5.5\times 10^{-5}$, $1.5\times 10^{-1}$]
Shear | 5 | [$6.0\times 10^{-5}$, $1.9\times 10^{-1}$]
Magnification | 1 | $[0.93,2.3]$
Magnification | 2 | $[0.88,3.2]$
Magnification | 3 | $[0.83,2.4]$
Magnification | 5 | $[0.72,2.9]$
Magnitude | 1 | $[-0.90,0.08]$
Magnitude | 2 | $[-1.26,0.14]$
Magnitude | 3 | $[-0.95,0.20]$
Magnitude | 5 | $[-1.16,0.36]$

###  2.4 $N$-body simulations

The $N$-body simulations used in this work were performed using the GADGET-2
code (Springel, 2005), following the configurations described in Table 2. The
GADGET-2 code uses the parallelized Tree–Particle Mesh (TreePM, see Xu, 1995)
algorithm to calculate the gravitational interactions between the simulation
particles. In this algorithm, these interactions are separated into a long-
range term, calculated using the Particle-Mesh methods, and a short-range term
for interactions of close neighbors that are described via the Tree algorithm
(Barnes and Hut, 1986).

The simulations contain only collisionless matter accounting for the amount of
cold dark matter and baryons, excluding any effect of baryonic feedback. The
dark matter particles were evolved from $z=99$ to $z=0$ in a comoving box with
a side equal to 150 Mpc/h and mass resolution of $2.6\times
10^{8}\,M_{\odot}/h$ (see Table 2). The choice for this box size is justified
by the fact that it has already been shown in Castro et al. (2018) that a box
of 128 Mpc/$h$ with a mass resolution of $10^{8}\,M_{\odot}/h$ is enough to
obtain convergent results for the lensing PDF.

###  2.5 Hydrodynamical simulations

The effect of baryons on lensing probability density functions and the
convergence power spectrum $P_{\kappa}(\ell)$ was thoroughly examined by
Castro et al. (2018) using hydrodynamic simulations. This study specifically
analyzed the multipole range of $P_{\kappa}(\ell)$ and the value ranges for
$\kappa$, $\gamma$, and $\mu$ across various redshifts, identifying thresholds
below which the baryonic contributions are minimal. Given the challenges
associated with extracting cosmological information in scenarios where
baryonic effects—such as AGN and stellar feedback, star formation, magnetic
fields, and thermal conduction—are significant, our analysis adopts a
simplified approach. We thus avoid the complexities of expensive
hydrodynamical simulations by focusing on the dark matter-dominated regime,
ensuring that baryonic effects remain below the 50% threshold. Although more
stringent thresholds could be employed, it is pertinent to note that for many
lensing statistics, like the moments of the PDFs (Marra et al., 2013), the
discrepancies between DM-only and baryonic cases are minimal, as our DM model
effectively captures the bulk of the distributions.

Table 3 displays the ranges within which the PDFs from the Hydro and DM
versions of Box 3 (hr) of the Magneticum simulation
suite111111http://www.magneticum.org differ by less than 50%, given a
resolution of 1.76 arcsec. To avoid overcrowding in the plots, only the ranges
corresponding to redshifts $z=1$ and $z=3$ will be presented.

Figure 1: Convergence PDFs for the different lensing methods considered in
this study, for an angular resolution of 1.76 arcsec.

##  3 Results

We now present our results for the PDFs of convergence, shear, magnification,
and magnitude. In the ray-tracing method employed, rays are calculated from
the source to the observer, resulting in statistics computed in the image
plane. In most cases, we subsequently translate these PDFs to the source
plane, which is relevant to observations. Hereafter, the PDFs in the image and
source planes are denoted by the subscripts I and S, respectively.

As detailed in Castro et al. (2018), the lensing PDFs intrinsically depend on
the angular resolution. This dependency is not directly related to the point
spread function of any particular telescope, but rather to the angular scale
used for computing averages. For example, for a fixed lens and source
position, the average magnification of galaxies with diameters of 1 and 5
arcsec may differ. As discussed in Section 2.5, our simulations yield reliable
results for a minimal angular resolution of 1.76 arcsec, which we adopt as our
fiducial resolution for most discussions. However, results for lower
resolutions are presented in Appendix A.

###  3.1 Lensing PDFs

Figure 2: Comparison of $\kappa$, $\gamma$ and $\mu$ PDFs (in the source
plane) for different redshifts for the $N$-body and PINOCCHIO case for the
angular resolution of 1.76 arcsec. $\texttt{\small
PINOCCHIO}_{\textrm{halos}}$ (not shown here for compactness) exhibit similar
behavior to PINOCCHIO. The shaded regions mark where baryonic corrections
exceed 50% for $z=1$ (blue) or for $z=3$ (green). As can be seen the two
methods agree reasonably well in the range where baryonic effects are under
control. See also Fig. 3.

In Fig. 1, we display the $\kappa$ PDFs derived using the three different
methods considered in this study: full $N$-body simulations, PINOCCHIO
(including the ‘halos’ variant), and the stochastic turboGL code, as described
in Section 2 and summarized in Table 1. These are depicted for $z=1$ and $z=3$
at our fiducial angular resolution of 1.76 arcsec. The methods generally
exhibit good agreement, particularly in the peak position. However, the
$N$-body simulations show a slightly larger variance, as indicated by the
lower PDF peak and normalization condition. Since these plots probe weak
convergences and, therefore, mildly nonlinear structures, the slight
discrepancies should be attributed to a lack of modeling precision at the
scale of filaments, which will be discussed subsequently.

Figure 3: Relative residual differences between approximate methods and the
full $N$-body results in the source plane, for $z=1$ (solid lines) and $z=3$
(dashed lines). The blue (green) shaded regions represent the interval in
which the baryonic correction is over 50% for $z=1$ ($z=3$).

However, differences in the tails of the distributions are obscured by the
linear scale used in Fig. 1. Therefore, in Fig. 2, we present a general
comparison between the results from $N$-body and PINOCCHIO for $\kappa$,
$\gamma$, and $\mu$ on a logarithmic scale, depicting four different redshifts
for each method. The turboGL code is excluded from this comparison as it does
not compute the shear PDFs so that $\mu$ can be inferred only in the weak-
lensing regime. As previously discussed, our analysis focuses on the region
where baryonic effects are maintained below the 50% level. The regimes where
this is not the case at $z=1$ and $z=3$ are highlighted with a blue and green
shade, respectively. It is evident that as the source distance increases, the
convergence, shear, and magnification PDFs become wider (indicating higher
variance), and the peaks of their distributions shift (indicating skewness).
We will further explore these trends for convergence in Section 3.3, noting
the general trend that lensing intensifies at higher redshifts due to the
increased number of lens encounters. We find that both PINOCCHIO and $N$-body
simulations are consistent with the theoretical model of randomly distributed
point-mass lenses (Schneider et al., 1992), which predicts that when
$\gamma<0.01$, the PDF is proportional to $\gamma^{2}$, and for $\mu\gg 1$,
the PDF is proportional to $\mu^{-1}$. However, while the asymptotic behavior
prediction of the randomly distributed points model for the $\gamma-$PDF
happens within the range where the baryonic effects are subdominant, the
prediction for the $\mu-$PDF occurs in the regime deeply dominated by baryonic
physics and is not shown in Fig. 2. The behavior of the magnification PDF at
this regime will be discussed in Sect. 3.2.

Figure 3 illustrates the differences between the PDFs from PINOCCHIO and
turboGL compared to the $N$-body results. Both approximate methods exhibit a
small excess of probability density at the peak of the distribution and a
deficiency in the tail. These discrepancies stem from the limitations of the
approximate methods in accurately reproducing the non-linear density field,
likely due to an inadequate representation of substructures and halo
concentrations which are prevalent in this non-linear regime and responsible
for high-convergence events. We find that PINOCCHIO provides a better
approximation than the stochastic method at $z=1$, but both methods perform
comparably at $z=3$. However, both methods gradually fail for larger $\kappa$,
exhibiting approximately 20% discrepancies for $\kappa\geq 0.2$, which falls
well within the DM-dominated regime. Moreover, the lower the redshift, the
smaller the $\kappa$ value at which this breakdown occurs. Interestingly, for
the shear, we observe much better agreement between PINOCCHIO and $N$-body in
the region of interest. Differences become apparent at $\gamma\gtrsim 0.1$, a
region dominated by baryonic effects. An explanation for this disagreement
will be proposed in Section 3.5. For the $\mu-$PDFs, we again observe that as
we venture into more non-linear regimes, specifically for $\mu>2$, notable
differences emerge between both methods. However, these differences are
proportionally smaller in the DM-dominated regime than in the case of
$\kappa$. Particularly, the agreement is consistently good for $\mu<1$—where
the PDF rapidly decreases, causing relative residuals to increase—and ranges
between $20\%$ and $40\%$ for larger $\mu$ until baryonic physics becomes
dominant. As we will explore later, the multiple-image regime begins beyond
these baryonic thresholds; thus, only the so-called type-I $\mu$-PDF is
relevant here (see, e.g., Takahashi et al., 2011; Castro et al., 2018). A
similar comparison for different angular resolutions is presented in Appendix
A.

###  3.2 Strong lensing PDFs

Figure 4: Convergence, shear and magnification PDFs at $z=5$ for different
angular resolutions for the different image types.

Here, we present the lensing PDFs in the strong-lensing regime, computed by
differentiating contributions from various image types defined as follows:

Type I: $\displaystyle|1-\kappa|>\gamma\;\;\rm{and}\;\;\kappa<1\,,$ Type II:
$\displaystyle|1-\kappa|<\gamma\,,$ Type III:
$\displaystyle|1-\kappa|>\gamma\;\;\rm{and}\;\;\kappa>1\,.$ (29)

Type II and III images indicate the strong-lensing regime, characterized by
multiple images (Schneider et al., 1992), while Type I involves only image
amplification.

Figure 4 displays the PDFs of $\kappa$, $\gamma$, and $\mu$ for $z=5$,
categorized into these three types. We see that the $N$-body simulations show
a higher probability of forming Types II and III images compared to structures
generated by the PINOCCHIO code, highlighting the limitations of the method in
accurately describing the strong-lensing regime (in the present implementation
the halo internal structure is described by a single spherical NFW profile
without substructure). The absence of most Type II and III images in PINOCCHIO
is also consistent with observations from the previous section that the
PINOCCHIO PDFs disagree from the $N$-body results for $\gamma\gtrsim 0.1$.
Furthermore, we find that the angular resolution does not significantly impact
the lensing PDFs in the strong-lensing regime, and comparisons between the
PINOCCHIO and PINOCCHIOhalos PDFs show nearly indistinguishable results, with
the exception of the Type II $\kappa$ PDFs.

Additionally, we corroborate several findings from Castro et al. (2018):

  * •

The magnification PDFs from the $N$-body simulations indicate that the
observed plateau in the demagnification region results from Type II images.

  * •

The emergence of secondary peaks in the high shear region is due to the
transition from Type I to Types II and III.

  * •

The plateau in the high-convergence tail of the convergence PDFs, starting at
$\kappa\sim 1$, stems from the presence of Type III images.

###  3.3 Dissecting the PDF

Figure 5 shows the halo-model convergence PDF from turboGL together with the
individual contributions from the 1- and 2-halo terms. We can see that the
total lensing variance is dominated, as expected, by the 1-halo contribution
and that the log-normal template features a mild skewness
$\tilde{\mu}_{3}=\mu_{3}/\sigma^{3}$ (the third standardized moment).

Next, we show, breaking the sums of Eqs. (16) and (17) into individual
elements, the contributions to the variance $\mu_{2}$ and skewness
$\tilde{\mu}_{3}$ of halos with a given mass (Fig. 6) and in a given lensing
plane redshift (Fig. 7), for sources at $z=1$ and $z=3$. It is also shown the
average number of lensing-halo hits per light ray. In the case of Fig. 7, a
value of 10 in a given $\Delta z$, for instance, means that the light ray hit
on average 10 halos in that redshift bin (bin widths are constant in comoving
distance), and similarly for Fig. 6.

Figure 5:  Image-plane convergence PDFs relative to the 1-halo and 2-halo
contributions together with the full PDF obtained by convolving the previous
two. One can see that the 1-halo contribution is more significant than the
2-halo one.

Figure 6:  Contribution to the variance $\mu_{2}$ (top) and skewness
$\tilde{\mu}_{3}=\mu_{3}/\sigma^{3}$ (bottom) of lensing halos with a given
mass for sources at $z=1$ and $z=3$ using Eqs. (16)–(17). Also shown (point
markers) is the average number of lensing-halo hits per light ray. See Section
3.3.

From Fig. 6 (top panel) we can see that our simulations consider the relevant
range of masses as the halos with mass $\approx 10^{13}h^{-1}M_{\odot}$
contribute to most of the variance, and halos with masses below
$10^{9}h^{-1}M_{\odot}$ or above $10^{15}h^{-1}M_{\odot}$ have a negligible
impact. Also, the lensing PDF for sources at $z=3$ has a variance which is
about 5 times greater than the PDF for sources at $z=1$ – the light ray
encounters more lenses in the former case. Regarding skewness, we observe that
more massive halos contribute significantly to the overall asymmetry of the
PDF. This is due to less massive halos being more numerous, thereby
distributing the same fraction of the Universe’s mass across a greater number
of objects. Consequently, a ray experiences a higher number of interactions
(hits). According to the central limit theorem, this results in a PDF that
more closely approximates a Gaussian distribution. Finally, a source at $z=3$
is further away, and a larger number of halos are encountered so that the
number of hits is larger and the skewness is lower. At the very high-mass end,
the skewness for sources at $z=3$ is slighter higher because the extra number
of hits with respect to $z=1$ does not compensate for the very skewed
distributions of those very massive halos.

The top panel of Fig. 7 shows the contributions of the lens planes to the
variance. The shape is modulated by the lensing efficiency of Eq. (7) and the
growth of structures. The largest contributions come at about half way between
source and observer. The bottom panel shows the contribution to the skewness.
The overall trend is that lens planes at lower redshifts contribute to a
higher skewness. From Fig. 6, we remember that smaller halos have a smaller
number of hits at lower redshifts (and produce relatively more skewness),
while this is the opposite for larger halos. The overall effect is dominated
by the smaller halos, which can be attributed to the progressive nonlinearity
in the Universe due to the growth of structures. We also see that the
contributions for different source redshifts are exactly one on top of the
other. It is useful to remember that the skewness in the matter and lensing
distributions is due to the fact that under-densities occupy more volume than
overdensities, so it is more likely to encounter an under-density or have a
de-magnified light beam. Therefore, the third standardized moment is only due
to geometry and does not depend on the lensing efficiency of Eq. (7), which
depends on the source redshift.

Figure 7:  Contribution to the variance (top) and skewness (bottom) of halos
of a given lensing plane redshift for sources at $z=1$ and $z=3$ using Eqs.
(16)-(17). Also shown (point markers) is the average number of lensing-halo
hits per light ray. See Section 3.3.

###  3.4 Momenta of the magnitude PDF

We now focus on analyzing the second, third, and fourth central moments of the
shift in magnitude $\Delta m$, which directly relates to magnification $\mu$
through:

$\Delta m\equiv 2.5\log_{10}\mu.$ (30)

Using magnitudes is advantageous because the magnitude PDF exhibits more
stable central moments than the $\mu$–PDF, which has a divergent variance in
cosmology. Nevertheless, truncating the magnitude PDF at a certain negative
$\Delta m$ value is prudent, as highly magnified objects could, in principle,
be individually detected and excluded from the statistics either a posteriori
or during target selection a priori.121212These methods accurately recover the
mean of the non-truncated PDFs, which is essential for unbiased cosmological
analyses. Since the mean is a non-central moment, we will not discuss here the
effects of truncation on it. Furthermore, as previously demonstrated, the
baryonic contribution becomes dominant for $\mu\gtrsim 3$, which is outside
the focus of this study.

Results for redshifts 1 and 3, computed in the source plane with a fixed
angular resolution of 1.76 arcsec, are depicted in Fig. 8. These PDFs are
assessed up to a cutoff magnitude, $\Delta m_{\rm cut}$, beyond which the PDF
is considered zero. The second central moment shows the strongest agreement,
indicating PINOCCHIO’s capability to accurately reproduce the body rather than
the tail of the PDF distribution – a consistency noted in earlier sections.
Particularly, at $z=3$, PINOCCHIO and PINOCCHIOhalos agree more closely with
$N$-body results compared to $z=1$, except for $\mu_{2,\rm{lens}}$ in
PINOCCHIOhalos. We interpret the latter disagreement as due to the higher
relative importance of non-halo particles at higher redshifts, at which halos
give a less complete description of the matter field (see next section). This
indicates that non-halo particles have a significant impact on the central
moments, suggesting that the discrepancies with the $N$-body results are not
simply due to the halo profiles and positioning.

Overall, PINOCCHIO well models statistical quantities within limits where
baryonic effects are negligible. This precision makes it a useful tool for
constructing simple models of central moments as functions of redshift and
cosmological parameters. These models are crucial for the MeMo method (Quartin
et al., 2014), which utilizes the weak lensing effect on the intrinsic
distribution of magnitudes to constrain key cosmological parameters, such as
the amplitude of the matter power spectrum and the matter density parameter
(Quartin et al., 2014; Castro and Quartin, 2014; Castro et al., 2016a).

Figure 8: Central moments’ ratio of the magnitude PDFs of the approximative
methods with respect to the $N$-body simulation, for different values of
$\Delta m_{\rm cut}$ and for $z=1$ and $z=3$. The blue (green) shaded regions
indicate the limit where the baryonic effects is dominant, above 50%, for
$z=1$ ($z=3$).

Figure 9: Projected density maps for the angular resolution of 1.76 arcsec.
The projected slice has a thickness of 150 Mpc/$h$.

###  3.5 Maps and power spectrum

Finally, to gain a deeper understanding of our findings, we will analyze the
maps of the lensing quantities that underpin the PDFs discussed in previous
sections, along with the convergence power spectra generated by the PINOCCHIO
code and the $N$-body simulations across various redshifts.

In Fig. 9, we display the projected density maps of matter distributed over
the lens plane at redshifts $z=1$ and $z=3$, reconstructed using the $N$-body
simulations and PINOCCHIO in its two configurations. These maps highlight the
distribution of halos and filaments, as well as the presence of galaxy
clusters at the intersections of filaments. We observe a better agreement
between PINOCCHIO and $N$-body at $z=3$, with PINOCCHIOhalos showing weaker
structures. This is consistent with our previous observation that non-halo
particles become more significant at higher redshifts, where the simulated
halos provide a less comprehensive depiction of the matter field. Indeed, the
absence of halos smaller than those defined by the mass resolution of the
simulation leads to a relatively greater importance of diffuse matter between
halos, responsible for power on large-scale modes crucial for replicating the
linear distribution of structures. Conversely, while PINOCCHIO and
PINOCCHIOhalos agree well at $z=1$, their correspondence with the $N$-body
simulation worsens, as halo profiles and positions do not accurately replicate
the small-scale patterns observed in the $N$-body simulations. These results
also highlight the limitations of the LPT approximation up to the 3rd order,
employed in the PINOCCHIO code, in capturing the distribution of matter on
smaller scales. Notably, the absence of more massive structures, which are
evident in the $N$-body simulation maps, confirms this limitation, aligning
with findings from early works on the PINOCCHIO code (Monaco, 2016; Munari et
al., 2017).

The convergence maps and power spectra for the $N$-body and PINOCCHIO
simulations are displayed in Fig. 10. The PINOCCHIO maps begin to diverge from
those of the $N$-body simulations at small scales, yet show good agreement on
larger scales. This observation is corroborated both through visual inspection
and by analyzing the angular convergence power spectrum, which shows a
1–2$\sigma$ consistency for $l\lesssim 10^{4}$, consistent with the
observations from Fig. 9. At small scales, the contribution of the halo
profile becomes dominant, where the lack of substructures in the profile
adopted here implies the disagreement observed in this limit.131313See Giocoli
et al. (2017) for the effects of substructures in the power spectrum and PDFs.
Additionally, Fig. 11 presents a comparison of the shear maps. There is broad
agreement between the maps for both $z=1$ and $z=3$. Distinctions between
PINOCCHIO, PINOCCHIOhalos, and the $N$-body simulations become noticeable only
at high-shear levels. Here, the filaments represented in PINOCCHIO appear more
pronounced than those in the $N$-body simulation, highlighting the limits of
agreement. This finding is consistent with the analysis of the PDFs detailed
in Sections 3.1 and 3.2, which highlighted disagreement for $\gamma\gtrsim
0.1$. In particular, we can infer that the correct modeling of filaments is
important for type II and III images.

Figure 10: Comparison of convergence maps and angular power spectra for the
angular resolution of 1.76 arcsec. The displayed ratio represents the average
across 20 PLCs, with error bands reflecting the propagated uncertainties of
this average.

Figure 11: Shear maps for the angular resolution of 1.76 arcsec.

##  4 Conclusions

In this paper, we investigated the effectiveness of the approximate methods
PINOCCHIO and turboGL in accurately reproducing lensing probability density
functions (PDFs). Our simulations assume a universe composed solely of dark
matter, excluding baryonic effects. We studied the impact of gravitational
lensing in both weak and strong regimes, differentiating the contributions
from the three types of images. We thoroughly compared the results of these
approximate methods with those obtained from $N$-body simulations. Key
findings include:

  1. 1.

Approximate methods are effective primarily within the weak and medium-lensing
region, which mainly arises from mildly nonlinear matter structures.
Accurately modeling the tails of the PDFs necessitates $N$-body simulations,
as demonstrated by our analysis of Types I, II, and III lensing PDFs.

  2. 2.

For supernovae and standard sirens these methods prove useful for statistical
analyses, given that strong lensing events are either rare or can be excluded.
Within this regime, baryonic effects are expected to have a minor impact.

  3. 3.

These approximate methods struggle to accurately capture the small-scale
nonlinear matter field. As the perturbation theory fails to accurately predict
small-scale halo clustering and intertwining matter distribution, enhancing
the modeling of halo internal structures could improve results. Recent work by
Berner et al. (2022) has initiated efforts to reconstruct halo substructure
from merger histories, which could be integrated into these methods to add
more small-scale power to the simulations. We also observed that non-halo
particles minimally impact lensing for $z<3$.

  4. 4.

Our evaluation of the second-to-fourth moments of the lensing PDFs indicated
that PINOCCHIO agrees with $N$-body results to within 8–15% in regimes where
baryonic effects are minimal.

  5. 5.

The stochastic simulations of turboGL are suitable for rapid likelihood
exploration due to their speed and flexibility in handling low-statistics PDFs
($\sim 10^{5}$ rays). In contrast, ray tracing is more effective at exploring
high-$\mu$ tails as the number of rays increases with the square of the field
of view, while stochastic methods scale linearly with sample size.

  6. 6.

The analyses of the moments in Figs. 6 and 7 and the maps (Figs. 9–11)
underscore the need for precise modeling of mass distribution and non-linear
clustering of halos with $M>10^{12}M_{\odot}$. The challenges are twofold:
($i$) while there is a moderate consensus on the impact of baryons at the
scale of galaxy clusters, their effects in individual galaxies and galaxy
groups remain far less understood and more uncertainties (see, e.g., Euclid
Collaboration: Castro, T., et al., 2023); ($ii$) lighter objects in filaments
require both accurate modeling of their individual profiles and a detailed
understanding of their interconnected mass distributions.

These insights are expected to guide future studies on gravitational lensing
of point sources, a topic of growing relevance in light of upcoming supernova
and gravitational wave datasets.

## Acknowledgements

TC is supported by the Agenzia Spaziale Italiana (ASI) under - Euclid-FASE D
Attivita’ scientifica per la missione - Accordo attuativo ASI-INAF n.
2018-23-HH.0, by the PRIN 2022 PNRR project "Space-based cosmology with
Euclid: the role of High-Performance Computing" (code no. P202259YAF), by the
Italian Research Center on High-Performance Computing Big Data and Quantum
Computing (ICSC), project funded by European Union - NextGenerationEU - and
National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2, within
the activities of Spoke 3, Astrophysics and Cosmos Observations, by the INFN
INDARK PD51 grant, and by the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA. VM
is supported by the Brazilian research agencies FAPES and Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq). MQ is supported by the
research agencies Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado
do Rio de Janeiro (FAPERJ) project E-26/201.237/2022 and CNPq. We acknowledge
the use of the Santos Dumont supercomputer of the National Laboratory of
Scientific Computing (LNCC, Brazil) through the project lenssims.

## Data and Software Availability

Codes will be shared on reasonable request to the corresponding author. Data
products and post-processing routines are available at (Castro et al.,
2024).141414github.com/TiagoBsCastro/LensingMethodsDeconstruction

## References

  * Abbott et al. (2018) Abbott T. M. C., et al. (2018) Dark Energy Survey year 1 results: Cosmological constraints from galaxy clustering and weak lensing. _Phys. Rev. D_ , 98(4):043526, 10.1103/PhysRevD.98.043526, 1708.01530
  * Amendola et al. (2015) Amendola L., Castro T., Marra V., Quartin M. (2015) Constraining the growth of perturbations with lensing of supernovae. _MNRAS_ , 449:2845–2852, 10.1093/mnras/stv497, URL https://academic.oup.com/mnras/article/449/3/2845/2892989, 1412.3703
  * Asgari et al. (2023) Asgari M., Mead A. J., Heymans C. (2023) The halo model for cosmology: a pedagogical review. 10.21105/astro.2303.08752, 2303.08752
  * Barnes and Hut (1986) Barnes J., Hut P. (1986) A Hierarchical 0 Nlogn Force Calculation Algorithm. _Nature_ , 324:446–449, 10.1038/324446a0
  * Bartelmann (2003) Bartelmann M. (2003) Numerical methods in gravitational lensing. _ArXiv e-prints_ , astro-ph/0304162
  * Bartelmann and Schneider (2001) Bartelmann M., Schneider P. (2001) Weak gravitational lensing. _Phys. Rep._ , 340:291–472, 10.1016/S0370-1573(00)00082-X, URL https://doi.org/10.1016/S0370-1573(00)00082-X, 9912508
  * Berner et al. (2022) Berner P., Refregier A., Sgier R., Kacprzak T., Tortorelli L., Monaco P. (2022) Rapid simulations of halo and subhalo clustering. _JCAP_ , 11:002, 10.1088/1475-7516/2022/11/002, 2112.08389
  * Bhattacharya et al. (2013) Bhattacharya S., Habib S., Heitmann K., Vikhlinin A. (2013) Dark Matter Halo Profiles of Massive Clusters: Theory versus Observations. _ApJ_ , 766(1):32, 10.1088/0004-637X/766/1/32, 1112.5479
  * Canevarolo and Chisari (2023) Canevarolo S., Chisari N. E. (2023) Lensing bias on cosmological parameters from bright standard sirens. 2310.12764
  * Canevarolo et al. (2024) Canevarolo S., van Vonderen L., Chisari N. E. (2024) Impact of lensing of gravitational waves on the observed distribution of neutron star masses. 2404.11480
  * Castro (2024) Castro T. (2024) Tiagobscastro/slicer: v2.0.1. 10.5281/zenodo.11048430, URL https://doi.org/10.5281/zenodo.11048430
  * Castro and Quartin (2014) Castro T., Quartin M. (2014) First measurement of $\sigma_{8}$ using supernova magnitudes only. _Mon. Not. Roy. Astron. Soc._ , 443:L6, 10.1093/mnrasl/slu071, 1403.0293
  * Castro et al. (2016a) Castro T., Marra V., Quartin M. (2016a) Constraining the halo mass function with observations. _Mon.Not.Roy.Astron.Soc._ , 463:1666–1677, 10.1093/mnras/stw2072, 1605.07548
  * Castro et al. (2016b) Castro T., Quartin M., Benitez-Herrera S. (2016b) Turning noise into signal: learning from the scatter in the Hubble diagram. _Phys. Dark Univ._ , 13:66–76, 10.1016/j.dark.2016.04.006, 1511.08695
  * Castro et al. (2018) Castro T., Quartin M., Giocoli C., Borgani S., Dolag K. (2018) The effect of baryons in the cosmological lensing PDFs. _MNRAS_ , 478:1305–1325, 10.1093/mnras/sty1117, URL https://doi.org/10.1093/mnras/sty1117, 1711.10017
  * Castro et al. (2024) Castro T., Marra V., Quartin M., Alfradique V. (2024) Tiagobscastro/lensingmethodsdeconstruction: v1.0.1. 10.5281/zenodo.11093020, URL https://doi.org/10.5281/zenodo.11093020
  * Euclid Collaboration: Castro, T., et al. (2023) Euclid Collaboration: Castro, T., et al. (2023) Euclid preparation. XXXIX. The effect of baryons on the Halo Mass Function. 2311.01465
  * Giocoli et al. (2015) Giocoli C., Metcalf R. B., Baldi M., Meneghetti M., Moscardini L., Petkova M. (2015) Disentangling dark sector models using weak lensing statistics. _Mon. Not. Roy. Astron. Soc._ , 452(3):2757–2772, 10.1093/mnras/stv1473, 1502.03442
  * Giocoli et al. (2017) Giocoli C., Di Meo S., Meneghetti M., Jullo E., de la Torre S., Moscardini L., Baldi M., Mazzotta P., et al. (2017) Fast weak-lensing simulations with halo model. _Mon. Not. Roy. Astron. Soc._ , 470(3):3574–3590, 10.1093/mnras/stx1399, 1701.02739
  * Hilbert et al. (2020) Hilbert S., Barreira A., Fabbian G., Fosalba P., Giocoli C., Bose S., Calabrese M., Carbone C., et al. (2020) The accuracy of weak lensing simulations. _MNRAS_ , 493:305–319, 10.1093/mnras/staa281, URL https://doi.org/10.1093/mnras/staa281, 1910.10625
  * Hockney and Eastwood (1988) Hockney R. W., Eastwood J. W. (1988) Computer simulation using particles. CRC Press
  * Kainulainen and Marra (2009) Kainulainen K., Marra V. (2009) New stochastic approach to cumulative weak lensing. _Phys. Rev. D_ , 80:123020, 10.1103/PhysRevD.80.123020", URL https://doi.org/10.1103/PhysRevD.80.123020, 0909.0822
  * Kainulainen and Marra (2011a) Kainulainen K., Marra V. (2011a) Accurate modeling of weak lensing with the stochastic gravitational lensing method. _Phys. Rev. D_ , 83:023009, 10.1103/PhysRevD.83.023009, URL https://doi.org/10.1103/PhysRevD.83.023009, 1011.0732
  * Kainulainen and Marra (2011b) Kainulainen K., Marra V. (2011b) Weak lensing observables in the halo model. _Phys. Rev. D_ , 84:063004, 10.1103/PhysRevD.84.063004, URL https://doi.org/10.1103/PhysRevD.84.063004, 1101.4143
  * Kilbinger (2015) Kilbinger M. (2015) Cosmology with cosmic shear observations: a review. _Rept. Prog. Phys._ , 78:086901, 10.1088/0034-4885/78/8/086901, 1411.0115
  * Kitaura and S. (2013) Kitaura F.-S., S. H. (2013) Cosmological structure formation with augmented Lagrangian perturbation theory. _Mon. Not. R. Astron. Soc.: Lett._ , 435:L78–L82, 10.1093/mnrasl/slt101, URL https://doi.org/10.1093/mnrasl/slt101, 1212.3514
  * Limber (1953) Limber D. N. (1953) The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field. _ApJ_ , 117:134, 10.1086/145672, URL https://articles.adsabs.harvard.edu/pdf/1953ApJ...117..134L
  * Macaulay et al. (2017) Macaulay E., Davis T. M., Scovacricchi D., Bacon D., Collett T. E., Nichol R. C. (2017) The effects of velocities and lensing on moments of the Hubble diagram. _Mon. Not. Roy. Astron. Soc._ , 467(1):259–272, 10.1093/mnras/stw3339, 1607.03966
  * Macaulay et al. (2020) Macaulay E., et al. (2020) Weak Lensing of Type Ia Supernovae from the Dark Energy Survey. _Mon. Not. Roy. Astron. Soc._ , 496(3):4051–4059, 10.1093/mnras/staa1852, 2007.07956
  * Manera et al. (2012) Manera M., Scoccimarro R., Percival W. J., Samushia L., McBride C. K., Ross A. J., Sheth R. K., White M., et al. (2012) The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: a large sample of mock galaxy catalogues. _MNRAS_ , 428:1036–1054, 10.1093/mnras/sts084, URL https://doi.org/10.1093/mnras/sts084, 1203.6609
  * Marra et al. (2013) Marra V., Quartin M., Amendola L. (2013) Accurate weak lensing of standard candles. I. Flexible cosmological fits. _Phys. Rev. D_ , 88:063004, 10.1103/PhysRevD.88.063004, 1304.7689
  * Monaco (2016) Monaco P. (2016) Approximate methods for the generation of dark matter halo catalogs in the age of precision cosmology. _Galaxies_ , 4(4):53, 10.3390/galaxies4040053, 1605.07752
  * Monaco et al. (2002) Monaco P., Theuns T., G. T. (2002) The PINOCCHIO algorithm: pinpointing orbit-crossing collapsed hierarchical objects in a linear density field. _MNRAS_ , 331:587–608, 10.1046/j.1365-8711.2002.05162.x, URL https://doi.org/10.1046/j.1365-8711.2002.05162.x, 0109323
  * Mpetha et al. (2024) Mpetha C. T., Congedo G., Taylor A., Hendry M. A. (2024) Impact of weak lensing on bright standard siren analyses. _ArXiv e-prints_ , 10.48550/arXiv.2402.19476, URL https://arxiv.org/abs/2402.19476, 2402.19476
  * Munari et al. (2016) Munari E., Monaco P., Sefusatti E., Castorina E., Mohammad F. G., Anselmi S., Borgani S. (2016) Improving fast generation of halo catalogues with higher order Lagrangian perturbation theory. _MNRAS_ , 465:4658–4677, 10.1093/mnras/stw3085, URL https://doi.org/10.1093/mnras/stw3085, 1605.04788
  * Munari et al. (2017) Munari E., Monaco P., Koda J., Kitaura F.-S., Sefusatti E., Borgani S. (2017) Testing approximate predictions of displacements of cosmological dark matter halos. _JCAP_ , 07:050, 10.1088/1475-7516/2017/07/050, URL https://doi.org/10.1088/1475-7516/2017/07/050, 1704.00920
  * Munshi et al. (2008) Munshi D., Valageas P., van Waerbeke L., Heavens A. (2008) Cosmology with weak lensing surveys. _Phys. Rep._ , 462:67–121, 10.1016/j.physrep.2008.02.003, URL https://doi.org/10.1016/j.physrep.2008.02.003, 0612667
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M. (1996) The Structure of cold dark matter halos. _Astrophys. J._ , 462:563–575, 10.1086/177173, astro-ph/9508025
  * Peacock and Smith (2000) Peacock J. A., Smith R. E. (2000) Halo occupation numbers and galaxy bias. _Mon. Not. Roy. Astron. Soc._ , 318:1144, 10.1046/j.1365-8711.2000.03779.x, astro-ph/0005010
  * Quartin et al. (2014) Quartin M., Marra V., Amendola L. (2014) Accurate weak lensing of standard candles. II. Measuring $\sigma_{8}$ with supernovae. _Phys. Rev. D_ , 89:023009, 10.1103/PhysRevD.89.023009, URL https://journals.aps.org/prd/abstract/10.1103/PhysRevD.89.023009, 1307.1155
  * Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E. (1992) Gravitational Lenses. Springer Berlin, 10.1007/978-3-662-03758-4
  * Scovacricchi et al. (2017) Scovacricchi D., Nichol R. C., Macaulay E., Bacon D. (2017) Measuring weak lensing correlations of Type Ia Supernovae. _Mon. Not. Roy. Astron. Soc._ , 465(3):2862–2872, 10.1093/mnras/stw2878, 1611.01315
  * Shan et al. (2021) Shan X., Wei C., Hu B. (2021) Lensing magnification: gravitational waves from coalescing stellar-mass binary black holes . _MNRAS_ , 508:1253–1261, 10.1093/mnras/stab2567, URL https://doi.org/10.1093/mnras/stab2567
  * Shapiro et al. (2010) Shapiro C., Bacon D. J., Hendry M., Hoyle B. (2010) Delensing gravitational wave standard sirens with shear and flexion maps. _MNRAS_ , 404:858–866, 10.1111/j.1365-2966.2010.16317.x, URL https://doi.org/10.1111/j.1365-2966.2010.16317.x
  * Springel (2005) Springel V. (2005) The cosmological simulation code GADGET-2. _MNRAS_ , 364:1105–1134, 10.1111/j.1365-2966.2005.09655.x, URL https://doi.org/10.1111/j.1365-2966.2005.09655.x, 0505010
  * Taffoni et al. (2002) Taffoni G., Monaco P., Theuns T. (2002) PINOCCHIO and the hierarchical build-up of dark matter haloes. _MNRAS_ , 333:623–632, 10.1046/j.1365-8711.2002.05441.x, URL https://doi.org/10.1046/j.1365-8711.2002.05441.x, 0109324
  * Takahashi et al. (2011) Takahashi R., Oguri M., Sato M., Hamana T. (2011) Probability distribution functions of cosmological lensing: convergence, shear, and magnification. _ApJ_ , 742:15, 10.1088/0004-637X/742/1/15, URL https://doi.org/10.1088/0004-637X/742/1/15, 1106.3823
  * Tassev et al. (2013) Tassev S., Zaldarriaga M., Eisenstein D. (2013) Solving large scale structure in ten easy steps with COLA. _JCAP_ , 06:036, 10.1088/1475-7516/2013/06/036, URL https://iopscience.iop.org/article/10.1088/1475-7516/2013/06/036, 1301.0322
  * Thiele et al. (2020) Thiele L., Hill J. C., Smith K. M. (2020) Accurate analytic model for the weak lensing convergence one-point probability distribution function and its autocovariance. _Phys. Rev. D_ , 102(12):123545, 10.1103/PhysRevD.102.123545, 2009.06547
  * Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G. (2013) The halo mass function through the cosmic ages. _MNRAS_ , 433(2):1230–1245, 10.1093/mnras/stt791, 1212.0095
  * Wu et al. (2023) Wu Z. F., Chan L. W. L., Hendry M., Hannuksela O. A. (2023) Reducing the impact of weak-lensing errors on gravitational-wave standard sirens. _MNRAS_ , 522:4059–4077, 10.1093/mnras/stad1194, URL https://doi.org/10.1093/mnras/stad1194
  * Xu (1995) Xu G.-H. (1995) A new parallel N body gravity solver: TPM. _Astrophys. J. Suppl._ , 98:355, 10.1086/192166, astro-ph/9409021

Figure 12: Similar to Figure 2, but for different angular resolutions at
$z=3$. We use the 1.76" resolution as reference for the shaded regions.

##  Appendix A PDFs at different angular resolutions

Figure 12 illustrates how changes in angular resolution affect the lensing
PDFs from PINOCCHIO and $N$-body simulations. Increasing the angular
resolution angle alters the peak and variance of the lensing PDFs in the
following ways:

  1. i)

the convergence peak shifts to lower $|\kappa|$ values and the variance
decreases;

  2. ii)

the shear peak shifts to lower $\gamma$ values without a noticeable change in
its variance;

  3. iii)

the magnification exhibits a combined effect of i) and ii), where the peak
shifts towards $\mu=1$ and the variance decreases.

This behavior is anticipated as using a larger resolution angle effectively
smooths the matter field, thus reducing the overall lensing effect. Notably,
the characteristic elbow at $\gamma\gtrsim 0.1$ vanishes with a larger
resolution angle. As discussed in Section 3.2, this feature is linked to type
II and III images, which are suppressed by the effective smoothing of the
matter field.