© The Authors 2025

1. Introduction

Despite the success of the Λ-cold dark matter (ΛCDM) concordance model (Heavens et al. 2017), which accurately describes observations from well-tested Solar System dynamics (Will 2014; Baker et al. 2021) to gravitational wave propagation (Creminelli & Vernizzi 2017; Abbott et al. 2017), some tension remains among the observations. In particular, there are discrepancies between cosmic microwave background (CMB) (Planck Collaboration VI 2020; Planck Collaboration V 2020) constraints and the lensing measurements of σ₈ (Abdalla et al. 2022; Di Valentino et al. 2021a) as well as the supernovae determinations of H₀ (Di Valentino et al. 2021b). These discrepancies have motivated tests of gravity on cosmic scales (Koyama 2016).

In response, modified gravity (MG) theories have been proposed to address observational tensions (Carroll 2001; Martin 2012; Moresco & Marulli 2017), including the possibility that general relativity (GR) fails on cosmological scales (Dolgov & Kawasaki 2003; Clifton et al. 2012; Joyce et al. 2015; Ishak 2019). These models introduce additional degrees of freedom, modifying structure formation and spacetime geometry. Among them, f(R) gravity mimics dark energy (DE) effects via an additional scalar field while employing screening mechanisms to restore GR predictions at small scales and in high-density regions (Joyce et al. 2015; Bertotti et al. 2003; Will 2005; Hinterbichler & Khoury 2010; Brax & Valageas 2013). Testing MG models is challenging due to their similarity to ΛCDM, but cosmic voids offer a promising avenue (Cautun et al. 2014). Their low densities weaken screening, making them ideal for detecting deviations from GR (Barreira et al. 2015; Baker et al. 2018).

A powerful approach to studying voids is through the weak lensing (WL) phenomenon, which probes the deflection of photons by large-scale structure (LSS). It introduces distinct distortions in background galaxies, quantified by cosmic shear γ and convergence κ (Bartelmann & Schneider 2001; Kilbinger 2015; Ishak 2019; Umetsu 2020) and is highly sensitive to the growth of structure and cosmic expansion (Bartelmann 2010; Troxel & Ishak 2015). The next generation of Stage-IV WL surveys, including Rubin (LSST Dark Energy Science Collaboration 2012), Euclid (Laureijs et al. 2011; Euclid Collaboration: Ajani et al. 2023; Euclid Collaboration: Congedo et al. 2024), and Roman (Spergel et al. 2015), will significantly improve sky coverage and angular resolution, enhancing constraining power. Euclid DR1 (Euclid Collaboration: Mellier et al. 2025), in particular, will increase WL source counts by a factor of 4-5, extending the redshift distribution (Kitching & Heavens 2017), with respect to ground-based surveys like KiDS (Hildebrandt et al. 2017; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023). To fully exploit these data, accurate numerical simulations and high-fidelity modeling techniques are essential.

In voids, the deflection of light occurs outward, opposite to the inward deflection observed in overdense areas. This phenomenon is due to the gravitational effects of spacetime curvature, driven by matter overdensities, which influence the trajectory of light. Additionally, the magnitudes of observed astronomical objects are influenced by the environment and line of sight inhomogeneities, leading to magnification or demagnification effects (see also Clarkson et al. 2012; Bolejko & Ferreira 2012). The unique reversal of gravitational lensing features, such as shear and convergence profiles, is commonly referred to as anti-lensing (Bolejko et al. 2013).

Cautun et al. (2018) compared the predictions of 3D and 2D voids for Euclid and LSST-like lensing surveys, finding that 2D voids are more sensitive to MG theories, showing a deeper lensing signal due to enhanced growth of cosmic structures. Moreover, WL studies applied to voids offer an advantage by reducing the reliance on luminous tracers needed for 3D void identification. In fact, WL is sensitive to the total matter distribution, enabling the detection of voids directly in the matter domain. Two primary methodologies are commonly employed in the literature for analyzing WL around 2D voids:

• WL tunnel voids: This approach measures tangential shear (γ_t) or convergence (κ), which reflect the distortions of background galaxy shapes caused by the projected profiles along the line of sight of extensive underdense structures, typically spanning hundreds of megaparsecs (Higuchi et al. 2013; Gruen et al. 2016; Davies et al. 2021a; Shimasue et al. 2024).

• Void lensing (VL): This method quantifies the excess surface mass density, representing the projection of total matter within a thin lens. This technique is suitable for voids with radii ≤50 h⁻¹ Mpc, where the thin lens approximation holds, and corresponds to studying the WL signal of voids using a tomographic approach (see Boschetti et al. 2024, for further details).

The first method provides less 3D information, but primarily leverages tangential shear for constraining alternative gravity models, as shown by Cautun et al. (2018) and Davies et al. (2019a). In contrast, VL seeks to retrieve dynamic and morphological information about voids, aligning more closely with their standard 3D definitions. However, tunnels exhibit the highest signal-to-noise ratio (S/N) of γ_t in WL maps, producing shear profiles up to ten times stronger (Davies et al. 2018) and offering superior constraints on MG models due to their small size, high abundance, and minimal overlap (Cautun et al. 2018; Davies et al. 2019a, 2021b). For these reasons, in this work we decided to adopt the WL tunnel approach.

This paper is organized as follows. In Sect. 2, we give a theoretical overview of weak gravitational lensing and f(R) gravity models. In Sect. 3, we describe the cosmological simulations employed in this work, followed by the construction of WL light cones and the implementation of ray-tracing techniques to generate convergence maps. In Sect. 4, we present a new code to identify WL tunnel voids, which we applied to the generated convergence maps. Next, in Sect. 5, we introduce and analyze the two primary void statistics explored in this work: the tunnel void size function and the stacked tunnel void tangential shear. In the final part of the paper, Sect. 6, we describe how we adopted two different approaches to model the tangential shear signal. Specifically, we first integrated two established 3D void surface density profiles along the line of sight and then introduced a new parametric model to reproduce the tangential shear signal. Finally, in Sect. 7, we summarize our conclusions.

2. Theory

In this section, we introduce the theoretical background of WL and f(R) MG models to describe the fundamental quantities utilized in this work.

2.1. Weak lensing formalism

The phenomenon of gravitational lensing occurs when light rays from distant sources pass through inhomogeneous matter distributions, experiencing deflections that distort the observed shapes of background sources (Bartelmann & Schneider 2001; Kilbinger 2015). To model this phenomenon, it is assumed that on scales relevant to gravitational light deflection, spacetime is described by a weakly perturbed Friedmann-Lemaître-Robertson-Walker metric.

Assuming the Newtonian gauge and defining Φ and Ψ as the Einstein-frame metric potentials, the perturbed metric is described as:

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=(1+2\mathrm{\Psi })\mathrm{d}{t}^2-a^2(t)(1-2\mathrm{\Phi })\mathrm{d}{x}_i\mathrm{d}{x}_j.\end{array}$$(1)

Here, Ψ represents the gravitational potential, which describes the temporal part of the metric perturbations, while Φ is the gravitational curvature potential, describing the spatial part of the metric perturbations.

Here, we let β represent the true angular position of a source at a comoving line of sight distance χ_s and redshift z_s = z(χ_s), with θ as its observed angular position on the sky. In this framework, the angular positions β and θ are connected via the lens equation, which traces the trajectory of the photon as it is deflected by an angle α through a lens or lens system. It is generalized as:

$$\begin{array}{c}\hfill \mathit{\beta }(\mathit{\theta },z_{\mathrm{s}})=\mathit{\theta }-\mathit{\alpha }=\mathit{\theta }-\frac{2}{c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{s}}}}\frac{f_{\mathrm{ls}}}{f_{\mathrm{l}}{f}_{\mathrm{s}}}{\mathrm{\nabla }}_{\mathit{\beta }}{\mathrm{\Phi }}_{\mathrm{len}}(\mathit{\beta }(\mathit{\theta },{\chi }_{\mathrm{l}}),{\chi }_{\mathrm{l}},z_{\mathrm{l}})\mathrm{d}{\chi }_{\mathrm{l}},\end{array}$$(2)

where c is the speed of light and f_ls = f_K(χ_s − χ_l), f_l = f_K(χ_l), f_s = f_K(χ_s), where we assume the case χ_s > χ_l with “s” and “l” representing the source and the lens, respectively. Here, we assume a generic curvature K, so f_K(χ) represents the comoving angular diameter distance, and f_K(χ)θ corresponds to the transverse comoving vector for a comoving line-of-sight distance χ = χ(z) (Misner et al. 1973; Hogg 1999; Umetsu 2020). Φ_len(β, χ_l, z_l) denotes the 3D lensing gravitational potential generated by a mass distribution ρ(β, χ_l) at redshift z_l (Hilbert et al. 2020).

Generally, Φ_len = (Φ + Ψ)/2, but if we also assume the negligibility of anisotropic stress in the perturbed metric, as is often the case for non-relativistic matter (e.g., cold dark matter), the two metric potentials become equal. Thus, we have:

$$\begin{array}{c}\hfill \mathrm{\Phi }=\mathrm{\Psi }={\mathrm{\Psi }}_{\mathrm{N}},\end{array}$$(3)

where Ψ_N is the Newtonian potential. This equality holds in the standard GR model, where the accelerated expansion of the Universe is driven by a cosmological constant, so we can use Φ_len = Φ under this assumption. Φ is related to the non-relativistic matter density contrast, ${\delta }_m=\delta \rho /\overline{\rho }$, through the Poisson equation:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Phi }(\mathit{\beta },{\chi }_{\mathrm{l}},z_{\mathrm{l}})=4\pi G{a}^2\delta \rho (\mathit{\theta },{\chi }_{\mathrm{l}},z_{\mathrm{l}})=\frac{3{\mathrm{\Omega }}_m{H}_0^2}{2a}{\delta }_m(\mathit{\theta },{\chi }_{\mathrm{l}},z_{\mathrm{l}}),\end{array}$$(4)

where G is the gravitational constant, $\overline{\rho }$ is the mean matter density of the Universe, Ω_m is the total matter density parameter, and H₀ represents the Hubble parameter at the present time.

The image distortion is described by the Jacobian distortion matrix A(θ, χ_s), which represents the relative positions of nearby light rays and is obtained by differentiating the lens equation with respect to θ. Assuming the flat-sky approximation, where β and θ are 2D Cartesian coordinate vectors (Becker 2013), we can express it as

$$\begin{array}{c}\hfill A_{\mathit{ij}}(\mathit{\theta },{\chi }_{\mathrm{s}})\equiv \frac{\partial {\beta }_i(\mathit{\theta },z_{\mathrm{s}})}{\partial {\theta }_j}={\delta }_{\mathit{ij}}-\frac{\partial {\alpha }_i(\mathit{\theta },z_{\mathrm{s}},z_{\mathrm{l}})}{\partial {\theta }_j},\end{array}$$(5)

and the distortion matrix A(θ, χ_s) takes the form:

$$\begin{array}{c}\hfill \mathbf{A}(\mathit{\theta },{\chi }_{\mathrm{s}})\approx \left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2-\omega \\ -{\gamma }_2+\omega & 1-\kappa +{\gamma }_1\\ \end{array}\right),\end{array}$$(6)

where δ_ij is the Kronecker delta. This defines the lensing convergence, κ, the lensing complex shear, γ = γ₁ + iγ₂, and the lensing rotation, ω, which is 𝒪(Φ²) (Petri et al. 2017). In the WL regime, we have |κ|,|γ|,|ω|≪1.

To solve the integral form of Eq. (5), a common approach is to expand it in a series of Φ, retaining only the first-order term, known as the Born approximation (Bartelmann & Schneider 2001). This approximation assumes negligible rotation and disregards lens-lens coupling, allowing the deflection angle α to be expressed as the gradient of a 2D lensing potential, ψ_len. Under this framework, the WL signal arises from the cumulative matter distribution along the line of sight (Davies et al. 2021b). The convergence, κ, follows from the 2D Poisson equation and can be written in terms of the matter density contrast, δ_m, as:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{Born}}(\mathit{\theta },z_{\mathrm{s}})=\frac{3H_0^2{\mathrm{\Omega }}_m}{2c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{s}}}}(1+z_{\mathrm{l}})W({\chi }_{\mathrm{l}},{\chi }_{\mathrm{s}}){\delta }_{\mathrm{m}}(\mathit{\theta },{\chi }_{\mathrm{l}},z_{\mathrm{l}})\mathrm{d}{\chi }_{\mathrm{l}},\end{array}$$(7)

where W(χ_l, χ_s)≡(f_lf_ls)/f_s is defined as the lensing kernel or lensing efficiency factor. This demonstrates that the measured WL convergence can be understood as the integration of the matter density along the unperturbed line of sight, weighted by the lensing kernel.

Moreover, the radial convergence profile of an object κ(r_p) is related to its radial tangential shear profile through the equation (Davies et al. 2018, 2021b):

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{t}}(r_p)=\overline{\kappa }(<r_p)-\kappa (r_p),\end{array}$$(8)

where $\overline{\kappa }(<r_p)$ is the mean enclosed convergence within the radius, r_p, defined as

$$\begin{array}{c}\hfill \overline{\kappa }(<r_p)=\frac{1}{\pi r_p^2}{\displaystyle {\int }_0^{r_p}}2\pi r_p^{\prime }\kappa (r_p^{\prime })d{r}_p^{\prime }.\end{array}$$(9)

Based on the assumption of the flat-sky approximation, here and throughout this work, r_p is used to denote the projected distance from the void center instead of θ.

2.2. f(R) modified gravity models

In this work, we analyzed a class of models called f(R) (Weyl 1918), which is a set of alternative fourth-order scalar-tensor theories of gravity. These theories have been extensively investigated as models deviating from GR (e.g., Utiyama & DeWitt 1962; Starobinsky 1980) due to their potential to explain both the early and late-time acceleration of the Universe (Starobinsky 2007; Hu & Sawicki 2007; Nojiri & Odintsov 2008, 2011). Specifically, they mimic the observed accelerated expansion of the Universe by introducing an additional scalar field that acts as a fifth force, repulsive on large scales. This scalar field replaces the cosmological constant Λ of the standard ΛCDM model and differs from GR in the evolution of density perturbations due to gravitational instability. They also use screening mechanisms to avoid Solar System constraints, aligning with well-tested GR predictions on small scales and in overdense regions.

f(R) theories of gravity extend the Einstein-Hilbert action by incorporating functions of the Ricci scalar, R, which can result in fourth-order field equations unless a constant term is added to the gravitational Lagrangian (Buchdahl 1970). The action for f(R) gravity is given by

$$\begin{array}{c}\hfill S=\frac{M_{\mathrm{P}}^2}{2}{\displaystyle \int {\mathrm{d}}^4}x\sqrt{-g}(R+f(R))+S_m[g_{\mu \nu },\psi ],\end{array}$$(10)

where $M_{\mathrm{P}}^2=1/(8\pi G)$ is the reduced Planck mass and S_m is the action of the matter field, ψ. The function f(R) controls deviations from GR, becoming significant in the low curvature regime (R → 0) and in low-density environments such as voids (Brans & Dicke 1961). GR can be recovered by setting f = −2Λ^GR for a more generalized cosmic acceleration, as given in Carroll et al. (2004).

A prominent f(R) model proposed by Hu & Sawicki (2007) is expressed as

$$\begin{array}{c}\hfill f(R)=-m^2\frac{c_1{\left(\frac{R}{m^2}\right)}^n}{c_2{\left(\frac{R}{m^2}\right)}^n+1},\end{array}$$(11)

where $m^2\equiv H_0^2{\mathrm{\Omega }}_{\mathrm{m}}$ defines the mass scale, while c₁, c₂, and n are non-negative parameters. For ΛCDM consistency, the condition c₁/c₂ = 6Ω_Λ/Ω_m must hold. In the limit c₂(R/m²)ⁿ ≫ 1, the scalar field, f_R ≡ df(R)/dR, can be approximated as

$$\begin{array}{c}\hfill f_R\approx -n\frac{c_1}{c_2^2}{\left(\frac{m^2}{R}\right)}^{n+1}.\end{array}$$(12)

For n = 1, the model simplifies, with f_R0 defined as:

$$\begin{array}{c}\hfill f_{R0}\equiv -\frac{1}{c_2}\frac{6{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}}{\left(\frac{m^2}{R_0}\right)}^2,\end{array}$$(13)

where R₀ is the background value of R at the present time.

By modifying the action in Eq. (10) with respect to g_μν, we obtain the modified Einstein equations for f_R as

$$\begin{array}{c}\hfill f_R{R}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }-{\mathrm{\nabla }}_{\mu }{\mathrm{\nabla }}_{\nu }{f}_R+g_{\mu \nu }\square f_R=8\pi G{T}_{\mu \nu },\end{array}$$(14)

where ∇ is the covariant derivative and □ ≡ g^μν∇_μ∇_ν is the D’Alembert operator. From its trace, we derive the motion equation for the scalar field:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\delta f_R=\frac{a^2}{3}[\delta R\left(f_R\right)-8\pi G\delta \rho ].\end{array}$$(15)

By extracting its time-time component and assuming small perturbations δf_R, δR, and δρ on a homogeneous background and the quasi-static field approximation (slow variation for f_R), we obtain the equivalent of the Poisson equation for scalar metric perturbations, 2ψ = δg₀₀/g₀₀, as

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\psi =\frac{16\pi G}{3}{a}^2\rho -\frac{a^2}{6}\delta R\left(f_R\right).\end{array}$$(16)

Combining Eqs. (15) and (16) allows us to derive exact solutions for extreme cases and analyze the scale and manner in which f(R) gravity diverges from GR.

Within this framework, the behavior of f(R) models varies depending on f_R0 values relative to the effective potential, the Newtonian potential, Ψ_N. If f_R0 ≪ Ψ_N, the behavior of the model closely recovers GR in high-curvature regions due to the screening mechanism. The most notable example of these mechanisms is the so-called Chameleon screening mechanism (Khoury & Weltman 2004a,b). The response of the scalar field to the effective potential, influenced by external matter sources, varies with density, increasing its effective mass in high-density regions and decreasing it in low-density ones (Chameleon field). Conversely, if f_R0 ≫ Ψ_N, gravity would be over-amplified. Therefore, f_R0 should be close in magnitude to Ψ_N, typically between 10⁻⁵ and 10⁻⁶. A value of 10⁻⁴ may still be acceptable when including the contribution of other components like massive neutrinos. In fact, it has been shown that according to the mass associated with these particles, the effect caused by LSS can be nearly counteractive to f(R) models. This happens because, for |f_R0| ≫ |ψ|, the scale size of the scalar field is comparable to the scale of neutrino thermal free-streaming. This creates a degeneracy visible in the matter power spectrum (Saito et al. 2008; Wagner et al. 2012; Cataneo et al. 2015), the halo mass function (Marulli et al. 2011; Villaescusa-Navarro et al. 2013), the clustering properties of CDM halos and redshift-space distortions (Zennaro et al. 2018; García-Farieta et al. 2020). We note that a possible hint towards the disentangling of the combined effects of f(R) models and massive neutrinos was found analyzing the abundance of large voids at high redshifts (Contarini et al. 2021). We investigate the impact of these separate components using WL tunnel voids in Sect. 5.

We go on to focus on the effect of MG on the WL phenomenon. Using the perturbed metric in Eq. (1), WL offers a direct connection to gravity theories through the 3D lensing potential Φ_len described in Sect. 2.1, because it is defined as the sum of the two metric potentials. In contrast to what happens in GR (see Eq. (3)), differentiating these potentials is often crucial for analyzing models, as the additional scalar field may affect photons and matter differently.

In f(R) models, using the general definition of Φ_len in Eq. (4), the two modified Poisson equations for the metric potentials can be derived from the variation of the Einstein equations as

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G{a}^2\delta \rho -\frac{c^2{\mathrm{\nabla }}^2}{2}\delta f_R,{\mathrm{\nabla }}^2\mathrm{\Psi }=4\pi G{a}^2\delta \rho +\frac{c^2{\mathrm{\nabla }}^2}{2}\delta f_R.\end{array}$$(17)

In terms of the Newtonian potential, these can be expressed as:

$$\begin{array}{c}\hfill \mathrm{\Phi }={\mathrm{\Psi }}_N-\frac{c^2}{2}\delta f_R,\mathrm{\Psi }={\mathrm{\Psi }}_N+\frac{c^2}{2}\delta f_R.\end{array}$$(18)

Consequently, the lensing potential remains unaffected as the corrections cancel out, preserving the standard GR form (Hu & Sawicki 2007; Giocoli et al. 2018a). We conclude that in f(R) models, deviations from GR arise not from changes in the lensing potential Φ_len, but only from alterations in the LSS matter distribution. Therefore, we would expect a different distribution of convergence and shear compared to GR predictions. In particular, in regions of low density, the scalar field δf_R contributes to an enhancement of the gravitational interaction, resulting in a more pronounced WL signal around voids than in ΛCDM models.

3. Methods

3.1. DUSTGRAIN-pathfinder simulations

In our analysis, we used data produced in Giocoli et al. (2018a) from the Dark Universe Simulations to Test Gravity In the Presence of Neutrinos (DUSTGRAIN) project, a suite of N-body cosmological simulations.

The DUSTGRAIN-pathfinder runs are cosmological collisionless simulations that track the evolution of 768³ cold dark matter (CDM) particles within a periodic box with a side length of 750 h⁻¹ Mpc. These particles have a mass of $m_{\mathrm{cdm}}^p=8.1\times {10}^{10}{h}^{-1}{M}_{\odot }$ and their gravitational softening has been set to ε_g = 25 h⁻¹ kpc, corresponding approximately to 1/40th of the mean inter-particle separation. All simulations assume a constant Ω_m to ensure matching power spectra on large scales, with and without neutrinos. The reference ΛCDM simulation assumes GR, M_ν = 0 eV. Its cosmological parameters follow Planck2015 constraints (Planck Collaboration XIII 2016), with Ω_m = Ω_cdm + Ω_b + Ω_ν = 0.31345, Ω_Λ = 0.68655, $H_0=67.31\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, n_s = 0.9658, and σ₈ = 0.842 at z = 0.

The alternative cosmological models included in this simulation suite are characterized by deviations from GR in the form of MG of class f(R). The strength of the scalar field is regulated by the parameter fR0, which is set to 10⁻⁴, 10⁻⁵ or 10⁻⁶. These simulations are referred to as fR4, fR5, and fR6, respectively.

Additionally, cosmological scenarios featuring massive neutrinos are also included in the DUSTGRAIN-pathfinder suite, with the goal of exploring the degeneracies between the effects of this hot dark matter component and f(R) gravity. However, in this paper, we want to focus on the isolated signature of neutrinos, leaving the more complex task of disentangling competing effects for future work. For this purpose, we consider a standard ΛCDM simulation with neutrinos of a mass M_ν = 0.15 eV, which we refer to as ΛCDM_0.15eV. In this case, we tracked the evolution of 2 × 768³ particles considering m_cdm^p = 8.01 × 10¹⁰ h⁻¹ M_⊙ and m_ν^p = 9.25 × 10⁸ h⁻¹ M_⊙ as dark matter and massive neutrino particle masses, respectively.

The initial conditions of the DUSTGRAIN-pathfinder simulations were generated using transfer functions calculated by the CAMB Boltzmann solver (Lewis et al. 2000) at the initial redshift z_i = 99. We clarify that they matched across all cosmological models, namely, they share identical random phases in their initial density fields. In the case of simulations including massive neutrinos, two correlated Gaussian random fields were used to initialize CDM and neutrino components from the same initial phase realization, following the methods described in Zennaro et al. (2018), Villaescusa-Navarro et al. (2018). This setup was deliberately chosen to ensure consistency in the large-scale modes and to allow a direct and faithful comparison across cosmologies, under controlled variance conditions.

To summarize, this work employed five simulation sets: a ΛCDM reference model, three MG scenarios with different fR0 values, and a standard cosmology featuring massive neutrinos. The main characteristics of these simulations are presented in Table 1. During each simulation, a sequence of 34 full comoving snapshots was stored, each representing the specified comoving volume at a particular cosmological epoch.

3.2. Weak lensing light cones

The construction of past light cones is crucial for accurately simulating WL effects produced by the total projected matter density distribution. In this work, we used a post-processing reconstruction method, assuming the flat-sky approximation to slice particle snapshots by their comoving distances from the observer (see Shirasaki et al. 2017). Specifically, we processed each DUSTGRAIN-pathfinder simulation set using the MAPSIM routine (Giocoli et al. 2014, 2017), which operates in two main steps, known as i-MAPSIM and ray-MAPSIM.

In the initial phase, i-MAPSIM, particle positions from simulation snapshots within the specified field of view (FOV) are projected onto various lens planes along the line of sight. We stacked 21 of the 34 available snapshots from z_min = 0 to z_max = 4. The light cones have a square base with sides of 5 deg, corresponding to the angular size of the simulation box at z_max, and cover a total sky area of 25 deg². Each light cone is shaped as a pyramid, with the observer at the vertex at z = z_min = 0 and the square base at the comoving distance χ(z_max), as shown in Fig. 1. To achieve high-resolution maps of the projected matter density, the stacked snapshots were divided into 27 lens planes. Since gravitational lensing depends on the projected matter density along the line of sight, each particle was assigned to the closest lens plane based on its moving distances, preserving its angular position within the specified FOV aperture. Simulations with massive neutrinos incorporate this extra component consistently.

Fig. 1. Example of light-cone construction with MAPSIM from z = 0 to z = 4 for a ΛCDM simulation with a box size of 750 h−1 Mpc, tiled to encompass the 5 × 5 deg2 light cone. The light cone geometry is depicted by the red solid line, expanding until its transverse comoving size matches that of the simulation box. Different colors represent different boxes, while the total of 21 different snapshots are defined by the dashed lines. The source plane (blue solid line) for ray tracing is placed at zs = 1. The snapshots within zobs < z < zs are the 13 utilized ones, while the gray region indicates the snapshots with z > zs that are not considered in this study. The yellow solid lines represent the 27 lens planes used, onto which the particles are projected. Among these, only those within the considered snapshots are utilized for constructing the convergence maps used in this work, positioned from z = 0.05 to 0.98, respectively.

The mass density was then interpolated from these projected particle positions to a 2D grid using a triangular-shaped cloud scheme. The grid pixels are designed to have the same angular size across all lens planes. The angular surface mass density Σ(θ) on the l-th lens plane perpendicular to the (0, 0, 1) direction, at a pixel with coordinate indices (i, j) is computed as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{(i,j)}^l(\mathit{\theta })=\frac{1}{A_{\mathrm{pix}}^l}{\displaystyle \sum _{k=1}^n}{m}_k,\end{array}$$(19)

where A_pix^l = 4π/N_pix is the comoving pixel area in steradians for the l-th lens plane, and m_k represents the masses of particles within the pixel, satisfying f_K(χ_l) < f_K(χ_l′), ∀l′≠l and χ_k < χ_s. In this way, a discretization of the 3D density distribution in several mass maps can be obtained up to the selected source plane z_s. The construction of light cones from simulations and the projection of the mass distribution onto individual planes was done at the same time by the MAPSIM tool in this first step.

To enhance the WL statistics, we generated 256 semi-independent light-cone realizations for each cosmological model. This was achieved by randomizing the comoving simulation boxes across the redshift range z_obs < z < z_s by applying the following transformations:

• Inverting the sign of the cartesian coordinates;

• Translating the position of the observer within the box;

• Permuting the coordinate axes.

This procedure preserves the clustering properties of the particle distribution in a snapshot, ensuring the statistical validity of the light cones (Roncarelli et al. 2007). MAPSIM algorithm enables the storage of halo and sub-halo catalogs associated with each randomization (Castro et al. 2018). The constructed light cones utilize approximately one-third of the simulation box volume, repurposing unused regions at lower redshifts and maximizing the utility of the simulation data (Jain et al. 2000). As shown in Euclid Collaboration: Ajani et al. (2023), 85 out of the 256 light cones can be considered fully independent, each sharing less than 1% of DM halos across different lines of sight at z_s = 1 (Sect. 3.3).

Reusing the same simulation box for multiple light cones may slightly underrepresent rare events in the convergence tails, especially in high- or low-κ observables, such as peak counts or voids (Garrison et al. 2019; Rasera et al. 2022; Chen & Yu 2024). Although randomization reduces sample variance, residual correlations persist but are largely mitigated by galaxy shape noise (Sect. 3.4), which suppresses systematic variance differences (Liu et al. 2015).

3.3. Convergence maps

To accurately simulate the weak gravitational lensing signal, it is essential to trace the paths of light rays across the multiple lens planes constructed from our light cones. The mass maps produced in the first step i-MAPSIM were used as input to perform multiplane ray-tracing calculations through ray-MAPSIM, as done in several studies such as Petri et al. (2016, 2017), Giocoli et al. (2017, 2018a), Hilbert et al. (2020). In this subsequent phase, the routine constructed the lensing convergence map using the Born approximation by summing the surface mass density from each lens plane along the line of sight, weighted according to the lensing kernels W(χ_l, χ_s). More details are given in Sect. 2.1.

In the Born approximation, deflections are integrated along a straight-line path rather than by iteratively displacing rays at each lens plane. This method keeps the rays naturally aligned with the simulation grid, avoiding the need to interpolate the projected matter density at the exact ray positions (Hilbert et al. 2020). Unlike the exact multiplane ray-tracing method, the Born approximation avoids solving the Poisson equation (Eq. 4), making it computationally more efficient (Giocoli et al. 2016). Studies such as Giocoli et al. (2017) and Castro et al. (2018) confirm that this method provides an accurate estimation of the convergence power spectrum and the probability distribution function (PDF) down to sub-arcminute angular scales. Moreover, Schäfer et al. (2012) showed through a perturbative expansion that the Born approximation remains highly reliable for WL analyses even at very small scales (l ≥ 10⁴) (Giocoli et al. 2018a). Although post-Born corrections improve multiplane ray-tracing accuracy in such cases, their impact on void statistics is negligible (Ferlito et al. 2024).

In this framework, to compute the convergence map for each lens plane, we used the discretized form of Eq. (7)¹. Introducing the critical surface density in comoving units for a lens plane at redshift z_l and sources at z_s > z_l, defined as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit},\mathrm{l},\mathrm{s}}^c(z_l,z_s)\equiv \frac{c^2}{4\pi G}\frac{f_{\mathrm{s}}a(z_l)}{f_{\mathrm{l}}{f}_{\mathrm{ls}}}=\frac{c^2}{4\pi G}\frac{a(z_l)}{W({\chi }_s,{\chi }_l)}.\end{array}$$(20)

we omputed the convergence on the l-th lens plane as²

$$\begin{array}{c}\hfill {\kappa }^l(\mathit{\theta },{\chi }_{\mathrm{s}})=\frac{{\mathrm{\Sigma }}^l(\mathit{\theta })}{{\mathrm{\Sigma }}_{\mathrm{crit},l,s}^c}.\end{array}$$(21)

At the end, the total convergence map is

$$\begin{array}{c}\hfill \kappa (\mathit{\theta },{\chi }_{\mathrm{s}})={\displaystyle \sum _l}{\kappa }^l(\mathit{\theta },{\chi }_{\mathrm{s}})={\displaystyle \sum _l}\frac{{\mathrm{\Sigma }}^l(\mathit{\theta })}{{\mathrm{\Sigma }}_{\mathrm{crit},l,s}^c}.\end{array}$$(22)

This ray-tracing configuration precisely simulates the gravitational lensing effects caused by underdense regions in the Universe, facilitating an in-depth examination of lensing profiles from cosmic voids in the context of MG models. As shown in Fig. 1, in this work, the convergence maps were constructed from the corresponding light cones, emulating observational conditions by placing the background sources plane at z_s = 1. This setup provides a continuous matter field over the specified FOV, segmented into 2048 × 2048 pixels, yielding a pixel resolution of approximately 9 arcsec. The origin of each map coordinate system (RA, Dec.) is positioned at the center such that each side spans from −150 to +150 arcmin. The redshift of the background source has been chosen because it is estimated to be the peak of the source probability distribution function, n(z_s), for the Euclid wide photometric survey (Euclid Collaboration: Scaramella et al. 2022; Euclid Collaboration: Ajani et al. 2023; Euclid Collaboration: Mellier et al. 2025).

3.4. Galaxy shape noise

Galaxy shape noise (GSN) represents a fundamental challenge in WL studies. It arises from the intrinsic random distribution in ellipticities and orientations of background galaxies, which heavily dominate the observed correlations in galaxy shapes caused by gravitational lensing (Kilbinger 2015). GSN introduces a stochastic noise component that causes significant variability in the reconstruction of convergence maps, which must be addressed to accurately identify under- (valleys) and over-densities (peaks) regions in WL fields (Lin & Kilbinger 2015). To simulate observational conditions, GSN is added to the simulated WL convergence maps as a Gaussian random field, specifically by superimposing random values drawn from a Gaussian distribution for each pixel (van Waerbeke 2000).

On each pixel of the map κ(θ), we added GSN, n(θ), modeled as a Gaussian random field with a top-hat filter with a size that corresponds to the pixel area $A_{\mathrm{pix}}={\theta }_{\mathrm{pix}}^2$ in arcmin⁻². According to Lin & Kilbinger (2015) its variance is given by

$$\begin{array}{c}\hfill {\sigma }_{\text{pix}}^2=\frac{{\sigma }_{ϵ}^2}{2}\frac{1}{n_{\mathrm{gal}}{A}_{\mathrm{pix}}},\end{array}$$(23)

where ${\sigma }_{ϵ}^2=\langle {ϵ}_1^2\rangle +$$\langle {ϵ}_2^2\rangle$ is the variance of the intrinsic ellipticity distribution of the source galaxies and n_gal is the source galaxy number density. We assumed a Euclid-like setup (Euclid Collaboration: Scaramella et al. 2022; Euclid Collaboration: Ajani et al. 2023; Euclid Collaboration: Mellier et al. 2025), so ${\sigma }_{ϵ}^2=0.3$ and $n_{\mathrm{gal}}=30{\mathrm{arcmin}}^{-2}$ at z_s = 1. In this way, we obtained a noised convergence map κ_n(θ) = κ(θ)+n(θ).

The inclusion of GSN significantly contaminates WL maps with noise (Lin & Kilbinger 2015). This affects key statistics such as void abundance and shear profiles (Davies et al. 2019b, 2021b). The GSN, however, can be reduced by applying a Gaussian smoothing filter characterized by a smoothing length θ_G (Kilbinger 2015). In this work, we used this approach to minimize GSN contamination in tunnel void statistics, ensuring consistency between the maps, with and without GSN.

Choosing an optimal value for θ_G depends on the specific analysis, but it usually ranges between 1 and 5 arcmin. Previous studies have shown that cosmological constraints derived from WL peaks can be improved by employing multiple smoothing scales (Liu et al. 2015). However, it has also been demonstrated that a single intermediate smoothing scale can provide a balanced approach for specific WL peak analyses while remaining robust against GSN contamination (Davies et al. 2019b). Selecting a small θ_G retains information from smaller scales within the WL map but leaves a significant level of GSN contamination. In contrast, a large θ_G suppresses noise more effectively but also erases small-scale features. Thus, selecting an optimal value requires balancing noise reduction with the preservation of meaningful WL signals.

In this work, we initially tested three smoothing scales: θ_G = 1, 2.5, and 5 arcmin. To evaluate the robustness of void detection under different noise conditions, we considered each case with and without GSN. The results align with Davies et al. (2019a) and Davies et al. (2021b), showing that small smoothing scales (θ_G = 1 arcmin) preserve fine structures but amplify noise, increasing variance in WL void statistics. Large scales (θ_G = 5 arcmin) reduce noise but erase meaningful features, lowering void detection sensitivity. Intermediate smoothing scales provide the best balance, minimizing noise while preserving statistical significance. Moreover, as demonstrated by Weiss et al. (2019), such intermediate scales help to reduce the discrepancies between the employed simulations (which, in our case, do not include baryonic physics) and fully hydrodynamic simulations. For these reasons, we finally decided to carry on the analysis using θ_G = 2.5 arcmin. Thus, we obtained the noised and smoothed convergence map as κ_2.5(θ)≡(κ_n*U)(θ) using

$$\begin{array}{c}\hfill U(\mathit{\theta })=\frac{1}{\pi {\theta }_{\mathrm{G}}}exp(-\frac{{\theta }^2}{{\theta }_{\mathrm{G}}^2}),\end{array}$$(24)

with U(θ) being the Gaussian window function with θ_G = 2.5 arcmin. This choice is particularly suitable for analyzing void shear profiles while retaining sufficient small-scale information within WL maps.

We assumed that intrinsic ellipticities are uncorrelated between source galaxies. In this way, the noise after the smoothing can be also described as a Gaussian random field (Bond & Efstathiou 1987). According to van Waerbeke (2000), its variance is related to the number of galaxies contained in the filter as

$$\begin{array}{c}\hfill {\sigma }_{\text{noise}}^2=\frac{{\sigma }_{ϵ}^2}{2}\frac{1}{2\pi n_{\mathrm{gal}}{\theta }_{\mathrm{G}}^2}.\end{array}$$(25)

This relation allows the derivation of the lensing S/N(θ)³ map of the noised and smoothed convergence map κ_2.5(θ) as

$$\begin{array}{c}\hfill \text{S/N}(\mathit{\theta })=\frac{{\kappa }_{2.5}(\mathit{\theta })}{{\sigma }_{\text{noise}}({\theta }_{\mathrm{G}})},\end{array}$$(26)

where σ_noise(θ_G) is the standard deviation of the smoothed GSN map (without contributions from the WL convergence map, i.e., noise only). It varies depending on the smoothing scale used to identify the WL peaks or minima, whereas it is constant only if we assume that galaxies are uniformly distributed. In our framework, with θ_G = 2.5 and z_s = 1, σ_noise ≈ 6.18 ⋅ 10⁻³.

Figure 2 illustrates the impact of applying the Euclid-like GSN and a Gaussian smoothing scale of θ_G = 2.5 arcmin to WL convergence maps. The top panel displays the noiseless convergence field for one of our 256 light cones, while the bottom panel presents the corresponding S/N map. The smoothing enhances the detectability of large-scale structures by reducing small-scale noise while preserving the primary features of the field. As observed, the convergence field exhibits a complex distribution of over- and underdensities, while the S/N map highlights significant structures by suppressing uncorrelated noise.

Fig. 2. Example of a realization of the ΛCDM light cone, having a FOV of 5 × 5 deg2 aperture and zs = 1. In the upper panel, we show the WL convergence map, and in the lower panel, the S/N of the same map with smoothing scale at θG = 2.5 arcmin. The color scale in both panels represents the convergence value for each pixel in the top panel and the S/N value in the bottom panel, as indicated by the corresponding color bars.

3.5. Probability density function and power spectrum

The statistical properties of the convergence field provide crucial insights into the underlying cosmological model. Therefore, we went on to analyze the one-point PDF and the angular power spectrum of the constructed maps.

The PDF of the convergence field quantifies the distribution of κ values across the sky and is particularly sensitive to non-Gaussian features induced by structure formation. Figure 3 shows the mean PDF computed from 256 light-cone realizations for all the analyzed cosmological scenarios, together with the mean PDF noiseless for the ΛCDM model. The noised and smoothed PDFs exhibit a subtle dependence on the cosmological model, while the difference is more pronounced when comparing the noised and smoothed ΛCDM PDF with its noiseless counterpart.

Fig. 3. Top panel: Average PDFs of the noised and smoothed convergence field from 256 light-cone random realizations with redshift zs = 1, for all the analyzed models: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red) and ΛCDM0.15 eV (purple). The black dashed line represents the average PDF without the addiction of GSN for the $\mathrm{\Lambda }\mathrm{CDM}$ model. Bottom panel: PDF noised and smoothed corresponding residuals computed with respect to the ΛCDM model. For the purposes of visualization, only the uncertainties of the noiseless and the reference ΛCDM measurements are reported in this plot as hatched gray and shaded gray areas, respectively.

The bottom panel of Fig. 3 illustrates the residuals of the noisy and smoothed PDFs with respect to the standard ΛCDM model. In particular, to calculate the residuals of any quantity X analyzed in this work, we adopted the following general formula

$$\begin{array}{c}\hfill \mathrm{\Delta }X=\frac{X^{\mathrm{model}}(r)-X^{\mathrm{ref}}(r)}{|X^{\mathrm{ref}}(r)|},\end{array}$$(27)

where X^model(r) represents the profile or observable under investigation, while X^ref(r) corresponds to a fixed reference model for comparison. Here, the errors are computed as the standard deviation of the measurements across different realizations of the maps. Throughout the paper, the errors are properly propagated during the calculation of the residuals.

Looking at the residuals, we can appreciate how the deviations from the ΛCDM model become more relevant for κ < 0. Effectively, the largest deviations occur in the negative κ range, where f(R) models display an excess probability compared to ΛCDM, while the massive neutrino model exhibits a suppression. The lower PDF values in underdense regions strengthen fractional deviations. MG models amplify structure formation, deepening underdensities and enhancing non-Gaussian features. Massive neutrinos, conversely, suppress structure growth, shifting the PDF toward lower contrast. Since these effects are more pronounced in the negative κ tail, we expect WL from tunnel voids to represent a key tool for testing alternative cosmologies as they can capture the differential impact of MG and neutrinos.

The angular power spectrum of the convergence field, P_κ(l), provides a complementary statistical measure by characterizing the variance of the κ maps fluctuations, describing how power is distributed across different angular scales. In Fig. 4, it is shown as a function of the multipole moment l, which corresponds to the inverse angular scale on the sky normalized for the field size (2π/θ_FOV). It is computed as the Fourier transform of the input map squared and binned in multipole space.

Fig. 4. Average angular power spectra of the noised and smoothed convergence field from 256 light-cone random realizations with redshift zs = 1, for all the analyzed models (top panel) and the correspondent residuals with respect to the reference ΛCDM model (bottom panel). Line styles and color schemes are the same as used in Fig. 3. The dotted vertical line marks the angular mode of the characteristic smoothing scale in the Fourier transform of the κ field.

In particular, Fig. 4 shows the mean power spectrum computed from 256 convergence maps for all the analyzed cosmological scenarios alongside the noiseless mean power spectrum for the ΛCDM model. The difference between the models can be better appreciated in the bottom panel, where the relative differences of the power spectra with respect to ΛCDM are presented. The largest deviations from the ΛCDM model occur at the intermediate scales, in the multipole range [300 − 1000]. This range corresponds to angular scales between [11.46, 3.44] arcmin. The impact of MG and massive neutrinos is maximized in this regime. f(R) models enhance structure formation, leading to an increase in power on these scales. Massive neutrinos suppress clustering, causing a damping of power. We can notice how these competing effects lead to significant deviations in the convergence power spectrum.

At larger scales (l ≲ 300), the Universe remains in the linear regime, where deviations are smaller due to the weaker impact of nonlinear modifications. At smaller scales (l ≳ 1000), nonlinear effects dominate, and Gaussian smoothing starts to suppress the structure-dependent variance.

In Fig. 4, the vertical line marks the characteristic smoothing scale in Fourier space (l = 1375.10). Beyond this limit (θ_G = 2.5 arcmin), the signal is attenuated and significantly deviates from the noiseless power spectrum, reducing sensitivity to cosmological model differences. This behavior motivates us to investigate in depth the effect of WL on larger scales and, given the findings from the analysis of the convergence PDF and power spectrum, to focus especially on voids.

4. 2D tunnel void finder

One of the main challenges in the study of cosmic voids lies in their definition and, consequently, in their identification. Cosmic voids are not luminous structures; instead, they can be defined only by the distribution of matter tracers, such as galaxies, which predominantly reside along their boundaries. This necessitates the use of ad-hoc techniques to reconstruct void shapes and centers from the tracers, called void finder algorithms. The absence of a universally accepted, first-principles definition of voids makes their identification via void finders highly dependent on the specific goals and type of analysis being performed (Colberg et al. 2008; Cautun et al. 2018; Davies et al. 2019a, 2021b).

In this work, we developed a novel 2D void finder algorithm for WL analyses, guided by these insights, incorporating features from existing methods in the literature while addressing their limitations to enhance both stability and effectiveness. Specifically, we implement a new tunnel finder based on WL absolute minima, adding GSN to our maps. We have made it publicly available (see Sect. ‘Data availability’). The structure and the main steps of the code are schematically represented in Fig. 5 and analyzed in detail in the following.

Fig. 5. Diagram showing the four main steps of the 2D void finder algorithm that allows the identification of voids from a WL convergence map.

4.1. PyTwinPeaks: Reconstruction of connected regions

The processed convergence maps serve as the input for our 2D void finder algorithm, which detects underdense regions based on their S/N. The algorithm called PyTwinPeaks is a Python version of TwinPeaks c++-based code (Giocoli et al. 2018b), a peaks/valleys finder that analyzes overdense and underdense S/N regions. In this work, we used PyTwinPeaks to reconstruct underdense lines-of-sight pixel regions at different convergence S/N thresholds, improving its efficiency. This algorithm was later expanded and transformed into a dedicated WL tunnel void finder.

The first phase of the algorithm consists of a reconstruction of the valleys: it begins with the creation of a mask that selects pixels in the S/N map below the chosen threshold. Next, the reconstruction is done via an affiliated component analysis, identifying pixels belonging to the same connected region. The algorithm builds up a total map as the sum of layers, where each layer represents a connected region with a S/N value that is below the selected threshold.

From each connected region, the code extracts 16 topological features, including perimeter, area, and eccentricity. It also records the (x, y) coordinates of the weighted centroid, computed using the S/N values of individual pixels as weights, and the coordinates of the local minimum, defined as the pixel with the lowest S/N value. The coordinates, originally in pixels, are converted to arcminutes using the relation

$$\begin{array}{c}\hfill x_{\mathrm{arcmin}}=(x_{\mathrm{pix}}-\frac{n_{\mathrm{pix}}}{2})·l_{\mathrm{arcmin}},\end{array}$$(28)

where x_pix represents the coordinate in the S/N map in pixel units, n_pix is the number of pixels per side of the map, and l_arcmin = FOV_arcmin/n_pix denotes the pixel size in arcminutes.

Once the underdense regions are reconstructed based on a specified threshold, overlapping cases between nearby regions are handled. Each connected region is initially approximated as a circle with a preliminary radius R_pr, which corresponds to the equivalent radius of the circle with the same area. Although the regions are disjoint by construction, their circular representation may lead to overlaps. Two regions merge into a larger structure if the distance between them satisfies the condition

$$\begin{array}{c}\hfill D_{\mathrm{pr}}<\frac{R_{\mathrm{pr},1}+R_{\mathrm{pr},2}}{2}.\end{array}$$(29)

The algorithm repeats this procedure over a given threshold range, generating a catalog of underdense regions. The number of thresholds tested is adjustable and affects both the accuracy of the search and the computation time. To identify the connected pixel regions resulting from the projection of underdense lines-of-sight, we analyzed 11 threshold levels from S/N = −5 to S/N = 0 with a step of 0.5. The resulting catalog contains the properties of each connected region for each threshold.

4.2. Void center assignment

This step is crucial, as it determines the final catalog of void centers in the projected S/N field and significantly affects the resulting statistics. As mentioned earlier, at the end of the first phase, the catalog includes the positions of the local minima, the weighted centroids, and the areas of each connected region across different thresholds. Each threshold produces a different selection of voids, and a center can be assigned to each connected region identified at a given threshold.

Among the many void finders in the literature, the two primary approaches for assigning the void center (including the ones explored in this work) are through the centroid (see Padilla et al. 2005; Platen et al. 2007; Neyrinck 2008; Sutter et al. 2015) and the local minimum (see Sánchez et al. 2017; Nadathur et al. 2019; Davies et al. 2021b; Hang et al. 2021; Kovács et al. 2022) of the connected region. In Fig. 6, we illustrate the hierarchical structure of voids identified across all thresholds, using underdense regions centered on weighted centroids (left panel) and local minima (right panel), colored accordingly for different S/N values. The void sizes shown here are not the final results of our finder; their dimensions are determined by the area of the region at the highest threshold, which does not fully capture the complex morphology of a tunnel void.

Fig. 6. Left panel: Graphical representation of all 2D voids with their preliminary radii of each threshold produced by the same map shown in Fig. 2, centered at the S/N-weighted centroids of the connected regions. Right panel: Same as the left panel, but with voids centered at the local minima of the connected regions. In both, each void of a given threshold is represented with the colored area of the same color, this color varies depending on the S/N threshold as shown in the colorbar.

Each approach has its strengths and limitations, and the choice between them depends on the specific use case and scientific objective. For example, although the centroid method weights the S/N values, the assigned center may occasionally fall outside the connected region. At higher (still negative) S/N values, connected regions can have elongated or irregular shapes, and in extreme cases (e.g., a crescent-shaped underdensity), the centroid may fall within or too close to a positive S/N (overdense) region. This misplacement contaminates the tangential shear profile and distorts the extracted signal. Moreover, as shown in Fig. 6, the centroid method introduces a continuous shift of void centers depending on the selected S/N threshold, which affects the total number of 2D voids. In fact, the positions of void centers directly influence how underdense regions merge. On the other hand, centroids help reduce the impact of GSN at lower S/N values. Using instead local minima as void centers maximizes the tangential shear signal and resolves the topological issues associated with centroids. However, voids centered on local minima are more susceptible to noise fluctuations (Davies et al. 2021b).

In light of these considerations, in the new version of PyTwinPeaks, we developed a method for void center assignment based on what we call absolute minima. This strategy is inspired by a common technique in 3D void finders known as the watershed algorithm (see Platen et al. 2007; Nadathur & Hotchkiss 2015), here adapted for 2D voids. The algorithm tracks local minima across multiple thresholds, starting from the lowest and checking whether each minimum persists (i.e., whether it falls within the area outlined by the next threshold) as the threshold increases up to a user-defined limit. A minimum that persists across all levels is classified as an absolute minimum.

The initial threshold sets the depth at which void centers are identified and directly impacts the depth of the final voids. By construction, this depth affects the amplitude of the tangential shear signal extracted from those regions. A balanced choice ensures a statistically significant number of 2D voids while preserving the strength of the WL signal. The final threshold is equally important. A low value causes absolute minima to coincide with local minima, creating small and fragmented voids. Conversely, a high value enlarges underdense regions along the line of sight, leading to the excessive merging of multiple structures. The final threshold must be carefully chosen to ensure reliable void statistics while avoiding overmerging. Contrary to the example shown in Fig. 6, after extensive testing, we found that the range S/N = [ − 4, −3] is optimal for setting the minimum and maximum selection thresholds: this technique produces tangential shear profiles with the deepest signals and the smallest uncertainties. The latter is computed from the scatter among different void profiles in the catalog, as we discuss in Sect. 5.

For the analysis presented in this work, the implementation of this method proved effective in resolving the limitations of the two previously mentioned approaches. Absolute minima in the S/N field define the void centers, ensuring an optimal shear signal. GSN contamination is simultaneously mitigated, together with map smoothing, by analyzing underdensities across multiple thresholds and checking for center persistence in higher S/N regions. The final catalog generated by this phase consists of underdense regions centered on a persistent minimum identified at the lowest threshold, with sizes determined by the area of the region at the highest threshold.

4.3. Radius assignment and cleaning

The next phase of the developed void finder algorithm assigns the final radius to the void centers identified in previous steps. To adapt this radius to the statistical signal of underdense lines-of-sight, such as tunnel voids, we impose each void to correspond to a circle enclosing the exact S/N final value. At the same time, the radius must be adjusted to account for the spatial resolution and boundary constraints of the convergence map.

The algorithm initializes a grid of 3 × 3 pixels around each void center, whose coordinates are converted from arcminutes to pixels using the inverse formula of Eq. (28). A circular mask⁴ is applied to this grid, selecting all pixels enclosed by or intersected by the mask. The mask radius is set according to the relationship

$$\begin{array}{c}\hfill r_{\mathrm{mask}}=\frac{l_{\mathrm{grid}}}{2}-\frac{l_{\mathrm{pix}}}{2},\end{array}$$(30)

where l_grid represents the grid size and l_pix denotes the pixel size, both measured in pixel units. The mask ensures high precision and computational efficiency by restricting the selection of relevant pixels.

The algorithm expands the grid by adding a row of pixels on each side until the average S/N within the circular mask exceeds the selected final threshold. The final radius is determined by interpolating between the last two values, refining the estimate of the radius that precisely corresponds to the considered threshold. As the final threshold for the expansion of the pixel grid in the ray assignment, we used the value ⟨S/N⟩= − 2.5 (see Sect. 5.2).

For voids located near the boundaries of the S/N map, removing them from the catalog would reduce statistical significance. To address this, the algorithm masks grid cells that extend beyond the boundaries of the map. Pixels outside the map are assigned a “not a number” (NaN) value, ensuring they do not influence the averaging process. This approach effectively extends the map while preserving the integrity of the void identification procedure.

We were then able to refine the void catalog by removing overlapping voids. With the final radius assigned to each underdensity, the last requirement is to eliminate voids that are not independent, as WL signals extracted from overlapping voids cannot be treated separately. Since an initial overlap check was performed earlier, this refinement is expected to have a limited impact on the final catalog.

The identified voids are sorted in descending order according to their size. The algorithm checks whether the distance between the centers is less than or equal to 75% of the sum of the two radii: D_centers ≤ 75%⋅(r_p, 1 + r_p, 2), where r_p, 1 and r_p, 2 represent the projected radius of the two voids. In this way, we can also correct for the void-in-void scenario, which happens when a void is completely contained in a larger one. The updated void list, ordered by size, prevents unnecessary removal of underdensities. In this way, when multiple voids overlap in sequence, only the central one is removed, preserving the largest and smallest, which remain non-overlapping.

The final output of the algorithm is an ASCII file containing three columns. Each row stores the void center positions in the X and Y directions and the void radius, all measured in arcminutes.

Figure 7 shows, as an example, a S/N map from one of the 256 convergence maps produced for the ΛCDM scenario, where circles mark the positions of the 2D voids. We can note that these voids correspond very well with the negative S/N regions. All negative regions on the map that are not associated with 2D voids are either considered part of a larger, spatially nearby void or are considered to be caused by S/N fluctuations. We underline how the chosen criterion for assigning the void radius is particularly effective in identifying voids as entirely negative regions surrounded by positive areas; namely, the void sizes result correctly expanded until they reach a WL positive signal. Additionally, this method optimizes the total computational time required. The void finder takes about 1 minute to produce the final void catalog from a single map and around 4 hours to analyze all 256 maps, running on a local machine with 12 cores, a maximum frequency of 4.7 GHz, and a CPU scaling⁵ at 13%.

Fig. 7. Example of the final tunnels catalog associated with the same S/N map of Fig. 2 and Fig. 6. Colors from blue to red represent areas with increasing S/N, as indicated by the colorbar on the right. The 2D voids identified by our algorithm are shown as areas delimited by black circles. The corresponding centers (absolute minima) are marked with yellow crosses.

We also emphasize that the tunnel voids identified by this new 2D void finder result from a specific combination of choices. These include the selection of the GSN, the smoothing scale, the threshold range for the watershed method to distinguish absolute minima, and the final threshold for the radius assignment. Despite the presence of free parameters in the algorithms, which makes the final catalog susceptible to user choices, this flexibility makes the algorithm particularly suitable for application to real data from ongoing and future surveys.

5. Tunnels void statistics

5.1. 2D void size function

Compared to 3D voids identified in galaxy surveys, generally covering scales of tens of megaparsecs (see, e.g., Contarini et al. 2023), WL tunnel voids extracted in this work from S/N maps are characterized by relatively small sizes. Their angular extension in the sky-plane spans from 2 to 15 arcmin.

The reduction in void size is expected due to the nature of WL tunnel voids. These voids form when photons from distant sources primarily pass through underdense regions, particularly at redshifts near the peak of the lensing kernel function. The observed underdensities result from multiple 3D voids that partially align with the observer’s line of sight. The projection onto the sky plane does not preserve the original size of 3D voids. Instead, it highlights only the areas where the signal is dominated by negative convergence.

In Fig. 8 we show the tunnel void size function (VSF), that is, the total number counts of 2D voids as a function of their projected radius R_v in arcminutes, averaged from the 256 S/N map realizations, for each cosmological model. We restrict our analysis to the range [2, 15] arcmin to avoid the statistically rarest objects, which are characterized by large uncertainties. The error associated with void number counts is estimated as the scatter across independent realizations, normalized by the square root of the number of maps. On the top panel, we present the five VSFs, and, in the bottom panel, the corresponding residuals with respect to the ΛCDM model, computed as in Eq. (27).

Fig. 8. Top panel: Average void size functions, per degree squared, of 2D WL voids from 256 light-cone random realizations with redshift zs = 1, for the different analyzed models: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), and ΛCDM0.15 eV (purple). Bottom panel: Corresponding residuals computed with respect to the ΛCDM model. In both panels, the shaded regions of the same colors around the curves indicate errors, computed as the standard deviation of counts in each bin divided by the number of maps.

The shape of all the VSFs is consistent with the one measured for 3D voids (Hamaus et al. 2016; Ronconi & Marulli 2017; Contarini et al. 2022, 2023), namely, it shows how small voids are more abundant with respect to the larger ones. The void radius, defined on the basis of the S/N in the 2D convergence map, reflects changes in the void interiors. Specifically, MG leads to larger, deeper voids, while massive neutrinos produce smaller, shallower ones (Contarini et al. 2021). This shifts the VSF accordingly, especially at large radii. At the small end, however, the peak remains fixed due to resolution limits: pixel size and smoothing set a lower bound on the detectable void size. A minimum radius cut and increased void mergers at large radii further impact the total tunnel void counts, leading to significant cosmology-dependent variations in VSF amplitude.

By analyzing different cosmological models, we note that all the size distributions are characterized by a similar shape and the main difference is in their amplitude. The observed trend for different cosmological scenarios aligns with expectations: in MG models, higher values of the parameter |f_R0| lead to a greater increase in void counts. In the presence of massive neutrinos, the abundance of voids is suppressed. This happens because a stronger scalar field action, associated with higher |f_R0| values, enhances the evolution of LSS, including cosmic voids. Conversely, larger neutrino masses suppress the growth of cosmic structures, reducing void formation. Consequently, the models fR6 and ΛCDM_0.15 eV are expected to be the most similar to the ΛCDM case, although they exhibit opposite trends. Specifically, up to R_v ∼ 10 arcmin, we observe an amplitude difference of ≤10% for the fR6 model, while for ΛCDM_0.15 eV, the difference is approximately 25%.

5.2. Stacked tangential shear profile

The tangential shear profile of each WL tunnel void is extracted by expanding 64 concentric circular shells around each void center in the noised and smoothed convergence map κ_2.5(θ), spanning from r_min = 0 to r_max = 10R_v. The tangential shear γ_t(r_p) is computed for each radial bin using Eq. (8). Finally, we constructed the stacked tangential shear profile by averaging the profiles of different voids in a chosen sample.

As done in Davies et al. (2019a) and Davies et al. (2021b), the uncertainty associated with the mean tangential shear profile of N void profiles in the i-th bin is evaluated as

$$\begin{array}{c}\hfill {\sigma }_i=\sqrt{\frac{C_{\mathit{ii}}}{N}},\end{array}$$(31)

where C_ii is the i-th diagonal element of the covariance matrix calculated using the tangential shear profiles extracted from all the 256 S/N maps. We refer to Appendix A for an in-depth analysis of the obtained covariance matrix.

In Fig. 9 we show the tunnel stacked tangential shear profiles for all the analyzed models. Firstly, we note that all tangential shear profiles are statistically negative, consistent with underdensities, and have an amplitude of the order of 10⁻³, matching values found in the literature for tunnel voids. In addition, we detect a γ_t(r_p) signal that remains negative out to the value of the radius assigned to our voids, resembling the behavior typically observed in WL minima-based finders (Davies et al. 2019a, 2021a,b). In particular, the observed trends are in agreement with theoretical expectations, namely, there are deeper tangential shear signals for MG models with stronger |f_R0| intensities and a shallower signal in the presence of massive neutrinos. This is explained by the fact that cosmic voids experience a depletion of their internal density profiles when their evolution is enhanced; conversely, they are less underdense when their growth is hampered.

Fig. 9. Top panel: Averaged tangential shear profiles extracted from 256 light-cone random realizations with redshift zs = 1, for the different analyzed models. Bottom panel: Corresponding residual with respect to the profile measured in the ΛCDM scenario. Both the plots utilize the following color scheme: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), and ΛCDM0.15 eV (purple). The shaded regions with the same colors around the curves represent the associated uncertainties, computed via the corrected covariance matrix (see Eqs. (31) and (A.2)). These are not displayed in the top panel for visualization reasons.

We note, however, that there is a pivot point at r_p/R_v ≃ 4 where all model results almost degenerate. Over this scale, there is an inversion in the trends of the profiles. This point marks the transition between regions influenced by underdensities and those affected by overdensities. Moving toward the outskirts of 2D voids, the tangential shear profiles begin to deviate from ΛCDM. This deviation follows the expected effects of overdensities, with an increased signal for MG models and a reduced signal for scenarios with massive neutrinos. The transition scale depends on the chosen definition of the 2D void radius. Although we tested several S/N thresholds in the range [ − 1, −3], we found that −2.5 offers the best trade-off between shear signal strength and void sample purity. Increasing the threshold shifts the pivot point to smaller r/R_v values, while lowering it moves the crossing outward. However, no clear predictive trend was observed. Higher thresholds also increase the incidence of 2D void mergers, reducing the void count and affecting the reliability of cosmological constraints. Therefore, using values near zero proved unfeasible.

Additionally, we focused on studying the dependency of the stacked tangential shear profile on the void size. This analysis allowed us to examine how the tangential shear properties vary with void radii, while providing valuable insights into the matter distribution around voids at different scales.

For each cosmological model, we divided the entire sample of our voids according to their size, into three equipopulated bins:

• Small voids with radii R_v in the range [1, 3.97] arcmin;

• Medium voids with radii R_v in the range ]3.97, 4.85] arcmin;

• Large voids with radii R_v in the range ]4.85, 15] arcmin.

Then we extracted the void tangential shear profiles from these equipopulated subsamples. The result is represented in Fig. 10.

Fig. 10. Top panel: Averaged tangential shear profiles measured from samples of voids with different sizes in the standard ΛCDM cosmology: small (in light-blue), medium (green) and large (red). We report with a black dashed line the one computed using the whole void sample. Bottom panel: Residuals between the profiles of each radius selection and the stacked void profile relative to the full sample. The shaded regions with same colors around the curves represent the uncertainty associated with the measures, analogously to Fig. 9.

From the top panel of this figure, we can identify a trend between the three bins. In particular, we note that the signal of small voids is deeper and rises more rapidly towards zero, while for the subsample of large voids, we observe a less profound signal and a slower rise. This behavior is similar to the one found for the density profiles of 3D voids (Nadathur et al. 2015; Hamaus et al. 2014; Voivodic et al. 2020), for which small voids show, in fact, deeper interiors and high compensation walls. This occurs because large voids have internal substructures and irregular shapes due to merger events. Furthermore, as the VSF analysis shows, small voids are much more numerous than large ones. As a result, the bin for large voids covers a broader range of radii. This causes the stacked signal for large 2D voids to be influenced by the varying tangential shear profiles, leading to a smoother average trend.

The bottom panel of Fig. 10 shows the residuals for the three cases, computed with respect to the tangential shear profile of the complete sample of voids, $\overline{{\gamma }_t}(r_p)$, which serves as our reference. We finally observe that the subsample of voids of intermediate size is the most similar to the average tangential shear profile computed with the totality of the void sample.

In Fig. 11, we present the stacked tangential shear profiles, divided into the three equipopulated bins for each cosmological model. This approach highlights the significant variations in the profile trends across the different models analyzed. It reveals a systematic trend in the amplitude and shape of the tangential shear profiles, with smaller voids (blue) exhibiting deeper signals compared to larger voids (red). This behavior reflects the distinct structural properties of voids at different scales and also in the alternative cosmological models analyzed. Furthermore, this observed dependency in void size and cosmological model will serve as a key parameter in Sect. 6, where we test different parameterizations to describe the tangential shear signal from WL tunnel voids. We finally highlight that, in the analyzed scenarios, the variations caused by differences in void size are more pronounced than those caused by changes in the cosmological model when considering voids of similar radius.

Fig. 11. Stacked tangential shear profiles of all the analyzed cosmologies, divided into subsamples of voids with different sizes. The split of the void sample is the same as presented in Fig. 10 and is used to distinguish voids of small (blue), medium (green), and large (red) sizes. The results for different cosmologies are shown using different line styles, as reported in the legend.

6. Modeling tangential shear profile of tunnel voids

In the final part of this work, we aim to model the tunnel void tangential shear profile using a parametric formula. Although no theoretical model based on first principles exists for the WL signal from voids, an empirical function can be used to fit the data to analyze the behavior and limits of its free parameters in different scenarios. To achieve this, we perform a Bayesian analysis by running a full Markov chain Monte Carlo (MCMC) procedure to sample the posterior distribution of the free parameters in the considered models. This approach does not directly constrain the cosmological parameters but represents a first step to understanding the physics behind the WL phenomenon in voids in order to test GR and MG models.

6.1. Tunnel void tangential shear from density profiles

The first method we followed is to start from a known parametric formula from the literature to represent the density profile of 3D voids. Following the approach of Pisani et al. (2014), Fang et al. (2019) and Boschetti et al. (2024), we can integrate the functional form of the density profile along the line of sight to obtain the projected surface mass density, Σ(r_p), which we can then relate to convergence through Eq. (21) and derive the tangential shear profile, γ(r_p).

The first functional form we consider to represent the density profile of 3D voids is the Hamaus-Sutter-Wandelt (HSW, Hamaus et al. 2014). It takes the form

$$\begin{array}{c}\hfill \frac{{\rho }_{\mathrm{v}}(r)}{\overline{\rho }}-1={\delta }_c\frac{1-{(r/r_{\mathrm{s}})}^{\alpha }}{1+{(r/r_{\mathrm{v}})}^{\beta }}.\end{array}$$(32)

It is characterized by four free parameters (δ_c, α, β, r_s), which represent, respectively, the central density contrast, the inner and outer slopes of the profile, and the characteristic scale where ${\rho }_{\mathrm{v}}=\overline{\rho }$. This functional form is commonly used to fit voids of different sizes and internal densities.

The second functional form we consider is the hyperbolic tangent (HT, Voivodic et al. 2020). It has only two free parameters, δ_c and s, and is expressed as:

$$\begin{array}{c}\hfill {\rho }_v(r/R_{\mathrm{v}})=1+|{\delta }_c|\overline{\rho }\{\frac{1}{2}[1+tanh\left(\frac{y-y_0}{s\left(R_{\mathrm{v}}\right)}\right)]-1\}.\end{array}$$(33)

In this parametrization, δ_c is the density contrast at the void center, y = ln(r/R_v) and y₀ = ln(r₀/r_v). The radius r₀ is determined by requiring that the integral of the profile up to R_v be Δ_v = 0.2. This allows us to express r₀(s) (in units of h⁻¹ Mpc ) as a second-order polynomial function: r₀(s) = 0.37s² + 0.25s + 0.89, where s represents the gradient of the profile. The parameter s works similarly to the concentration parameter in the NFW profile, governing the rate at which density increases as we move away from the center of the void.

As shown in Fig. 12, neither of these functions leads to a good fit of our data. Despite the large degrees of freedom of the considered parametric forms and the very large prior assigned to their coefficients, the complex shape of the tangential shear profiles could not be reproduced accurately. For both cases, the reduced χ² obtained with the best-fit parameters is indeed very far from unity, namely, ${\stackrel{\sim }{\chi }}_{\mathrm{HSW}}^2\simeq 2031$ and ${\stackrel{\sim }{\chi }}_{\mathrm{HT}}^2\simeq 7750$.

Fig. 12. Best fits of the averaged stacked tunnel voids tangential shear profiles, γt(rp), measured in ΛCDM cosmology adopting three different parametric functions. The data are derived as the mean of the 256 ΛCDM mock light cones and represented with black markers having error bars computed as specified in Sect. 5.2. The first two models considered are derived from the line-of-sight integration of two 3D void density profile functions: HSW (red dashed line) and HT (green dotted line). The last one is the new formula developed in this work for WL tunnel voids (blue solid line), reported in Eq. (34). The shaded area, in some cases barely visible, represents the 1σ uncertainty associated with each fit.

The motivation for this inconsistency must be sought from the very nature of our voids. In our case, we are trying to model the WL signal generated by tunnel voids, which are the result of projecting numerous underdensities along the line of sight. Following the approach described in this section, we are instead imposing the modeling of the tangential shear through the projection of the typical density profile of isolated 3D voids. This assumption cannot be valid for our tunnel voids, which instead derive from the projection of a complex distribution of underdensities with different sizes and positions. Therefore, we emphasize that in the analyses of Fang et al. (2019) and Boschetti et al. (2024), the usage of the 3D density profiles was effective for the following reasons. Both authors make use of a tomographic approach, considering only the WL signal generated by the matter distribution present in relatively thin redshift slices. In this way, they exploit a lens plane that follows the simplified condition of having only one, or in any case a few, 3D voids along the line of sight. In a future study, we will analyze the connection between the average tangential shear of individual 3D voids and that of tunnel voids identified with our finder using a tomographic approach.

6.2. A new parametric formula

Here, we present and validate a new parametric function suitable for the modeling of the tunnel voids’ tangential shear profiles extracted from our void catalogs. Our goal is to accurately reproduce the shape and amplitude of our shear profiles, employing the smallest possible number of free parameters. Moreover, we want this new functional form to be flexible enough to represent the shear profiles both for different void sizes and the five cosmological models analyzed in this work. A useful example to understand the degree of “flexibility” we need to reach with this function is shown in Fig. 11. To the best of the authors’ knowledge, this is the first and so far only parametric formula proposed in the literature specifically developed to model the tangential shear signal of WL tunnel voids.

The functional form found that best met these requirements can be expressed as

$$\begin{array}{c}\hfill {\gamma }_t(r_p)=a·r_p(\frac{1-r_p^b}{1+r_p^c}-\frac{exp(d·r_p)}{1+exp(e·r_p-(d+e))}).\end{array}$$(34)

The coefficients a, b, c, d, and e represent the five free parameters of our model. For now, these parameters are not intended to have a physical meaning but rather to regulate the shape and amplitude of the profile across different scales. Specifically, the parameter a is the normalization of the curve, while b shifts the position of its minimum. The steepness of the rise in the outer part of the profile is modulated by c, whereas d controls how deep the profile declines at its minimum. Finally, e defines the scale at which the exponential growth begins to dominate. To understand in detail their role, in Appendix B, we illustrate and discuss the effect that varying singularly each parameter produces on the shear profile.

In Fig. 12, we present one of the most important results of this work, namely the application of the new function for fitting the tunnel voids stacked tangential shear profile extracted from the ΛCDM simulations. It is easy to notice that there is an excellent match with the data across all scales. The reduced χ² is indeed now closer to the unity, namely, ${\stackrel{\sim }{\chi }}_{\mathrm{M}25}^2\simeq 0.37$. The low chi-squared value might suggest that the model could have more free parameters than necessary. However, it is important to note that the same model will also be used to fit the shear profile for voids selected by size and in MG scenarios, which shows a large and complex variation at all scales analyzed. Moreover, the reduced χ² obtained neglecting the off-diagonal terms is ${\stackrel{\sim }{\chi }}_{\mathrm{M}25}^2\simeq 4.44$, indicating a strong correlation between the radial bins (see Appendix A for details).

As previously shown in Fig. 10, there is a high dependency of the tangential shear profiles on the void size. According to the selected void radius range, the depth and the position of the profile minimum change, as well as the slope and the height of the outer part. The parametric formula proposed in Eq. (34) was built with enough degrees of freedom to capture these variations. This is analyzed in Appendix C, where we provide the best-fit values of the coefficients for the ΛCDM scenario, splitting the void sample into the same size bins shown in Fig. 10.

Now, we want to focus on the cosmological discriminating power of the tunnel void tangential shear profiles, and to do this, we will examine the posterior distribution of the parameters of the function in Eq. (34). Notice that, in this analysis, firstly we make no distinctions based on the size of the voids and perform the stacking of all the voids extracted for a given cosmology.

We report in Fig. 13 the confidence contours obtained from fitting the shear profiles extracted in simulations featuring the five cosmological models considered in this work: ΛCDM, fR4, fR5, fR6, ΛCDM_0.15 eV. Overall, we note that the confidence contours are mostly ordered according to the level of enhancement or damping of the growth of LSS, with fR4 and ΛCDM_0.15 eV showing the largest separation from the ΛCDM case, but in opposite directions. We note, however, that the confidence contours result in partially overlapping, especially when considering the models characterized by the smaller deviations from the ΛCDM case. The area of the confidence contours is in fact significantly enlarged by the correlations between radial bins, which decrease the statistical disentangling between the scenarios. Nevertheless, the separation between the most extreme cases – fR4 and ΛCDM_0.15 eV – from the reference ΛCDM is over 1σ for all the parameters of the model.

Fig. 13. 68% and 95% confidence levels for the five parameters of the functional form reported in Eq. (34), used to model the tunnel void tangential shear profiles extracted from simulations featuring different cosmological scenarios: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), ΛCDM0.15 eV (purple).

It is also apparent from this plot that notable positive correlations are present in the (a, c) and (d, e) planes, while negative correlations are observed in the (a, e) and (a, d) ones. The latter, in particular, is characterized by tilted and very elongated ellipses, indicating a strong correlation between these parameters. We refer to Appendix B for a deeper understanding of the behavior of the model coefficients. In Fig. 14, we present the same confidence contours derived for all analyzed cosmological models, but restricted to the sub-sample of tunnel voids with R_v in the range ]3.79, 4.85] arcmin.

Fig. 14. Same as Fig. 13, but for a selection of voids with radius in the range ]3.79, 4.85] arcmin.

The effect of the void radius selection is not very strong, showing the robustness of the methodology. However, we note that a particular selection in size may help the distinction between specific scenarios. For example, the constraints computed for the ΛCDM_0.15 eV case appear now more detached from the reference ΛCDM. This might suggest that selecting certain spatial scales could enhance the effects caused by the presence of massive neutrinos, making it more favorable for studying cosmological scenarios with these characteristics. Additionally, selecting by radius prevents the averaging of very different profile behaviors. This implies that WL analyses of voids can be further optimized by carefully selecting void size ranges to enhance sensitivity to both modifications of gravity and neutrino mass effects.

For the sake of conciseness, we do not show the results for the other two radius selections (i.e., small and large voids). However, we clarify that the results are consistent with those presented for medium-sized voids. Finally, we emphasize once again that these filters in void size can be used to leverage specific differences between the analyzed cosmological models. For example, we found that applying the selection for large voids (]4.85, 15] arcmin) makes it possible to highlight the differences between the fR6 and ΛCDM cases.

We finally underline that the confidence contours presented in this section are meant to be representative of a photometric survey like the one currently being carried out by Euclid (see Sects. 3.2 and 3.3). In our case, however, each of our 256 maps has a sky area of 25 deg². Such a configuration can be interpreted as a mosaic of 256 semi-independent observations, leading to an effective sky area not greater than 6400 deg², which is more than a factor of 2 smaller than the one planned for the final data release of Euclid (i.e., 14 000 deg²). We believe, therefore, that the constraining power showcased by this analysis to be close to one half (45.7%) of the one expected for the Euclid wide survey.

To summarize this section, we can conclude that our identified tunnel voids are sensitive to the effects of MG and massive neutrinos. Furthermore, the modeling of their tangential shear profiles via the newly proposed parametric formula can offer new insights to analyze the variations induced by these cosmological models.

7. Summary and conclusions

In this work, we investigated the WL signal of tunnel voids in the convergence field, focusing on MG models within the f(R) framework (Hu & Sawicki 2007) and their effects with respect to those of massive neutrinos. These alternative f(R) models are considered viable as they reproduce both early- and late-time cosmic acceleration without a cosmological constant, while employing screening mechanisms to remain consistent with GR on small scales and in high-density regions. We also included a massive neutrino model to disentangle its degenerate effects from those of MG on LSS formation.

Previous studies have demonstrated that cosmic voids are powerful cosmological probes for testing deviations from GR, particularly in the context of f(R) models. For instance, 3D void statistics have been used to constrain f(R) gravity and massive neutrino scenarios (Contarini et al. 2021), demonstrating their ability to break degeneracies between these components. Similarly, WL statistics of both 2D and tunnel voids have exhibited enhanced sensitivity to f(R)-induced modifications of structure formation in lensing surveys similar to Euclid and LSST (Cautun et al. 2018).

Using an innovative pipeline, we analyzed the S/N of noised and smoothed convergence maps from the DUSTGRAIN-pathfinderN-body simulations (Giocoli et al. 2018a; Peel et al. 2018), systematically comparing their behavior with the standard ΛCDM model. We summarize our main results below.

• We presented a 2D WL tunnel void finder called PyTwinPeaks that allows computationally efficient identification of underdense regions in convergence maps, providing precise measurements of void centers, sizes, and geometries without requiring 3D galaxy positions.

• We studied the void size function, which revealed that models with stronger MG effects display an increased number of tunnels, while the presence of massive neutrinos suppresses their abundance.

• We analyzed the stacked tangential shear profile of identified tunnel voids. We extracted the average stacked tangential shear profile across the five cosmological models, finding that voids in f(R) MG scenarios exhibit a more pronounced shear signal owing to their deeper density contrasts. Conversely, voids in ΛCDM_0.15 eV display a weaker signal due to the suppression of LSS growth.

• We investigated the behavior of the tangential shear profile according to the size of voids using three equipopulated subsamples. Smaller voids show a deeper and steeper signal, while larger voids have a shallower and flatter signal.

• We modeled the void shear profiles, developing a new five-parameter empirical parametric formula. Importantly, it fits all the stacked shear profiles accurately, making it a valuable tool for analyzing WL signals in tunnel voids.

It is important to underline that the differences in the shear measurements can be used to discriminate between various cosmological models, revealing the presence of MG effects. This is especially important at the minimum of the profile, where this approach enables a clear statistical separation. A Bayesian analysis was performed to constrain the free parameters of the proposed formula. We first applied MCMC sampling to the ΛCDM shear profiles, analyzing the corresponding confidence contours to examine parameter correlations. Scrutinizing different cosmologies, we observed that confidence contours shift according to structure growth rates, with fR4 and ΛCDM_0.15 eV displaying the most extreme behavior.

Given our encouraging findings, we plan to further expand this project in several ways. First, the current version of PyTwinPeaks, which we have already made publicly available, will be restructured to handle both WL underdensities and overdensities. To strengthen the statistical robustness of our results, we also plan to apply the pipeline to larger-volume simulations (Euclid Collaboration: Castander et al. 2025) and explore a broader range of cosmological models (Euclid Collaboration: Ajani et al. 2023).

This development is supported by the modular structure of our pipeline, where each component can be independently adapted, making the methodology applicable across different cosmological suites of mock catalogs. The core of the method is the 2D void finder, which can be easily configured to process input maps produced with various lensing reconstruction techniques, including full-sky multiplane ray-tracing or tomographic approaches. Additionally, shear-based lensing estimators can replace convergence to improve compatibility with real data affected by masking or anisotropic noise in κ reconstruction. Moreover, the new formula for WL tunnel voids is expected to generalize across a wide range of scenarios.

An extended WL tomographic analysis will be performed to track signal variations across redshifts and identify the optimal range for studying void tangential shear profiles. By following the evolution of cosmic structures in redshift bins, we aim to extract higher-order WL statistics that capture key aspects of structure formation and growth. This will support forecasts for Stage-IV LSS surveys and enable tests of dynamical dark energy and MG. Ultimately, combining WL statistics from both voids and clusters will yield tight cosmological constraints and offer a unified framework for probing fundamental questions in cosmology.

Furthermore, we will explore the feasibility of applying our analysis to real data catalogs. Our new tunnel void pipeline will be applied to the first data release of the ESA-Euclid mission, allowing for a direct observational test of our methodology. The integration of standard cosmological probes alongside WL measurements from voids will maximize the information extracted, leveraging the complementarity of constraints from under- and overdensities of LSS (Bayer et al. 2021; Kreisch et al. 2022; Contarini et al. 2024; Pelliciari et al. 2023). The photometric galaxy catalog from Euclid will allow us also to analyze the cross-correlation between WL and underdense regions, namely, the void-lensing cross-correlation (see e.g., Bonici et al. 2023).

Finally, our long-term objective is to develop a theoretical model for the tangential shear profile of cosmic voids from first principles, a step that remains an open challenge. This will require exploring semi-analytical models to provide a more comprehensive framework for using WL tunnel voids and 3D void lensing as a cosmological probe. By deriving a physically motivated model, we aim to directly constrain cosmological parameters and their evolution over time, improving the predictive power of void-lensing studies.

Data availability

To our knowledge, no 2D void finder has been publicly released to date. To support reproducibility and further development, we provide our WL tunnel void finder implementation as open-source code at https://github.com/LeonardoMaggiore/PyTwinPeaks.

Acknowledgments

We thank the anonymous referee for his/her comments that helped improve the presentation of our results. We would like to thank A. Pisani, N. Hamaus, N. Schuster, and C.T. Davies for their help, useful discussions, and fruitful collaborations. We also acknowledge the use of computational resources from the parallel computing cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at the Physics and Astronomy Department in Bologna. LM and CG acknowledge the financial contribution from the PRIN-MUR 2022 20227RNLY3 grant ’The concordance cosmological model: stress-tests with galaxy clusters’ supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 “Attività scientifiche per la missione Euclid – fase E”. GC thanks the support from INAF theory Grant 2022: Illuminating Dark Matter using Weak Lensing by Cluster Satellites. We acknowledge the use of the Python libraries NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), COBAYA (Torrado & Lewis 2021) and GETDIST (Lewis 2019).

References

Appendix A: Covariance matrix

In this appendix, we present the tangential shear covariance matrix for the WL tunnel void sample of ΛCDM scenario. The diagonal of this matrix is used to calculate the errors in the tunnel voids stacked tangential shear profile through Eq. (31), while its inverse is inserted in the likelihood to compute the fit of the data through Eq. (34).

Figure A.1 shows the measured covariance matrix, C_ij, for WL tunnel voids identified through the developed finder, in WL maps with GSN included and with a smoothing scale of θ_s = 2.5 arcmin. It is computed via the formula

$$\begin{array}{c}\hfill C_{\mathit{ij}}=\frac{1}{N-1}{\displaystyle \sum _{k=1}^N}[{\gamma }_{\mathrm{t}}(i)-{\overline{\gamma }}_{\mathrm{t}}(i)][{\gamma }_{\mathrm{t}}(j)-{\overline{\gamma }}_{\mathrm{t}}(j)],\end{array}$$(A.1)

where i and j are radial bin indices and N = 22203 is the total number of samples of WL tunnel voids in the ΛCDM scenario. γ_t is the tangential shear and an over-bar denotes the mean from N WL tunnel voids. Since we considered our light cones to be independent realizations, N is not the number of our maps (as in Davies et al. 2021b) but the number of voids used to calculate the covariance matrix. The size of this matrix is of 64 × 64 and represents the statistical relationship between the radial bins used to calculate the stacked tangential shear profile.

For the calculation of the log-likelihood used in the Bayesian analysis presented in Sect. 6.2, we used the precision matrix obtained by inverting the tangential shear covariance matrix. This inverse was corrected using the Anderson-Hartlap factor (Anderson 2003; Hartlap et al. 2007), α, through the relation $C_{\text{corr}}^{-1}=\alpha C^{-1}$, where C⁻¹ is the inverse of the sample covariance matrix. The α factor is given by

$$\begin{array}{c}\hfill \alpha =\frac{N-N_{\mathrm{bin}}-2}{N-1},\end{array}$$(A.2)

where N is the total number of voids in the sample, and N_bin = 64 is the number of radial bins. We use this method⁶ to correct for the bias introduced by inverting a covariance matrix of profiles extracted from noised maps.

The covariance matrix shown in Fig. A.1 exhibits a dominant diagonal structure, indicating a strong correlation between the adjacent bins. This effect is caused by the overlap between the void shear profiles of the same S/N map, which occurs especially at larger distances from the center. Covariance matrices computed from alternative cosmologies, while remaining of the same order of magnitude, exhibit variations at the 10⁻⁶ level relative to the ΛCDM case, reflecting changes in LSS.

The smooth gradient observed in the off-diagonal elements highlights the progressive nature of correlations, with stronger covariance at large separations. This suggests that the uncertainty in the tangential shear profile does not behave as purely uncorrelated Gaussian noise but instead encodes features intrinsic to the void environment. Additionally, the overall amplitude of the covariance matrix is a direct indication of the sample variance in the WL tunnels measurement, implying that increasing the number of detected voids would reduce uncertainties proportionally.

Fig. A.1. Tunnel voids tangential shear profiles covariance matrix measured for the entire sample of voids in the ΛCDM scenario, in each of the 64 radial bins.

A key observation is that correlations remain relatively strong for large radial distances, where the signal transitions into the background noise. The transition from positive to negative correlations in off-diagonal elements also suggests complex interactions between shear contributions from neighboring voids, which could be further investigated in future studies.

Overall, these matrices provide critical insights into the statistical properties of the measured WL tunnel void signal. The observed correlations highlight the necessity of including the full covariance structure in parameter estimation and model fitting, as neglecting these correlations could lead to biased constraints on cosmological parameters. Future improvements in the void detection algorithm and larger survey volumes, such as those expected from Euclid, will be crucial in reducing the sample variance and enhancing the robustness of WL tunnel voids and overdense structures as a cosmological probe.

Appendix B: Free parameters of the model

Here, we focus on the interpretation of the coefficients of Eq. (34), providing a brief description of the main effects of each parameter:

• a)

  Amplitude or normalization. It modifies the amplitude of the function, i.e., its normalization, and it is related to the overall depth of the tunnel void.

• b)

  Position of the minimum. It plays a role in determining the position of the absolute minimum of the profile, also influencing the amplitude of the profile at large distances from the void center.

• c)

  Exponential growth. It regulates the exponential rise after the minimum, so the slope of the outer part of the profile.

• d)

  Depth of the minimum. It determines the minimum value of the profile, and it is related to the central density contrast of the tunnel void.

• e)

  Starting point of the exponential term. It influences the scale at which the exponential part of the function begins to become dominant.

Figure B.1 illustrates the influence of each parameter on the tangential shear profile. Variations in parameter a primarily rescale the entire function, while parameter b shifts the location of the minimum. Parameters c and e control the steepness of the transition at large scales, whereas d directly affects the depth of the void-induced shear signal.

Fig. B.1. Effects of varying individually the coefficients of the function reported in Eq. (34). For each parameter, ten curves are represented and colored according to the value of the selected parameter. The colorbar located on the right of each subplot shows the association of the color used with the parameter value. The different panels refer to the parameters a (top left), b (top right), c (central left), d (central right), and e (bottom left).

Regarding the physical meaning of the parameters, we can currently only hypothesize that it is linked to the complex morphology and density profile of tunnel voids, particularly around the peak of the lensing kernel that marks the segment of the line of sight where the density distribution carries the most weight. (see Sect. 2.1).

The interplay between these parameters provides a novel approach to characterizing cosmic voids in WL maps. By fitting void shear profiles across different cosmological models, we can extract information on the underlying physics governing LSS formation and evolution. The response of these parameters to varying gravitational models and neutrino masses enables a deeper understanding of how modifications to gravity and the presence of hot dark matter affect the large-scale distribution of mass in the Universe.

The results obtained through this parametric modeling demonstrate that WL tunnel voids offer a robust tool for probing the nature of gravity on cosmic scales. By refining the interpretation of these parameters and linking them to fundamental cosmological quantities, future analyses can enhance the constraining power of void lensing studies and strengthen the role of cosmic voids as a key probe for testing deviations from General Relativity. Extensions of this analysis could explore whether these constraints tighten further with increased void statistics or if systematic uncertainties, such as shape noise or projection effects, influence the parameter degeneracies.

Appendix C: Void size dependency

Here we focus on the dependence of the best-fit parameters of Eq. (12) on the size of the tunnel voids. The data used for the fit are the averaged void shear profiles extracted from the 256 realizations of the ΛCDM S/N maps, splitting the void sample into the three equipopulated subsamples of voids described in Sect. 5. We report in Table C.1 the best-fit values of each of the five free parameters of the model, together with their associated 1σ uncertainty. The best fit of our model is obtained by extracting the median value of the individual parameter distributions, which are well approximated by a Gaussian.

In this case, we decided not to represent the confidence contours in a corner plot because the size of these contours would have been too small. In fact, the error associated with these best fits is, in most cases, between 0.4% and 2.7%, with an average relative uncertainty around 2%. The only exception is the parameter d, which for the medium-size tunnel voids shows a relative uncertainty of approximately 9.5%. The parameter c appears instead to be the most constrained one, with a percentage error ranging from 0.43% to 0.78%. We believe that the parameter values shown in Table C.1 can be used as a reference for future studies on tunnel voids.

The degeneracy between these parameters is also dependent on the void size, although to an almost negligible extent. We refer to Figs. 13 and 14 for a more in-depth study of the covariance between these parameters with different void size selections.

All Tables

All Figures

Fig. 1. Example of light-cone construction with MAPSIM from z = 0 to z = 4 for a ΛCDM simulation with a box size of 750 h−1 Mpc, tiled to encompass the 5 × 5 deg2 light cone. The light cone geometry is depicted by the red solid line, expanding until its transverse comoving size matches that of the simulation box. Different colors represent different boxes, while the total of 21 different snapshots are defined by the dashed lines. The source plane (blue solid line) for ray tracing is placed at zs = 1. The snapshots within zobs < z < zs are the 13 utilized ones, while the gray region indicates the snapshots with z > zs that are not considered in this study. The yellow solid lines represent the 27 lens planes used, onto which the particles are projected. Among these, only those within the considered snapshots are utilized for constructing the convergence maps used in this work, positioned from z = 0.05 to 0.98, respectively.

In the text

Fig. 2. Example of a realization of the ΛCDM light cone, having a FOV of 5 × 5 deg2 aperture and zs = 1. In the upper panel, we show the WL convergence map, and in the lower panel, the S/N of the same map with smoothing scale at θG = 2.5 arcmin. The color scale in both panels represents the convergence value for each pixel in the top panel and the S/N value in the bottom panel, as indicated by the corresponding color bars.

In the text

Fig. 3. Top panel: Average PDFs of the noised and smoothed convergence field from 256 light-cone random realizations with redshift zs = 1, for all the analyzed models: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red) and ΛCDM0.15 eV (purple). The black dashed line represents the average PDF without the addiction of GSN for the $\mathrm{\Lambda }\mathrm{CDM}$ model. Bottom panel: PDF noised and smoothed corresponding residuals computed with respect to the ΛCDM model. For the purposes of visualization, only the uncertainties of the noiseless and the reference ΛCDM measurements are reported in this plot as hatched gray and shaded gray areas, respectively.

In the text

Fig. 4. Average angular power spectra of the noised and smoothed convergence field from 256 light-cone random realizations with redshift zs = 1, for all the analyzed models (top panel) and the correspondent residuals with respect to the reference ΛCDM model (bottom panel). Line styles and color schemes are the same as used in Fig. 3. The dotted vertical line marks the angular mode of the characteristic smoothing scale in the Fourier transform of the κ field.

In the text

	Fig. 5. Diagram showing the four main steps of the 2D void finder algorithm that allows the identification of voids from a WL convergence map.
In the text

Fig. 6. Left panel: Graphical representation of all 2D voids with their preliminary radii of each threshold produced by the same map shown in Fig. 2, centered at the S/N-weighted centroids of the connected regions. Right panel: Same as the left panel, but with voids centered at the local minima of the connected regions. In both, each void of a given threshold is represented with the colored area of the same color, this color varies depending on the S/N threshold as shown in the colorbar.

In the text

	Fig. 7. Example of the final tunnels catalog associated with the same S/N map of Fig. 2 and Fig. 6. Colors from blue to red represent areas with increasing S/N, as indicated by the colorbar on the right. The 2D voids identified by our algorithm are shown as areas delimited by black circles. The corresponding centers (absolute minima) are marked with yellow crosses.
In the text

Fig. 8. Top panel: Average void size functions, per degree squared, of 2D WL voids from 256 light-cone random realizations with redshift zs = 1, for the different analyzed models: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), and ΛCDM0.15 eV (purple). Bottom panel: Corresponding residuals computed with respect to the ΛCDM model. In both panels, the shaded regions of the same colors around the curves indicate errors, computed as the standard deviation of counts in each bin divided by the number of maps.

In the text

Fig. 9. Top panel: Averaged tangential shear profiles extracted from 256 light-cone random realizations with redshift zs = 1, for the different analyzed models. Bottom panel: Corresponding residual with respect to the profile measured in the ΛCDM scenario. Both the plots utilize the following color scheme: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), and ΛCDM0.15 eV (purple). The shaded regions with the same colors around the curves represent the associated uncertainties, computed via the corrected covariance matrix (see Eqs. (31) and (A.2)). These are not displayed in the top panel for visualization reasons.

In the text

Fig. 10. Top panel: Averaged tangential shear profiles measured from samples of voids with different sizes in the standard ΛCDM cosmology: small (in light-blue), medium (green) and large (red). We report with a black dashed line the one computed using the whole void sample. Bottom panel: Residuals between the profiles of each radius selection and the stacked void profile relative to the full sample. The shaded regions with same colors around the curves represent the uncertainty associated with the measures, analogously to Fig. 9.

In the text

	Fig. 11. Stacked tangential shear profiles of all the analyzed cosmologies, divided into subsamples of voids with different sizes. The split of the void sample is the same as presented in Fig. 10 and is used to distinguish voids of small (blue), medium (green), and large (red) sizes. The results for different cosmologies are shown using different line styles, as reported in the legend.
In the text

Fig. 12. Best fits of the averaged stacked tunnel voids tangential shear profiles, γt(rp), measured in ΛCDM cosmology adopting three different parametric functions. The data are derived as the mean of the 256 ΛCDM mock light cones and represented with black markers having error bars computed as specified in Sect. 5.2. The first two models considered are derived from the line-of-sight integration of two 3D void density profile functions: HSW (red dashed line) and HT (green dotted line). The last one is the new formula developed in this work for WL tunnel voids (blue solid line), reported in Eq. (34). The shaded area, in some cases barely visible, represents the 1σ uncertainty associated with each fit.

In the text

	Fig. 13. 68% and 95% confidence levels for the five parameters of the functional form reported in Eq. (34), used to model the tunnel void tangential shear profiles extracted from simulations featuring different cosmological scenarios: ΛCDM (black), fR4 (green), fR5 (blue), fR6 (red), ΛCDM0.15 eV (purple).
In the text

	Fig. 14. Same as Fig. 13, but for a selection of voids with radius in the range ]3.79, 4.85] arcmin.
In the text

	Fig. A.1. Tunnel voids tangential shear profiles covariance matrix measured for the entire sample of voids in the ΛCDM scenario, in each of the 64 radial bins.
In the text

Fig. B.1. Effects of varying individually the coefficients of the function reported in Eq. (34). For each parameter, ten curves are represented and colored according to the value of the selected parameter. The colorbar located on the right of each subplot shows the association of the color used with the parameter value. The different panels refer to the parameters a (top left), b (top right), c (central left), d (central right), and e (bottom left).

In the text