Weak-lensing tunnel voids in simulated light-cones: a new pipeline to investigate modified gravity and massive neutrinos signatures

Cosmic voids offer a unique opportunity to explore modified gravity (MG) models. Their low-density nature and vast extent make them especially sensitive to cosmological scenarios of the class $f(R)$, which incorporate screening mechanisms in dense, compact regions. Weak lensing (WL) by voids, in particular, provides a direct probe for testing MG scenarios. While traditional voids are identified from 3D galaxy positions, 2D voids detected in WL maps trace underdense regions along the line of sight and are sensitive to unbiased matter distribution. To investigate this, we developed an efficient pipeline for identifying and analyzing tunnel voids, i.e., underdensities detected in WL maps, specifically in the signal-to-noise ratio (SNR) of the convergence. In this work, we used this pipeline to generate realistic SNR maps from cosmological simulations featuring different $f(R)$ scenarios and massive neutrinos, comparing their effects against the standard $\Lambda$CDM model. Using the convergence maps and the 2D void catalogs, we analyzed various statistics, including the probability density function, angular power spectrum, and void size function. We then focused on the tangential shear profile around 2D voids, demonstrating how the proposed void-finding algorithm maximizes the signal. We showed that MG leads to deeper void shear profiles due to the enhanced evolution of cosmic structures, while massive neutrinos have the opposite effect. Furthermore, we found that parametric functions typically applied to 3D void density profiles are not suitable for deriving the shear profiles of tunnel voids. We, therefore, proposed a new parametric formula that provides an excellent fit to the void shear profiles across different void sizes and cosmological models. Finally, we tested the sensitivity of the free parameters of this new formula to the cosmological model, revealing its potential as a probe for detecting the effects of MG models and the presence of massive neutrinos.

1 Introduction

Despite the success of the $\Lambda$-cold dark matter ($\Lambda$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 tensions between observations remain. In particular, discrepancies exist between Cosmic microwave background (CMB) (Planck Collaboration et al. 2020a, b) constraints and both lensing measurements of $\sigma_{8}$ (Abdalla et al. 2022; Di Valentino et al. 2021a) and supernovae determinations of $H_{0}$ (Di Valentino et al. 2021b), motivating 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 $\Lambda$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). WL introduces distinct distortions in background galaxies, quantified by cosmic shear $\gamma$ and convergence $\kappa$ (Bartelmann & Schneider 2001; Kilbinger 2015; Ishak 2019; Umetsu 2020) and it 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. 2024), in particular, will increase WL source counts by a $4-5$ factor, extending the redshift distribution (Kitching & Heavens 2017), with respect to ground-based surveys like KiDS (Hildebrandt et al. 2017; Energy Survey et al. 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 ($\gamma_{t}$) or convergence ($\kappa$), 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 Mpc (Higuchi et al. 2013; Gruen et al. 2015; 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 $\leq 50\ h^{-1}$ Mpc, where the thin lens approximation holds, and corresponds to studying the WL signal of voids using a tomographic approach (see Boschetti et al. 2023, 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 (SNR) of $\gamma_{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, and we apply it 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 adopt two different approaches to model the tangential shear signal; first, we integrate two established 3D void surface density profiles along the line of sight, and then, we introduce 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 weak lensing 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 $\Phi$ and $\Psi$ as the Einstein-frame metric potentials, the perturbed metric is described as:

$$\mathrm{d}s^{2}=(1+2\Psi)\mathrm{d}t^{2}-a^{2}(t)(1-2\Phi)\mathrm{d}x_{i}\mathrm{\leavevmode\nobreak\ d}x_{j}\,.$$ (1)

Here, $\Psi$ represents the gravitational potential, which describes the temporal part of the metric perturbations, while $\Phi$ is the gravitational curvature potential, describing the spatial part of the metric perturbations.

Let $\bm{\beta}$ represent the true angular position of a source at a comoving line of sight distance $\chi_{\mathrm{s}}$ and redshift $z_{\mathrm{s}}=z(\chi_{\mathrm{s}})$, and $\bm{\theta}$ its observed angular position on the sky. In this framework, the angular positions $\bm{\beta}$ and $\bm{\theta}$ are connected via the lens equation, which traces the trajectory of the photon as it is deflected by an angle $\bm{\alpha}$ through a lens or lens system. It is generalized as:

$$\bm{\beta}(\bm{\theta},z_{\mathrm{s}})=\bm{\theta}-\bm{\alpha}=\bm{\theta}-\frac{2}{c^{2}}\int_{0}^{\chi_{\mathrm{s}}}\frac{f_{\mathrm{ls}}}{f_{\mathrm{l}}f_{\mathrm{s}}}\nabla_{\bm{\beta}}\Phi_{\mathrm{len}}\bigl{(}\bm{\beta}(\bm{\theta},\chi_{\mathrm{l}}),\chi_{\mathrm{l}},z_{\mathrm{l}}\bigr{)}\,\mathrm{d}\chi_{\mathrm{l}}\,,$$ (2)

where $c$ is the speed of light and $f_{\mathrm{ls}}=f_{K}(\chi_{\mathrm{s}}-\chi_{\mathrm{l}})$, $f_{\mathrm{l}}=f_{K}(\chi_{\mathrm{l}})$, $f_{\mathrm{s}}=f_{K}(\chi_{\mathrm{s}})$, where we assume the case $\chi_{\mathrm{s}}>\chi_{\mathrm{l}}$ with “s” and “l” representing the source and the lens, respectively. Here, we assume a generic curvature $K$, so $f_{K}(\chi)$ represents the comoving angular diameter distance, and $f_{K}(\chi)\bm{\theta}$ corresponds to the transverse comoving vector for a comoving line-of-sight distance $\chi=\chi(z)$ (Misner et al. 1973; Hogg 1999; Umetsu 2020). $\Phi_{\mathrm{len}}\bigl{(}\bm{\beta},\chi_{\mathrm{l}},z_{\mathrm{l}}\bigr{)}$ denotes the 3D lensing gravitational potential generated by a mass distribution $\rho(\bm{\beta},\chi_{\mathrm{l}})$ at redshift $z_{\mathrm{l}}$ (Hilbert et al. 2020).

Generally, $\Phi_{\rm len}=(\Phi+\Psi)/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:

$$\Phi=\Psi=\Psi_{\rm N}\,,$$ (3)

where $\Psi_{\rm 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 $\Phi_{\rm len}=\Phi$ under this assumption. $\Phi$ is related to the non-relativistic matter density contrast, $\delta_{m}=\delta\rho/\bar{\rho}$, through the Poisson equation:

$$\nabla^{2}\Phi\bigl{(}\bm{\beta},\chi_{\mathrm{l}},z_{\mathrm{l}}\bigr{)}=4\pi Ga^{2}\delta\rho(\bm{\theta},\chi_{\mathrm{l}},z_{\mathrm{l}})=\frac{3\Omega_{m}H_{0}^{2}}{2a}\delta_{m}(\bm{\theta},\chi_{\mathrm{l}},z_{\mathrm{l}})\,,$$ (4)

where $G$ is the gravitational constant, $\bar{\rho}$ is the mean matter density of the Universe, $\Omega_{m}$ is the total matter density parameter, and $H_{0}$ represents the Hubble parameter at the present time.

The image distortion is described by the Jacobian distortion matrix $\mathbf{A}(\bm{\theta},\chi_{\mathrm{s}})$, which represents the relative positions of nearby light rays and is obtained by differentiating the lens equation with respect to $\bm{\theta}$. Assuming the flat-sky approximation, where $\bm{\beta}$ and $\bm{\theta}$ are 2D Cartesian coordinate vectors (Becker 2013), we can write:

$$A_{ij}(\bm{\theta},\chi_{\mathrm{s}})\equiv\frac{\partial\beta_{i}(\bm{\theta},z_{\mathrm{s}})}{\partial\theta_{j}}=\delta_{ij}-\frac{\partial\alpha_{i}(\bm{\theta},z_{\mathrm{s}},z_{\mathrm{l}})}{\partial\theta_{j}}\,,$$ (5)

and the distortion matrix $\mathbf{A}(\bm{\theta},\chi_{\mathrm{s}})$ takes the form:

$$\mathbf{A}(\bm{\theta},\chi_{\mathrm{s}})\approx\begin{pmatrix}1-\kappa-\gamma_{1}&-\gamma_{2}-\omega\\ -\gamma_{2}+\omega&1-\kappa+\gamma_{1}\\ \end{pmatrix}\,,$$ (6)

where $\delta_{ij}$ is the Kronecker delta. This defines the lensing convergence $\kappa$, the lensing complex shear $\gamma=\gamma_{1}+\mathrm{i}\gamma_{2}$, and the lensing rotation $\omega$, which is $\mathcal{O}(\Phi^{2})$ (Petri et al. 2017). In the WL regime $|\kappa|,|\gamma|,|\omega|\ll 1$.

To solve the integral form of Eq. (5), a common approach is to expand it in a series of $\Phi$, 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 $\bm{\alpha}$ to be expressed as the gradient of a 2D lensing potential, $\psi_{\rm len}$. Under this framework, the WL signal arises from the cumulative matter distribution along the line of sight (Davies et al. 2021b). The convergence $\kappa$ follows from the 2D Poisson equation and can be written in terms of the matter density contrast $\delta_{\rm m}$ as:

$$\kappa_{\rm Born}(\bm{\theta},z_{\mathrm{s}})=\frac{3H_{0}^{2}\Omega_{m}}{2c^{2}}\int^{\chi_{\mathrm{s}}}_{0}(1+z_{\mathrm{l}})W(\chi_{\mathrm{l}},\chi_{\mathrm{s}})\delta_{\rm m}(\bm{\theta},\chi_{\mathrm{l}},z_{\mathrm{l}})\,\mathrm{d}\chi_{\mathrm{l}}\,,$$ (7)

where $W(\chi_{\mathrm{l}},\chi_{\mathrm{s}})\equiv(f_{\mathrm{l}}f_{\mathrm{ls}})/f_{\mathrm{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 $\kappa(r_{p})$ is related to its radial tangential shear profile through the equation (Davies et al. 2018, 2021b):

$$\gamma_{\rm{t}}(r_{p})=\bar{\kappa}(<r_{p})-\kappa(r_{p})\,,$$ (8)

where $\bar{\kappa}(<r_{p})$ is the mean enclosed convergence within radius $r_{p}$ defined as:

$$\bar{\kappa}(<r_{p})=\frac{1}{\pi r_{p}^{2}}\int_{0}^{r_{p}}2\pi r_{p}^{\prime}\kappa(r_{p}^{\prime})dr_{p}^{\prime}\,.$$ (9)

By 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 $\theta$.

2.2 $f(R)$ modified gravity models

In this work, we analyze 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 $\Lambda$ of the standard $\Lambda$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:

$$S=\frac{M_{\mathrm{P}}^{2}}{2}\int\mathrm{d}^{4}x\sqrt{-g}(R+f(R))+S_{m}\left[g_{\mu\nu},\psi\right],$$ (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 $\psi$. The function $f(R)$ controls deviations from GR, becoming significant in the low curvature regime ($R\rightarrow 0$) and in low-density environments such as voids (Brans & Dicke 1961). GR can be recovered by setting $f=-2\Lambda^{\mathrm{GR}}$ for a more generalized cosmic acceleration, see Carroll et al. (2004).

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

$$f(R)=-m^{2}\frac{c_{1}\left(\frac{R}{m^{2}}\right)^{n}}{c_{2}\left(\frac{R}{m^{2}}\right)^{n}+1},$$ (11)

where $m^{2}\equiv H_{0}^{2}\Omega_{\mathrm{m}}$ defines the mass scale, and $c_{1},c_{2},$ and $n$ are non-negative parameters. For $\Lambda$CDM consistency, the condition $c_{1}/c_{2}=6\Omega_{\Lambda}/\Omega_{\mathrm{m}}$ must hold. In the limit $c_{2}\left(R/m^{2}\right)^{n}\gg 1$, the scalar field $f_{R}\equiv\mathrm{d}f(R)/\mathrm{d}R$ approximates to:

$$f_{R}\approx-n\frac{c_{1}}{c_{2}^{2}}\left(\frac{m^{2}}{R}\right)^{n+1}.$$ (12)

For $n=1$, the model simplifies, with $f_{R0}$ defined as:

$$f_{R0}\equiv-\frac{1}{c_{2}}\frac{6\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}\left(\frac{m^{2}}{R_{0}}\right)^{2},$$ (13)

where $R_{0}$ is the background value of $R$ at the present time.

By modifying the action in Eq. (10) with respect to $g_{\mu\nu}$, we obtain the modified Einstein equations for $f_{R}$:

$$f_{R}R_{\mu\nu}-\frac{1}{2}fg_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}f_{R}+g_{\mu\nu}\Box f_{R}=8\pi GT_{\mu\nu}\,,$$ (14)

where $\nabla$ is the covariant derivative and $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the D’Alembert operator. From its trace, we derive the motion equation for the scalar field:

$$\nabla^{2}\delta f_{R}=\frac{a^{2}}{3}\left[\delta R\left(f_{R}\right)-8\pi G\delta\rho\right],$$ (15)

and by extracting its time-time component, assuming small perturbations $\delta f_{R}$, $\delta R$, $\delta\rho$ 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\psi=\delta g_{00}/g_{00}$:

$$\nabla^{2}\psi=\frac{16\pi G}{3}a^{2}\rho-\frac{a^{2}}{6}\delta R\left(f_{R}\right).$$ (16)

Combining Eqs. (15) and (16) allows deriving exact solutions for extreme cases and analyzing the scale and manner in which $f(R)$ gravity diverges from GR.

In this framework, the behavior of $f(R)$ models varies depending on $f_{R0}$ values relative to the effective potential, the Newtonian potential $\Psi_{N}$. If $f_{R0}\ll\Psi_{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}\gg\Psi_{N}$, gravity would be over-amplified. Therefore, $f_{R0}$ should be close in magnitude to $\Psi_{N}$, typically between $10^{-5}$ and $10^{-6}$. A value of $10^{-4}$ 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 $\left|f_{R0}\right|\gg|\psi|$, 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 will investigate the impact of these separate components using WL tunnel voids in Sect. 5.

We now 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 $\Phi_{\mathrm{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 $\Phi_{\mathrm{len}}$ in Eq. (4), the two modified Poisson equations for the metric potentials are derived from the variation of the Einstein equations:

$$\nabla^{2}\Phi=4\pi Ga^{2}\delta\rho-\frac{c^{2}\nabla^{2}}{2}\delta f_{R},\,\,\,\,\nabla^{2}\Psi=4\pi Ga^{2}\delta\rho+\frac{c^{2}\nabla^{2}}{2}\delta f_{R}\,,$$ (17)

which, in terms of the Newtonian potential, are expressed as:

$$\Phi=\Psi_{N}-\frac{c^{2}}{2}\delta f_{R},\quad\Psi=\Psi_{N}+\frac{c^{2}}{2}\delta f_{R}\,.$$ (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 $\Phi_{\mathrm{len}}$ but only from alterations in the LSS matter distribution. Therefore, we expect a different distribution of convergence and shear compared to GR predictions. In particular, in regions of low density, the scalar field $\delta f_{R}$ contributes to an enhancement of the gravitational interaction, resulting in a more pronounced WL signal around voids than in $\Lambda$CDM models.

3 Methods

3.1 DUSTGRAIN-pathfinder simulations

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

The DUSTGRAIN-pathfinder runs are cosmological collisionless simulations that track the evolution of $768^{3}$ cold dark matter (CDM) particles within a periodic box of side length $750\,h^{-1}\,\mathrm{Mpc}$. These particles have a mass of $m^{p}_{\mathrm{cdm}}=8.1\times 10^{10}\,h^{-1}M_{\odot}$ and their gravitational softening has been set to $\varepsilon_{g}=25\,h^{-1}\mathrm{kpc}$, corresponding approximately to $1/40$-th of the mean inter-particle separation. All simulations assume a constant $\Omega_{\mathrm{m}}$ to ensure matching power spectra on large scales, with and without neutrinos. The reference $\Lambda$CDM simulation assumes GR, $M_{\nu}=0\,\mathrm{eV}$. Its cosmological parameters follow Planck2015 constraints (Planck Collaboration et al. 2016), with $\Omega_{\mathrm{m}}=\Omega_{\mathrm{cdm}}+\Omega_{\mathrm{b}}+\Omega_{\nu}=0.31345$, $\Omega_{\Lambda}=0.68655$, $H_{0}=67.31\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$, $n_{s}=0.9658$, and $\sigma_{8}=0.842$ at $z=0$.

Table 1: Summary of the main numerical and cosmological parameters of the subset of the DUSTGRAIN-pathfinder simulations considered in this work. Here, $f_{R0}$ represents the MG parameter, $M_{\nu}$ the neutrino mass, $m^{p}_{\rm cdm}$ the CDM particle mass, and $\Omega_{\rm cdm}$ and $\Omega_{\nu}$ the $\mathrm{CDM}$ and neutrino density parameters, respectively.

Simulation	Gravity	$f_{R0}$	$M_{\nu}$ [eV]	$\Omega_{\mathrm{cdm}}$	$\Omega_{\nu}$	$m^{p}_{\mathrm{cdm}}\ [h^{-1}\ M_{\odot}]$
$\Lambda$CDM	GR	-	0	0.31345	0	$8.1\times 10^{10}$
$fR4$	$f(R)$	$-1\times 10^{-4}$	0	0.31345	0	$8.1\times 10^{10}$
$fR5$	$f(R)$	$-1\times 10^{-5}$	0	0.31345	0	$8.1\times 10^{10}$
$fR6$	$f(R)$	$-1\times 10^{-6}$	0	0.31345	0	$8.1\times 10^{10}$
$\Lambda$CDM0.15	GR	-	0.15	0.30987	0.00358	$8.01\times 10^{10}$

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^{-4}$, $10^{-5}$ or $10^{-6}$. 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 $\Lambda$CDM simulation with neutrinos of mass $M_{\nu}=0.15\ \mathrm{eV}$, which we refer to as $\Lambda$CDM_0.15eV. In this case, we track the evolution of $2\times 768^{3}$ particles considering $m^{p}_{\rm cdm}=8.01\times 10^{10}\,h^{-1}M_{\odot}$ and $m^{p}_{\nu}=9.25\times 10^{8}\,h^{-1}M_{\odot}$ as dark matter and massive neutrino particle masses, respectively.

To summarize, this work employs five simulation sets: a $\Lambda$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 has been 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 use 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 process 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 stack $21$ of the $34$ available snapshots from $z_{\mathrm{min}}=0$ to $z_{\mathrm{max}}=4$. The light-cones have a square base with sides of $5\,\mathrm{deg}$, corresponding to the angular size of the simulation box at $z_{\mathrm{max}}$, and cover a total sky area of $25\,\mathrm{deg}^{2}$. 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 $\chi(z_{\mathrm{max}})$, as shown in Fig. 1. To achieve high-resolution maps of the projected matter density, the stacked snapshots are divided into $27$ lens planes. Since gravitational lensing depends on the projected matter density along the line of sight, each particle is 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.

The mass density is then interpolated from these projected particle positions to a two-dimensional 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 $\Sigma(\bm{\theta})$ on the $l$-th lens plane perpendicular to the $(0,0,1)$ direction, at pixel with coordinate indices $(i,j)$ is computed as:

$$\Sigma^{l}_{(i,j)}(\bm{\theta})=\frac{1}{A_{pix}^{l}}\sum_{k=1}^{n}m_{k}\,,$$ (19)

where $A_{pix}^{l}=4\pi/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}(\chi_{l})<f_{K}(\chi_{l^{\prime}}),\quad\forall l^{\prime}\neq l\quad\text{and}\quad\chi_{k}<\chi_{s}$. In this way, a discretization of the 3D density distribution in several mass maps is 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 is done at the same time by the MapSim tool in this first step.

To enhance the WL statistics, we generate $256$ semi-independent light-cone realizations for each cosmological model. This is achieved by randomizing the comoving simulation boxes across the redshift range $z_{\text{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, maximizing the utility of the simulation data (Jain et al. 2000).

Figure 1: Example of light-cone construction with MapSim from $z=0$ to $z=4$ for a $\Lambda$CDM simulation with a box size of $750\,h^{-1}\,\mathrm{Mpc}$, tiled to encompass the $5\times 5\,\mathrm{deg}^{2}$ 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 $z_{s}=1$. The snapshots within $z_{\text{obs}}<z<z_{s}$ are the $13$ utilized ones, while the gray region indicates the snapshots with $z>z_{s}$ 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.

3.3 Convergence maps

In order 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 are 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 constructs 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(\chi_{\mathrm{l}},\chi_{\mathrm{s}})$ (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\geq 10^{4}$) (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 use 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:

$$\Sigma^{c}_{\rm{crit},l,s}(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})}\,.$$ (20)

we can compute the convergence on the $l$-th lens plane as :

$$\kappa^{l}(\bm{\theta},\chi_{\mathrm{s}})=\frac{\Sigma^{l}(\bm{\theta})}{\Sigma^{c}_{\rm{crit},l,s}}\,.$$ (21)

At the end, the total convergence map is:

$$\kappa(\bm{\theta},\chi_{\mathrm{s}})=\sum_{l}\kappa^{l}(\bm{\theta},\chi_{\mathrm{s}})=\sum_{l}\dfrac{\Sigma^{l}(\bm{\theta})}{\Sigma^{c}_{\rm{crit},l,s}}\,.$$ (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\times 2048$ pixels, yielding a pixel resolution of approximately $9\,\mathrm{arcsec}$. The origin of each map coordinate system (R.A., Dec.) is positioned at the center such that each side spans from $-150$ to $+150$ $\mathrm{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. 2024).

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 $\kappa(\bm{\theta})$ we added GSN, $n(\bm{\theta})$, modeled as a Gaussian random field with a top-hat filter with a size that corresponds to the pixel area $A_{\mathrm{pix}}=\theta^{2}_{\mathrm{pix}}$ in $\mathrm{arcmin}^{-2}$. According to Lin & Kilbinger (2015) its variance is given by

$$\sigma_{\text{pix }}^{2}=\frac{\sigma_{\epsilon}^{2}}{2}\frac{1}{n_{\mathrm{gal}}A_{\mathrm{pix}}}\,,$$ (23)

where $\sigma_{\epsilon}^{2}=\left\langle\epsilon_{1}^{2}\right\rangle+$ $\left\langle\epsilon_{2}^{2}\right\rangle$ is the variance of the intrinsic ellipticity distribution of the source galaxies and $n_{\mathrm{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. 2024), so $\sigma_{\epsilon}^{2}=0.3$ and $n_{\mathrm{gal}}=30\ \mathrm{arcmin}^{-2}$ at $z_{s}=1$. In this way, we obtained a noised convergence map $\kappa_{n}(\bm{\theta})=\kappa(\bm{\theta})+n(\bm{\theta})$.

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 $\theta_{\mathrm{G}}$ (Kilbinger 2015). In this work, we use this approach to minimize GSN contamination in tunnel void statistics, ensuring consistency between maps with and without GSN.

Choosing an optimal value for $\theta_{\mathrm{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 $\theta_{\mathrm{G}}$ retains information from smaller scales within the WL map but leaves a significant level of GSN contamination. In contrast, a large $\theta_{\mathrm{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.

Figure 2: Example of a realization of the $\Lambda$CDM light cone, having a FOV of $5\times 5$ deg² aperture and $z_{s}=1$. In the upper panel, we show the WL convergence map, and in the lower panel, the SNR of the same map with smoothing scale at $\theta_{\mathrm{G}}=2.5\,\text{arcmin}$. The color scale in both panels represents the convergence value for each pixel in the top panel and the SNR value in the bottom panel, as indicated by the corresponding color bars.

In this work, we initially tested three smoothing scales: $\theta_{\mathrm{G}}=1$, $2.5$, and $5\,\text{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 ($\theta_{\mathrm{G}}=1$ arcmin) preserve fine structures but amplify noise, increasing variance in WL void statistics. Large scales ($\theta_{\mathrm{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 $\theta_{\mathrm{G}}=2.5$ arcmin. So, we obtain the noised and smoothed converged map as $\kappa_{2.5}(\bm{\theta})\equiv\left(\kappa_{n}*U\right)(\bm{\theta})$:

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

with $U(\bm{\theta})$ being the Gaussian window function with $\theta_{\mathrm{G}}=2.5\ \mathrm{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:

$$\sigma_{\text{noise }}^{2}=\frac{\sigma_{\epsilon}^{2}}{2}\frac{1}{2\pi n_{\mathrm{gal}}\theta_{\mathrm{G}}^{2}}\,.$$ (25)

This relation allows the derivation of the lensing SNR($\bm{\theta}$) map of the noised and smoothed convergence map $\kappa_{2.5}(\bm{\theta})$ as:

$$\text{SNR}(\bm{\theta})=\frac{\kappa_{2.5}(\bm{\theta})}{\sigma_{\text{noise}}(\theta_{\mathrm{G}})}\,,$$ (26)

where $\sigma_{\text{noise}}(\theta_{\mathrm{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 with which the WL peak or minima are identified, while it is constant only if we assume that galaxies are uniformly distributed. In our framework, with $\theta_{\mathrm{G}}=2.5$ and $z_{s}=1$, $\sigma_{\text{noise}}\approx 6.18\cdot 10^{-3}$.

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

Figure 2 illustrates the impact of applying the Euclid-like GSN and a Gaussian smoothing scale of $\theta_{\mathrm{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 SNR 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 SNR map highlights significant structures by suppressing uncorrelated noise.

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 analyze now the one-point PDF and the angular power spectrum of the constructed maps.

The PDF of the convergence field quantifies the distribution of $\kappa$ 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 $\Lambda$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 $\Lambda$CDM PDF with its noiseless counterpart.

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

$$\Delta X=\frac{X^{\rm model}(r)-X^{\rm ref}(r)}{|X^{\rm ref}(r)|}\,,$$ (27)

where $X^{\rm model}(r)$ represents the profile or observable under investigation, while $X^{\rm 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 $\Lambda$CDM model become more relevant for $\kappa<0$. Effectively, the largest deviations occur in the negative $\kappa$ range, where $f(R)$ models display an excess probability compared to $\Lambda$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 $\kappa$ 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_{\kappa}(l)$, provides a complementary statistical measure by characterizing the variance of the $\kappa$ 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\pi/\theta_{FOV}$). It is computed as the Fourier transform of the input map squared and binned in multipole space.

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 $\Lambda$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 $\Lambda$CDM are presented. The largest deviations from the $\Lambda$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\lesssim 300$), the Universe remains in the linear regime, where deviations are smaller due to the weaker impact of nonlinear modifications. At smaller scales ($l\gtrsim 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 ($\theta_{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.

Figure 4: Average angular power spectra of the noised and smoothed convergence field from $256$ light-cone random realizations with redshift $z_{s}=1$, for all the analyzed models (top panel) and the correspondent residuals with respect to the reference $\Lambda$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 $\kappa$ field.

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).

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 at https://github.com/LeonardoMaggiore/PyTwinPeaks. The structure and the main steps of the code are schematically represented in Fig. 5 and analyzed in detail in the following.

4.1 PyTwinPeaks: reconstruction of connected regions

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

The processed convergence maps serve as the input for our 2D void finder algorithm, which detects underdense regions based on their SNR. 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 SNR regions. In this work, we use PyTwinPeaks to reconstruct underdense lines-of-sight pixel regions at different convergence SNR 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 the valleys reconstruction: it begins with the creation of a mask that selects pixels in the SNR map below the chosen threshold. Next, reconstruction occurs through 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 SNR 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 SNR values of individual pixels as weights, and the coordinates of the local minimum, defined as the pixel with the lowest SNR value. The coordinates, originally in pixels, are converted to arcminutes using the relation:

$$x_{\mathrm{arcmin}}=\left(x_{\mathrm{pix}}-\frac{n_{\mathrm{pix}}}{2}\right)\cdot l_{\mathrm{arcmin}},$$ (28)

where $x_{\mathrm{pix}}$ represents the coordinate in the SNR map in pixel units, $n_{\mathrm{pix}}$ is the number of pixels per side of the map, and $l_{\mathrm{arcmin}}=\mathrm{FOV}_{\mathrm{arcmin}}/n_{\mathrm{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_{\mathrm{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:

$$D_{\mathrm{pr}}<\frac{R_{\mathrm{pr},1}+R_{\mathrm{pr},2}}{2}.$$ (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 $\text{SNR}=-5$ to $\text{SNR}=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 SNR 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—also the ones explored in this work—are through the centroid and the local minimum 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 SNR 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.

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 SNR values, the assigned center may occasionally fall outside the connected region. At higher (still negative) SNR values, connected regions can have elongated or irregular shapes, and in extreme cases—such as a crescent-shaped underdensity—the centroid may fall within or too close to a positive SNR (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 SNR 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 SNR 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), here adapted for 2D voids. The algorithm tracks local minima across multiple thresholds, starting from the lowest and checking whether each minimum persists—that is, 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 $\text{SNR}=[-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 will 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 SNR 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 SNR regions.

Figure 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 SNR weighted centroids of the connected regions. Right panel: the same 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 SNR threshold as shown in the colorbar.

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 SNR 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\times 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:

$$r_{\mathrm{mask}}=\frac{l_{\mathrm{grid}}}{2}-\frac{l_{\mathrm{pix}}}{2},$$ (30)

where $l_{\mathrm{grid}}$ represents the grid size and $l_{\mathrm{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.

Figure 7: Example of the final tunnels catalog associated with the same SNR map of Fig. 2 and Fig. 6. Colors from blue to red represent areas with increasing SNR, 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.

The algorithm expands the grid by adding a row of pixels on each side until the average SNR 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 $\langle\text{SNR}\rangle=-2.5$.
For voids located near the boundaries of the SNR 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 then 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_{\mathrm{centers}}\leq 75\%\cdot(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 also correct for the void-in-void scenario, happening 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, an SNR map from one of the $256$ convergence maps produced for the $\Lambda$CDM scenario, where circles mark the positions of the 2D voids. We can note that these voids correspond very well with the negative SNR 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 SNR 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, i.e. 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\%$.

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. 2022a), WL tunnel voids extracted in this work from SNR maps are characterized by relatively small sizes. Their angular extension in the sky-plane spans from $2$ to $15\ \mathrm{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), i.e. the total number counts of 2D voids as a function of their projected radius $R_{\rm v}$ in arcminutes, averaged from the $256$ SNR map realizations, for each cosmological model. We restrict our analysis to the range $[2,15]\ \mathrm{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 $\Lambda$CDM model, computed as in Eq. (27).

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. 2022b, a), i.e. it shows how small voids are more abundant with respect to the larger ones. The reduction in void counts at very small scales is influenced by spatial resolution. Pixel size and the smoothing applied to the SNR maps set a lower limit on the detectable void size.

Figure 8: Top panel: average void size functions, per degree squared, of 2D WL voids from $256$ light-cone random realizations with redshift $z_{s}=1$, for the different analyzed models: $\Lambda$CDM (black), $fR4$ (green), $fR5$ (blue), $fR6$ (red) and $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ (purple). Bottom panel: the corresponding residuals computed with respect to the $\Lambda$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.

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 $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ are expected to be the most similar to the $\Lambda$CDM case, although they exhibit opposite trends. Specifically, up to $R_{\rm v}\sim 10\ \mathrm{arcmin}$, we observe an amplitude difference of $\leq 10\%$ for the $fR6$ model, while for $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$, the difference is approximately $25\%$.

5.2 Stacked tangential shear profile

Figure 9: Top panel: averaged tangential shear profiles extracted from $256$ light-cone random realizations with redshift $z_{s}=1$, for the different analyzed models. Bottom panel: corresponding residual with respect to the profile measured in the $\Lambda$CDM scenario. Both the plots utilize the following color scheme: $\Lambda$CDM (black), $fR4$ (green), $fR5$ (blue), $fR6$ (red) and $\Lambda$CDM${}_{0.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 36). These are not displayed in the top panel for visualization reasons.

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 $\kappa_{2.5}(\bm{\theta})$, spanning from $r_{\rm min}=0$ to $r_{\rm max}=10R_{\mathrm{v}}$. The tangential shear $\gamma_{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 to the mean tangential shear profile of $N$ void profiles in the $i$-th bin is evaluated as:

$$\sigma_{i}=\sqrt{\frac{C_{ii}}{N}}\,,$$ (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 SNR maps. We refer the reader 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^{-3}$, matching values found in the literature for tunnel voids. In addition, we detect a $\gamma_{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, 2021b, 2021a). In particular, the observed trends are in agreement with theoretical expectations, i.e. 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.

We note, however, that there is a pivot point at $r_{p}/R_{v}\simeq 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 $\Lambda$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.

Figure 10: Top panel: averaged tangential shear profiles measured from samples of voids with different sizes in the standard $\Lambda$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 to the measures, analogously to Fig. 9.

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

For each cosmological model, we divide the entire sample of our voids according to their size, into three equi-populated bins:

• small voids with radii $R_{\rm v}$ in the range $[1,3.97]\ \mathrm{arcmin}$;

• medium voids with radii $R_{\rm v}$ in the range $]3.97,4.85]\ \mathrm{arcmin}$;

• large voids with radii $R_{\rm v}$ in the range $]4.85,15]\ \mathrm{arcmin}$.

Then we extract the void tangential shear profiles from these equi-populated sub-samples. The result is represented in Fig. 10.

Figure 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.

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 equi-populated 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.

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 follow 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. (2023), we can integrate the functional form of the density profile along the line of sight to obtain the projected surface mass density, $\Sigma(r_{p})$, that we can relate to convergence through Eq. (21) and then derive the tangential shear profile, $\gamma(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 has the form:

$$\frac{\rho_{\mathrm{v}}(r)}{\bar{\rho}}-1=\delta_{c}\frac{1-\left(r/r_{\mathrm{s}}\right)^{\alpha}}{1+\left(r/r_{\mathrm{v}}\right)^{\beta}}\,.$$ (32)

It is characterized by four free parameters ($\delta_{c}$, $\alpha$, $\beta$, $r_{\mathrm{s}}$), which represent, respectively, the central density contrast, the inner and outer slopes of the profile, and the characteristic scale where $\rho_{\mathrm{v}}=\bar{\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, $\delta_{c}$ and $s$, and is expressed as:

$$\rho_{v}\left(r/R_{\rm v}\right)=1+|\delta_{c}|\,\bar{\rho}\left\{\frac{1}{2}\left[1+\tanh\left(\frac{y-y_{0}}{s\left(R_{\rm v}\right)}\right)\right]-1\right\}\,.$$ (33)

In this parametrization, $\delta_{c}$ is the density contrast at the void center, $y=\ln\left(r/R_{\rm v}\right)$ and $y_{0}=\ln\left(r_{0}/r_{v}\right)$. The radius $r_{0}$ is determined by requiring that the integral of the profile up to $R_{\rm v}$ be $\Delta_{v}=0.2$. This allows us to express $r_{0}(s)$ (in units of $h^{-1}\ \mathrm{Mpc}$ ) as a second-order polynomial function: $r_{0}(s)=0.37s^{2}+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 $\chi^{2}$ obtained with the best-fit parameters is indeed very far from unity, i.e. $\tilde{\chi}_{\rm HSW}^{2}\simeq 2031$ and $\tilde{\chi}_{\mathrm{HT}}^{2}\simeq 7750$.

Figure 12: Best fits of the averaged stacked tunnel voids tangential shear profiles, $\gamma_{t}(r_{p})$, measured in $\Lambda$CDM cosmology adopting three different parametric functions. The data are derived as the mean of the $256$ $\Lambda$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\sigma$ 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. (2023), 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.

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

$$\gamma_{t}(r_{p})=a\cdot r_{p}\left(\frac{1-r_{p}^{b}}{1+r_{p}^{c}}-\frac{\exp(d\cdot r_{p})}{1+\exp(e\cdot r_{p}-(d+e))}\right)\ .$$ (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 $\Lambda$CDM simulations. It is easy to notice that there is an excellent match with the data across all scales. The reduced $\chi^{2}$ is indeed now closer to the unity, i.e. $\tilde{\chi}_{\rm M25}^{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 $\chi^{2}$ obtained neglecting the off-diagonal terms is $\tilde{\chi}_{\rm M25}^{2}\simeq 4.44$, indicating a strong correlation between the radial bins (see Appendix A for details).

Figure 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: $\Lambda$CDM (black), $fR4$ (green), $fR5$ (blue), $fR6$ (red), $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ (purple).

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 changes, 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 $\Lambda$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.

Figure 14: As Fig. 13 but for a selection of voids with radius in the range $]3.79,4.85]$ arcmin.

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: $\Lambda$CDM, $fR4$, $fR5$, $fR6$, $\Lambda$CDM${}_{0.15\,\mathrm{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 $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ showing the largest separation from the $\Lambda$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 $\Lambda$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 $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ – from the reference $\Lambda$CDM is over 1$\sigma$ for all the parameters of the model.
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 the reader 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_{\rm v}$ 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 $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ case appear now more detached from the reference $\Lambda$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 $\Lambda$CDM cases.
To sum up this section, we can conclude that our identified tunnel voids are sensitive to the effects of MG and massive neutrinos, and the modeling of their tangential shear profiles via the new 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 weak lensing signal of tunnel voids in the convergence field, focusing on modified gravity (MG) models within the $f(R)$ framework (Hu & Sawicki 2007) and their effects with respect to those of massive neutrinos. Using an innovative pipeline, we analyzed the signal-to-noise ratio (SNR) of noised and smoothed convergence maps from the DUSTGRAIN-pathfinder $N$-body simulations (Giocoli et al. 2018a; Peel et al. 2018), systematically comparing their behavior with the standard $\Lambda$CDM model.

In what follows, we summarize our main results.

• 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 show 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 $\Lambda$CDM${}_{0.15\,\mathrm{eV}}$ display a weaker signal due to the suppression of large-scale structure growth.

• We investigated the behavior of the tangential shear profile according to the size of voids using three equi-populated sub-samples. 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 accurately fits all stacked shear profiles, 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, especially 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 $\Lambda$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 $\Lambda$CDM${}_{0.15\,\mathrm{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. 2024) and explore a broader range of cosmological models (Euclid Collaboration: Ajani et al. 2023).
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 precise 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, i.e. 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 involve 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.

Appendix A Covariance matrix

In this Appendix, we present the tangential shear covariance matrix for the WL tunnel void sample of $\Lambda$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).

Fig. 15 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 $\theta_{\rm{s}}=2.5$ arcmin. It is computed via the formula:

$$C_{ij}=\frac{1}{N-1}\sum_{k=1}^{N}[\gamma_{\rm{t}}(i)-\bar{\gamma}_{\rm{t}}(i)][\gamma_{\rm{t}}(j)-\bar{\gamma}_{\rm{t}}(j)]\,,$$ (35)

where $i$ and $j$ are radial bin indices and $N=22203$ is the total number of samples of WL tunnel voids in the $\Lambda$CDM scenario. $\gamma_{\rm{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\times 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), $\alpha$, through the relation $C_{\text{corr}}^{-1}=\alpha C^{-1}$, where $C^{-1}$ is the inverse of the sample covariance matrix. The $\alpha$ factor is given by

$$\alpha=\frac{N-N_{\rm{bin}}-2}{N-1}\,,$$