© The Authors 2025

1. Introduction

Galaxy clusters serve as powerful tools in cosmology, acting as indicators of density peaks in the large-scale matter distribution (see e.g. Allen et al. 2011 for a comprehensive review). Their relevance lies in the potential to constrain cosmological parameters through their number density (see, for instance, Rosati et al. 2002; Vikhlinin et al. 2009; Rozo et al. 2010; Pacaud et al. 2018; Costanzi et al. 2019; Lesci et al. 2022a; Ghirardini et al. 2024; Seppi et al. 2024) and spatial distribution (as explored by Veropalumbo et al. 2014; Marulli et al. 2018; To et al. 2021; Lesci et al. 2022b; Romanello et al. 2024). Moreover, galaxy clusters can be used to study astrophysical processes influencing galaxy formation and evolution. The impact of clusters and galactic environments on the development of galaxy properties is evident in observed differences between field galaxies and those located in denser environments (e.g. Dressler 1980; Kuchner et al. 2017; George et al. 2012). The combined study of cosmology and the astrophysical properties of clusters is needed, as the two fields are deeply interconnected.

Galaxy clusters are predominantly composed of dark matter (approximately 80%), which forms haloes. About 15% of the cluster's mass is ionised hot gas within the intracluster medium (ICM), while stars and gas within the galaxies contribute to only 5% of the total mass.

Large-scale surveys have identified thousands of galaxy clusters through various methods. X-ray observations have detected the hot ICM emission (e.g. Böhringer et al. 2001; Pacaud et al. 2016; Piffaretti et al. 2011; Sadibekova et al. 2024), while the Sunyaev-Zel’dovich (SZE) effect has been observed at millimetre wavelengths (Schaffer et al. 2011; Hasselfield et al. 2013; Bleem et al. 2020; Planck Collaboration IX 2011; Aiola et al. 2020; Koulouridis et al. 2021). Optical and near-infrared surveys have also been successful in identifying clusters, either by detecting overdensities of galaxies (Kepner et al. 1999; Ramella et al. 2001; Sarron et al. 2018; Werner et al. 2023) or by exploiting known properties of the clusters, such as the luminosity or colour distribution of their galaxies (Koester et al. 2007; Gilbank et al. 2011; Rykoff et al. 2014; Farrens 2011; Bellagamba et al. 2018; Maturi et al. 2019; Eisenhardt et al. 2008; Andreon et al. 2009; Rettura et al. 2014; Euclid Collaboration: Adam et al. 2019).

However, these methods primarily rely on baryonic observables (X-ray luminosity, SZE flux, optical richness), which only account for a small fraction of the total cluster mass. This reliance on baryonic matter can introduce selection biases, such as favouring cool-core clusters with centrally concentrated X-ray emission (Planck Collaboration I 2011; Andreon & Moretti 2011; Eckert et al. 2011; Andreon et al. 2024). Brighter clusters are more easily detected, while fainter ones may be missed. By using complementary methods that incorporate independent observables, such as weak lensing, these biases can be mitigated, enabling the detection of both brighter and fainter clusters with greater reliability.

Clusters induce gravitational lensing effects on the light emitted by background galaxies (see e.g. Schneider 1992; Pyne 1995). In the weak lensing regime, the distortion of background galaxy shapes is detectable only statistically (see e.g. Schneider et al. 2002; Pires et al. 2020). Weak-lensing cluster detection stands out by using matter concentration as a tracer, independent of the physical state of the baryonic component. This approach allows the detection of clusters regardless of their emission (X-rays, optical) or scattering (SZ) properties. However, the sensitivity of cluster detection through weak lensing critically relies on the presence of background sources available for lensing.

In the early stages of optical surveys, one of the first techniques for the detection of galaxy clusters through weak lensing signals was the aperture mass (AM) method, introduced by Schneider et al. (1998). This approach convolved the lensing signal with a filter function that depends on a specific scale, providing a foundation for subsequent advancements in cluster detection techniques. Around the same time, a straightforward Gaussian filter applied to the reconstructed convergence also gained popularity due to its simplicity and effectiveness in detecting clusters (e.g. White et al. 2002; Miyazaki et al. 2002, 2007; Hamana et al. 2004; Tang & Fan 2005; Gavazzi & Soucail 2007; Fan et al. 2010; Shan et al. 2012, 2018). These early methods led to the development of more refined approaches, including optimisations of the AM method through the design of various filters. These filters have since been rigorously tested on simulations and applied to numerous optical surveys (e.g. Schneider 1996; Schneider et al. 1998; Jarvis et al. 2004; Schirmer et al. 2004, 2007; Hetterscheidt et al. 2005; Hennawi & Spergel 2005; Maturi et al. 2005, 2007; Pace et al. 2007; Dietrich & Hartlap 2010; Hamana et al. 2012, 2020; Lin & Kilbinger 2015; Miyazaki et al. 2018; Oguri et al. 2021; Chen et al. 2025). More advanced filtering techniques included the full information on the expected shape of the halo profile in the kernel of the filter. This improvement aimed to reduce the contamination by shape noise and the lensing contribution of the large-scale structure (LSS) to the error (e.g. Hennawi & Spergel 2005; Maturi et al. 2005; Wittman et al. 2006). Furthermore, given the diverse properties and scales of the clusters to be detected, multi-scale filtering techniques have also been explored, yielding promising results (e.g. Starck et al. 2006; Leonard et al. 2015; Lanusse et al. 2016; Leroy et al. 2023).

In this study, we present and test the performance of AMICO-WL, an augmentation of the Adaptive Matched Identifier of Clustered Objects (AMICO) (Bellagamba et al. 2018; Maturi et al. 2019) that includes a weak-gravitational lensing branch for the detection of galaxy clusters. AMICO is a representative of the class of linear optimal matched filters; this algorithm is effective in extracting signals with a maximised signal-to-noise ratio (S/N). AMICO was officially adopted as one of two algorithms for optical cluster detection by the ESA Euclid mission¹ (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2025). Its performance stood out in the Euclid Cluster Finder Challenge (Euclid Collaboration: Adam et al. 2019), demonstrating good completeness and purity as a function of redshift when applied to mock galaxy photometric catalogues replicating the expected properties of Euclid photometric catalogues. The AMICO algorithm has already demonstrated success in wide-field and/or deep surveys such as the Kilo-Degree Survey (KiDS²; de Jong et al. 2016), J-PAS (Maturi et al. 2023) and COSMOS (Toni et al. 2024). The analysis of the KiDS data resulted in a cluster sample (Maturi et al. 2019) that has been used in various cosmological studies (e.g. Giocoli et al. 2021; Ingoglia et al. 2022; Lesci et al. 2022a, b) and investigations of galaxy properties in cluster environments (e.g. Radovich et al. 2020; Puddu et al. 2021; Castignani et al. 2022, 2023).

AMICO-WL implements an adjusted version of the optimal linear matched filter defined in Maturi et al. (2005). The filter assumes a mean radial profile of the halo shear pattern and a power spectrum for the noise. The efficient AMICO infrastructure is adapted for weak lensing data, for example, by changing the cleaning procedure. Moreover, the pipeline is complemented with new tools developed to increase the reliability of the detection, such as a thresholding method based on the S/N analysis of the map peaks and an effective foreground removal procedure for the reduction of noise from galaxies between the lens and the observer.

In this paper, we use as a case study the application of AMICO-WL to mock data from dedicated light-cone simulations, replicating the expected source density of the Euclid wide-field survey (Euclid Collaboration: Scaramella et al. 2022). The Euclid mission, with its wide-field survey covering approximately 14 000 deg², will measure the shapes of billions of background galaxies, making the weak-lensing probe sensitive to the detection of clusters. For this reason, the Euclid Collaboration initiated the Weak Lensing Selected Cluster Challenge, a platform to compare different analysis methods. AMICO-WL stands as a key contender among the participating algorithms. In the second phase of the challenge, using twelve 100 deg² fields of simulated Euclid data, AMICO-WL demonstrated good performance, achieving significant results in purity and completeness. The comprehensive results of this challenge will be presented in Manjón-García et al. (in prep.).

In this work, we quantify the performances of AMICO-WL by implementing a new matching procedure based on the ‘blinking’ of redshift planes in the simulated data to investigate the nature of the detections and the functioning of the detection method in more detail.

This paper is organised as follows. In Sect. 2, we summarise the fundamental points of gravitational lensing theory. The simulated dataset and reference catalogue are presented in Sect. 3. We describe the optimal filtering technique on which AMICO-WL is based and the foreground removal cleaning applied to the dataset in Sect. 4. Section 5 shows the steps of the AMICO-WL algorithm, from the creation of the map to the computation of the S/N threshold and the cleaning procedure. The results of the application of AMICO-WL to the mock data are evaluated in Sect. 6. In Sect. 7, we discuss the nature of the spurious detections. Finally, in Sect. 8, we outline our conclusions.

2. Weak gravitational lensing

2.1. Theory

In this section, we briefly summarise the aspects of weak gravitational lensing relevant to the present study. For a more detailed description of the lensing theory, we refer to Bartelmann & Schneider (2001).

Let us consider a lens at the angular position vector θ on the sky. We can write the surface mass density as

$$\mathrm{\Sigma }(\theta )={\displaystyle {{\displaystyle \int }}_{-\infty }^{+\infty }}\rho (\theta ,z)\mathit{dz},$$(1)

where ρ is the mass density, from which we can define the lensing convergence as

$$\kappa (\theta )\equiv \frac{\mathrm{\Sigma }(\theta )}{{\mathrm{\Sigma }}_{\mathrm{crit}}},$$(2)

with the critical surface density Σ_crit defined as

$${\mathrm{\Sigma }}_{\mathrm{crit}}\equiv \frac{c^2}{4\pi G}\frac{D_{\mathit{ls}}}{D_l{D}_s}·$$(3)

The angular diameter distances D_s, D_l, and D_ls refer to the distances between the observer and the source, the observer and the lens, and between the lens and the source, respectively.

Using the two-dimensional Poisson equation, the relation between the lensing potential ψ, a scalar field encoding the light deflection, and the convergence is given by

$${\nabla }_{\theta }^2\psi (\theta )=2\kappa .$$(4)

The deflection angle,

$$\alpha (\theta )\equiv {\nabla }_{\theta }\psi (\theta ),$$(5)

expresses the difference between the true angular position in the source plane β, and the observed position θ, through the lens equation

$$\beta =\theta -\alpha (\theta ).$$(6)

Because of this angular distortion, the shapes of the lensed images differ from the intrinsic shape of the source because light bundles are deflected differentially. If a source is much smaller than the angular scale on which the lensing properties change, the lens mapping can be linearised. The lens equation in the first-order lens mapping can be written as

$$\beta =\mathcal{A}(\theta )\theta ,$$(7)

where the Jacobian matrix $\mathcal{A}$ describes the distortion of the images

$$\mathcal{A}(\theta )=\frac{\partial \beta }{\partial \theta }=\left({\delta }_{\mathit{ij}}-\frac{{\partial }^2\psi (\theta )}{\partial {\theta }_i\partial {\theta }_j}\right)=\left(\begin{array}{cc}\hfill 1-\kappa -{\gamma }_1\hfill & \hfill -{\gamma }_2\hfill \\ \hfill -{\gamma }_2\hfill & \hfill 1-\kappa +{\gamma }_1\hfill \end{array}\right).$$(8)

The shear components γ₁ and γ₂ are directly related to the second derivatives of the lensing potential Ψ through

$${\gamma }_1\equiv \frac{1}{2}({\partial }_1^2-{\partial }_2^2)\mathrm{\Psi }$$(9)

and

$${\gamma }_2\equiv {\partial }_1{\partial }_2\mathrm{\Psi }.$$(10)

The convergence κ can also be expressed in terms of the lensing potential as

$$\kappa \equiv \frac{1}{2}\left({\partial }_1^2+{\partial }_2^2\right)\mathrm{\Psi }.$$(11)

From the previous equations, where we used the Einstein derivative notation, it is clear that κ describes the isotropic contraction or dilatation of an image while γ is the anisotropic distortion. Moreover, it is evident that the shear γ=γ₁+iγ₂ can be directly derived from the potential Ψ. The reverse does not always hold, and a general shear field γ might not yield a physical lensing potential. This leads to the decomposition of the shear field into two independent modes, E and B, associated with scalar potentials Ψ_E and Ψ_B. The shear components are then written as

$${\gamma }_1=\frac{1}{2}({\mathrm{\Psi }}_{E,11}-{\mathrm{\Psi }}_{E,22})-{\mathrm{\Psi }}_{B,12},$$(12)

and

$${\gamma }_2=\frac{1}{2}({\mathrm{\Psi }}_{B,11}-{\mathrm{\Psi }}_{B,22})+{\mathrm{\Psi }}_{E,12}.$$(13)

Comparing the two expressions above with Eqs. (9) and (10) we note that for gravitational lensing, the E-mode potential is equivalent to the lensing potential (Ψ_E(θ) = Ψ(θ)) and the B-mode potential is zero (Ψ_B(θ) = 0), i.e. for gravitational lensing the E-mode field is related to the shear γ, while the B-mode field is equal to zero.

The spin-2 shear components, γ₁ and γ₂ are defined with respect to a reference Cartesian coordinate frame. Therefore, it is useful to consider these components with respect to a given reference point, such as the cluster centre, to define the tangential γ₊(θ′;θ) and cross γ_×(θ′;θ) components of the shear at a given position θ′ with respect to θ. These quantities are directly observable in the weak-lensing regime, where |κ|≪1, |γ|≪1. Averaging the measured values within a circle of radius ϑ around θ, we can write

$$\langle {\gamma }_{+}\rangle (\vartheta ;\theta )=\frac{\mathrm{\Delta }\mathrm{\Sigma }}{{\mathrm{\Sigma }}_{\mathrm{crit}}}(\vartheta ;\theta ),$$(14)

and

$$\langle {\gamma }_{\times }\rangle (\vartheta ;\theta )=0,$$(15)

where ΔΣ(ϑ;θ) represents the excess surface mass density (Umetsu 2020). Equation (14) shows that the average tangential component of the shear is directly related to the excess surface density. Equation (15) highlights that for a pure gravitational lensing signal (where B-modes are absent), the average cross-component of the shear is expected to be statistically consistent with zero. The presence of non-zero 〈γ_×〉 can indicate the presence of systematic uncertainties or non-lensing effects that introduce a B-mode component to the shear field. Thus, the measurement of the B-mode signal provides a useful null test against the presence of systematic uncertainties due to correlated noise.

2.2. Weak lensing of galaxies

In the weak lensing regime, it is common to use the complex ellipticity ϵ=ϵ₁+iϵ₂=(a−b)/(a+b)e^i2ϕ, where a and b are the semi-major and semi-minor axes of the galaxy image, and ϕ is the orientation of the major axis with respect to a reference frame. A source with intrinsic ellipticity ϵ_s displays an image with ellipticity

$$ϵ=\frac{{ϵ}_s+g}{1+g^{*}{ϵ}_s},$$(16)

where g represents the reduced shear, defined as

$$g(\theta )\equiv \frac{\gamma (\theta )}{1-\kappa (\theta )}·$$(17)

Equation (16) shows that the reduced shear determines the observed elongation, combining the effect of lensing and the intrinsic shape of the sources. Since background galaxies are expected to be randomly oriented, the average value of intrinsic source ellipticity 〈ϵ_s〉 is expected to be close to zero. Thus, the average value of the ellipticity of the images represents an estimate of the reduced shear: 〈ϵ〉≃g. Furthermore, since we are in the weak-lensing regime and κ≪1, thus g≃γ, we can conclude that the ellipticity of a galaxy image is an unbiased estimate of the local shear. We can then define the tangential ϵ₊ and cross ϵ_× components of the ellipticity in the same reference frame as γ₊ and cross γ_×.

3. Weak lensing simulations

To evaluate the performance of our methodology in identifying galaxy clusters via weak gravitational lensing (see the next section), we used dedicated light-cone simulations.

We built mass density planes piling the boxes of a cosmological simulation in a light-cone using the MapSim algorithm (Giocoli et al. 2015) and shot the rays to construct lensing maps. This method was already tested and used on a variety of cosmological simulations (Tessore et al. 2015; Castro et al. 2018; Giocoli et al. 2018a) and recently compared and validated with other algorithms (Hilbert et al. 2020), finding only per cent level differences for both cosmic shear two-point and peak statistics. We employed the MapSim pipeline to read the snapshots from an N-body simulation realisation and build a light-cone up to a redshift z = 4. For the N-body run we made use of the reference ΛCDM simulation of the DUSTGRAIN-pathfinder suite (Giocoli et al. 2018b) which considers a volume of (750 h⁻¹ Mpc)³ filled with N_p = 768³ particles, corresponding to a particle mass resolution m_p = 8×10¹⁰ h⁻¹ M_⊙. The initial conditions have been generated considering cosmological parameters consistent with the Planck-2015 data analysis (Planck Collaboration I 2016), namely total matter density Ω_m = 0.31345, baryon density Ω_b = 0.0491, Hubble constant H₀ = 67.31 km s⁻¹ Mpc⁻¹, scalar spectral index n_s = 0.9658, and mean amplitude of linear density fluctuations on the 8 h⁻¹ Mpc scale σ₈ = 0.842. In particular, we used 21 snapshots in the redshift range between z = 0 and z = 4. Given the box length of 750 h⁻¹ Mpc, 7 boxes are needed to cover the comoving distance of about 5 h⁻¹ Gpc to a source redshift z_s = 4. The volume required to construct the light-cone is divided along the line of sight into multiple contiguous redshift slices, combining the individual snapshots to obtain better redshift sampling. If the redshift slice extends beyond the boundary of a single box, two lens planes are constructed from a single snapshot. The total number of lens planes up to z_s = 4 is 27. To avoid replicating the same structure along the line of sight, the periodic boxes are randomised by redefining the centre of the box and changing the ordering and signs of the axes (Roncarelli et al. 2007). From the line of sight plane realisations, we then shot the rays using the Born approximation and constructed convergence and shear maps on a square patch of 5×5 deg² for different source redshifts, considering a pixel resolution 2024×2024 on a grid.

In Fig. 1, we show the convergence map obtained assuming a source redshift z_s = 1. The left panel displays the convergence map obtained considering all matter along the line of sight, while in the central and right panels, we show the same map where we discard the matter distribution present between z = 0.2–0.3 and z = 0.3–0.4, respectively. The coloured squares highlight some particular features that are not present in the central and right panels. In particular, the red squares are related to the ones missing in the z = 0.2–0.3 plane, while the yellow squares are the ones missing between z = 0.3–0.4. The most noticeable case is the central cluster that is located between z = 0.2–0.3 and is not present in the middle panel. This shows that, by removing individual planes one at a time, we can, therefore, use a ‘blinking’ approach to gain some information about the redshift: this will be useful when matching dark matter haloes with weak lensing detections.

Fig. 1. Convergence map (5 deg on a side) for fixed source redshift zs = 1. In the left panel, we consider all the matter present along the line of sight from z = 0 to z = 1, while in the central and right panels, we remove the matter located in a redshift slice as indicated in the legend. The red and yellow squares highlight some particular high-density regions that disappear in the central and right panels when discarding the matter present in the considered redshift slices.

We then populated the light-cone using a Euclid-like source redshift distribution (see Fig. 2); we assumed a distribution normalised to the total value of 30 galaxies per square arcmin, with a peak around z = 1 (Boldrin et al. 2012, 2016; Euclid Collaboration: Giocoli et al. 2024), and we randomly distribute them in the map, neglecting source clustering. This gave us the possibility to construct a shear catalogue that will be the input for our AMICO-WL code. In computing the lensing quantities, we interpolated between two consecutive shear planes. To display the impact of the considered source redshift distribution on the lensing quantities, in Fig. 3, we show the resulting convergence map. In the left panel, we consider all matter between z = 0 and z = 4, while in the central and right panels, we remove matter between z = 0.2–0.3 and z = 0.4–0.5, respectively. The effect of the missing structures in the central and right panels is, in this case, modulated by the source redshift distribution.

Fig. 2. Redshift distribution of the weak lensing simulated sources. The black empty histogram refers to the COMPLETE catalogue, while the hatched histograms represent the foreground-removed catalogues FR02, FR04 and FR06 (described in Sect. 4.2), shown as green, light-blue and magenta, respectively.

Fig. 3. Reconstructed convergence map from the source redshift distribution catalogue. The left panel refers to the full light-cone up to z = 4, while in the central and right panels, we show the maps where we exclude the matter located between z = 0.2−0.3 and z = 0.4−0.5, respectively.

We also included the expected shear noise in the simulation to mimic realistic surveys. Uncertainties in the shear – referred to as shape noise – arise from a combination of the intrinsic shape of the galaxies and measurement errors, including, among other factors, uncertainties in the galaxy shape estimation and the point spread function correction. With the intrinsic shape of the galaxies being the dominant component and the galaxies being randomly distributed, the shape noise can be modelled as additive noise. This can be well approximated by a Gaussian distribution with a mean of μ = 0 and a standard deviation of σ_ϵ = 0.26 (e.g. Leauthaud et al. 2007; Schrabback et al. 2015, 2018). For each galaxy, the shape noise was, therefore, included in the catalogue by adding Gaussian noise to the two components of the shear. The final catalogue of galaxies is composed of the angular positions (i.e. right ascension and declination), the two perturbed components of the shear, and the true redshift – which means that we do not consider any photometric redshift uncertainty in our analysis.

Our MapSim pipeline complements the particle light-cone with the line of sight corresponding halo and sub-halo catalogues that were created at each simulation snapshot by applying a friend-of-friends algorithm (Davis et al. 1985) to identify the collapsed dark matter haloes. The algorithm used a linking distance λ = 0.16×d, where d represents the mean separation distance between particles. Subsequently, the SUBFIND algorithm was applied to determine for each friend-of-friends-identified halo the mass M₂₀₀ and the radius R₂₀₀, which correspond to the mass and the radius of a spherical region around the fiducial centre of each halo, enclosing 200 times the critical density of the Universe, and the corresponding sub-halo population. The resulting catalogue comprises the angular positions (i.e. right ascension and declination) of the identified dark matter haloes and sub-haloes, along with their estimated redshift (z_l), mass, and radius. This is the catalogue used in the matching procedure to test the performance of the AMICO-WL algorithm (see Sect. 6.1).

4. Optimal filtering technique

In this section, we describe the technique on which AMICO-WL is based, the optimal matched filtering, and a linear weak-lensing estimator.

4.1. Definition of the optimal filter

The optimal matched filter implemented in AMICO-WL was first presented in Maturi et al. (2005) and Maturi et al. (2007) for an unbiased detection of the weak lensing signal of dark matter haloes aiming at a minimum variance in the estimates, achieved by suppressing the contamination of the large-scale structures and thus reducing the number of spurious detections. The performance of the weak lensing application of the optimal filter was then compared with other weak lensing estimators in Pace et al. (2007). In contrast to the standard implementation of AM, the shape of the optimal filter is determined not only by the shear profile of the lens but also by the correlation properties of the noise affecting the weak lensing measurements. Here, we model the contribution of virialised objects by assuming a radial profile for the shear signal and the properties of the LSS by the corresponding power spectrum. In particular, for galaxy clusters, we assume that the shear signal is produced by a Navarro-Frenk-White (NFW) profile (Navarro et al. 2004, 1996).

The weak gravitational lensing signal of a dark matter halo, S(θ), is given by an amplitude, A, and an angular shape, τ(θ). The data we measure, D(θ), are contaminated by the noise, N(θ), so that they can be modelled as

$$D(\theta )=S(\theta )+N(\theta )=A\tau (\theta )+N(\theta ).$$(18)

The linear filter is defined as the convolution kernel, Ψ(θ), applied to the data, D(θ), in order to yield an optimal estimate, A_est, of the amplitude of the signal at the position θ

$$A_{\mathrm{est}}(\theta )={\displaystyle {\int }_{{ℝ}^2}}D({\theta }^{\prime })\mathrm{\Psi }(\theta -{\theta }^{\prime }){\text{d}}^2{\theta }^{\prime }.$$(19)

The filter used in this linear estimator must satisfy two constraints. First, A_est must be unbiased, i.e. its average b over many realisations has to vanish

$$b\equiv \langle A_{\mathrm{est}}-A\rangle =A\left[{\displaystyle {\int }_{{ℝ}^2}}\tau ({\theta }^{\prime })\mathrm{\Psi }({\theta }^{\prime }){\text{d}}^2{\theta }^{\prime }\right]=0.$$(20)

Second, the measurement noise σ², determined by the mean-squared deviation of the estimate from its true value,

$$\begin{array}{cc}\hfill {\sigma }^2=\langle (A_{\mathrm{est}}-A{)}^2\rangle =b^2+\frac{1}{(2\pi {)}^2}{\int }_{{ℝ}^2}|\hat{\mathrm{\Psi }}(k){|}^2{P}_N(k){\text{d}}^{2}k& \hfill \end{array}$$(21)

has to be minimal. To find the function Ψ that satisfies these two conditions, we combine them by means of a Lagrangian multiplier λ, we carry out the variation of L=σ²+λb with respect to Ψ and thus find the filter Ψ minimising L. The solution to this constrained minimisation in the Fourier domain is given by

$$\hat{\mathrm{\Psi }}(k)=\frac{1}{(2\pi {)}^2}{\left[{\displaystyle {\int }_{{ℝ}^2}}\frac{|\hat{\tau }(k){|}^2}{P_N(k)}{\text{d}}^{2}k\right]}^{-1}\frac{\hat{\tau }(k)}{P_N(k)},$$(22)

where $\hat{\tau }(k)$ is the Fourier transform of the expected shear profile of an NFW halo and P_N is the noise power spectrum. The filter is therefore constructed to be most sensitive for those spatial frequencies where the signal $\hat{\tau }$ is large and the noise P_N is small. Finally, we back-Fourier transform $\hat{\mathrm{\Psi }}(k)$ into the real space to perform the estimates.

For weak lensing observations, the signal can be either the convergence or the shear. However, since the convergence is not directly measurable, it is convenient to use the shear and, specifically, the tangential shear component. In this case Eq. (19) takes the form

$$A_{\mathrm{est}}(\theta )={\displaystyle {\int }_{{ℝ}^2}}{\gamma }_{\mathrm{obs}}({\theta }^{\prime };\theta )\mathrm{\Psi }(\theta -{\theta }^{\prime }){\text{d}}^2{\theta }^{\prime }.$$(23)

The noise is assumed to be given by three contributions, namely the contributions from the finite number of background sources, the intrinsic random ellipticities and orientations, and the correlated weak-lensing signal due to the large-scale structure of the Universe. The first two are characterised by the power spectrum

$$P_{ϵ}=\frac{1}{2}\frac{{\sigma }_{{ϵ}_s}^2}{n_g},$$(24)

which depends on the dispersion of the intrinsic ellipticities of the sources σ_ϵs, and on the number density of background galaxies n_g. As seen in Sect. 3, the σ_ϵs also takes into account the measurement noise. The contribution from the large-scale structures can be described by the shear power spectrum P_γ(l), thus the total noise power spectrum can be written as

$$P_N(k)=P_{ϵ}(k)+P_{\gamma }(k),$$(25)

where P_γ is computed assuming non-linear theory models (Lewis et al. 2000; Smith et al. 2003; Takahashi et al. 2012).

4.2. Filter implementation and foreground removal

In our analysis, we implemented four different realisations of the optimal filter: together with the whole galaxy catalogue, we considered three redshift cuts, corresponding to all galaxies with redshift z>z_min, with z_min = 0.2, z_min = 0.4 and z_min = 0.6. These redshift distributions are also shown in Fig. 2. This foreground removal strategy is adopted in order to clean the dataset from the noise that foreground galaxy ellipticities produce on the weak lensing measures on the lens planes located at a higher redshift. The redshift-cut shear catalogues allow us to select galaxies that are primarily background sources, thus serving as ideal targets for weak lensing detection of dark matter haloes. Galaxies below the redshift thresholds z_min could be background sources for some clusters but often act as foreground objects for the majority of clusters of interest. This foreground contamination can introduce noise into the weak lensing signal, hindering our ability to detect and accurately characterise dark matter haloes.

The optimal filter must be optimised according to the specific properties of each truncated catalogue. For this purpose, the cosmic shear power spectrum and the effective density of galaxies have to be adopted. In Table 1, we summarise the properties of the four filters. For the template halo, we assume a standard set of parameters already used, for instance, in Maturi et al. (2007). For the COMPLETE catalogue, the parameters of the template halo are meant to detect high-mass haloes that have an intermediate distance between the observer and the sources, a condition for which the lensing effect is maximised due to the geometrical dependencies of the lensing strength (Pace et al. 2007). For the filters applied to the foreground removed catalogues, the template lens redshift is set instead to z_min, the redshift value of the threshold used for the cut. The computation of the filters needs the contribution of the power spectrum given by the ellipticity of the galaxies, Eq. (24). Last row of Table 1 represents the instrumental beam FWHM, i.e. the resolution of the filter, computed as the inverse of the surface number density of galaxies. In Fig. 4, we show the profile of the four filters in real space. From the figure, it is clear that the lower the redshift of the template model, the shallower the central peak of the filter, i.e. the broader the filter. This is clear by noting the similarity in the profile of the COMPLETE and FR04 optimal filters, where the lens redshift is set tos z_l = 0.4.

Fig. 4. Radial profiles (in real space) of the optimal filters employed for the four AMICO-WL runs: COMPLETE (black solid line), FR02 (green dashed line), FR04 (cyan dash-dotted line) and FR06 (pink dotted line).

From here on, we will refer to the application of the first optimal filter shown in Table 1 to the catalogue of galaxy ellipticity without cuts as COMPLETE, while the application of the other three optimal filters with lens redshift cut equal to z_min = 0.2, z_min = 0.4 and z_min = 0.6 to the correspondent foreground removed sub-catalogue as FR02, FR04 and FR06, respectively.

5. AMICO-WL algorithm

The AMICO algorithm (Bellagamba et al. 2011, 2018) was developed to identify galaxy clusters using photometric galaxy catalogues. After a dedicated challenge presented in Euclid Collaboration: Adam et al. (2019), it has been selected by the Euclid Collaboration as one of the algorithms that will be used to produce the official Euclid galaxy cluster catalogues. The original version of the code, assuming a model for the density profile and for the luminosity function of galaxies in clusters, identifies galaxy groups and clusters by finding overdensities and matching them with the defined model based on the optimal filter technique, as presented in Sect. 4. The versatility of Optimal Matched Filtering offers the possibility to implement other observables, as, in our case, the ones based on weak-gravitational lensing.

In this work, we modified AMICO to process the gravitational lensing signal of galaxy clusters derived from the ellipticity measurements of background galaxies. We implemented a new weak-lensing branch in AMICO by exploiting the polymorphism of C++ to minimise the modifications of the well-tested existing code. We call this modified version AMICO-WL. Our algorithm follows the layout of the original code, which can be broken down into two main steps:

1. Map creation: the optimal filter technique is applied to the data, producing amplitude and variance maps of the signal.

2. Detection: the S/N map is computed, and the clusters are detected as peaks, also employing a ‘cleaning’ procedure.

In this section, we describe AMICO-WL in detail, first explaining in Sect. 5.1 the method used to create maps and then describing in Sect. 5.2 the detection algorithm, including the implemented cleaning procedure and in Sect. 5.3 the method for the computation of the optimal S/N detection threshold.

5.1. Signal map creation

In order to compute the signal maps, the algorithm needs a catalogue of galaxies reporting their positions (Cartesian or sky coordinates), the two components of the ellipticity, the redshift z_g, necessary for the foreground removal, and their corresponding weak lensing weight w (if available).

The algorithm identifies all galaxies that fall in each pixel, then computes the lensing contribution related to those galaxies in the surrounding region according to the filter value and finally adds to the amplitude map, taking into account the weight provided by the filter Ψ. This is done by replacing the integral in Eq. (19) with a sum over all galaxies, computed using the tangential component of the ellipticity ϵ₊, which is a good estimator of γ₊ (see Sect. 2.2), as

$$E_{\mathrm{map}}(\theta )\equiv A_{\mathrm{est}}(\theta )=\frac{1}{n_g}{\displaystyle \sum _k}{ϵ}_{+}(\theta ;{\theta }_k)w_k\mathrm{\Psi }(|{\theta }_k-\theta |),$$(26)

where the sum is extended to all galaxies within a circular area centred on θ and extending up to a cut radius, which is a user-defined parameter. If the lensing weights w are not used, they are internally set to unity. The normalisation factor n_g is the effective number density of galaxies computed as the sum of the lensing normalised statistical weights of the galaxies. Along the amplitude map, the code produces the variance map, which can be written as

$${\sigma }_{E_{\mathrm{mode}}}^2(\theta )\equiv {\sigma }_{A_{\mathrm{est}}}^2(\theta )=\frac{1}{2n_g^2}{\displaystyle \sum _k}|{ϵ}_{+}(\theta ;{\theta }_k){|}^2{w}_k^2{\mathrm{\Psi }}^2(|{\theta }_k-\theta |).$$(27)

Note that this noise estimate only accounts for the galaxy-shape noise. To derive the total noise, we add in quadrature the LSS noise contribution given by

$${\sigma }_{\mathrm{LSS}}^2=\frac{1}{2\pi }{\displaystyle {\int }_0^{\infty }}|\hat{\mathrm{\Psi }}(k){|}^2{P}_{\gamma }(k)k\text{d}k.$$(28)

The B-mode is evaluated by exploiting the same algorithmic infrastructure by simply replacing the ellipticity ϵ₊ in Eq. (26) with the cross component of the ellipticity ϵ_×. This B-mode can be used to estimate the noise fluctuations and possible systematics in the shape measurement pipeline, as they do not contain any lensing contribution by construction (see Eq. (15)). The variance map for the B-mode is then obtained with the sum of the ${\sigma }_{B_{\mathrm{mode}}}^2$ and the ${\sigma }_{\mathrm{LSS}}^2$. This is done to ensure consistency in the statistical interpretation of the filtered signal, the filter is constructed assuming a power spectrum of the noise also for the LSS. This allows the B-mode to be used not only as a systematics check on the signal map, but also as a statistical null test for the expected variance structure under the same filtering assumptions.

In Fig. 5, we show the E-mode amplitude (top-left panel) and the corresponding variance map (bottom-left panel) returned by AMICO-WL in the total 5° × 5° field of view for the COMPLETE case, which considers all galaxies with no redshift cuts. The haloes with the largest signal are visible as clear peaks in the top panel. Comparing the amplitude map obtained by AMICO-WL with the corresponding simulated convergence map (Fig. 1) we notice that the largest peaks in the E-mode map clearly match the structures of the convergence map. An accurate definition of the detected peaks and the matching with the structures in the simulations will be discussed in Sect. 6.1. In Fig. 5 we also show the B-mode amplitude (top right panel) and variance maps (bottom right panel). The lensing signal peaks present in the E-mode cannot be found in the B-mode signal map, which is only produced by noise and thus appears as the background map of the E-mode. On the other hand, from the comparison of the variance maps, we see that they are very similar, as expected from their consistent construction.

Fig. 5. The E-mode amplitude and its variance lensing signal maps are shown in the top-left and bottom-left panels, respectively, while the B-mode amplitude and its variance lensing signal maps are shown in the top-right and bottom-right panels. The maps cover a field of view 5° × 5° and refer to the COMPLETE run.

5.2. S/N threshold selection

Once the amplitude and variance maps have been created, the code procedes with the detection of the peaks. The detection is not carried out on the amplitude map but on the S/N map, which AMICO-WL computes on each pixel as the ratio between the amplitude value and the square root of the median value of the variance map. AMICO-WL also provides the possibility to compute S/N locally, i.e. by computing the noise using the variance map value in the pixel instead of the median value. Since the variance is highly homogeneous, as shown in the bottom left panel of Fig. 5, using the non-local value of S/N improves the stability of the detection procedure. However, this approach is reliable whenever the uncertainties are homogeneous on large scales.

To determine the minimum S/N used to define a reliable sample of detections, we developed a method that exploits the statistical properties of the positive and negative E-mode peaks in these maps. As a technical note, we use SExtractor (Bertin & Arnouts 1996) to detect the peaks in the maps down to very low S/N, and when looking for the negative peaks we simply invert the sign of the map. The positive peaks in the E-mode map are generated by (i) lensing produced by galaxy clusters, (ii) LSS overdensities or (iii) noise. In contrast, the negative minima of the E-mode map are generated by (i) noise and (ii) lensing signal from LSS underdensities³ but not by galaxy clusters. We exploit this key difference between the positive and negative peaks to evaluate the sample purity by comparing the S/N cumulative distribution obtained from the original E-mode amplitude map (E_S/N) and that of the absolute values of the negative peaks (${\hat{E}}_{S/N}$). The S/N cumulative distribution of the peaks for the COMPLETE run is shown in the upper panel of Fig. 6.

Fig. 6. S/N threshold determination analysis. Top panel: Cumulative S/N distributions of positive (ES/N) and absolute negative peaks (${\hat{E}}_{S/N}$) in the E-mode map detected with SExtractor. Bottom panel: Difference in cumulative counts of positive and negative peaks as a function of the S/N. The S/N threshold achieving a 70% purity is defined where the cumulative positive peak count reaches 70% of the total (horizontal line) and is marked across both panels (red vertical dashed line).

As the statistics of the negative peaks reflect the presence of spurious detections due to noise and the LSS contribution, as already discussed, we use it to define the sample purity as

$$\mathcal{P}=\frac{E_{S/N}-{\hat{E}}_{S/N}}{E_{S/N}}·$$(29)

The purity values as a function of S/N are shown in the bottom panel of Fig. 6: according to the desired value of purity from this distribution it is possible to select a S/N threshold that AMICO-WL uses as the lower detection threshold. For example, a sample purity of 70% is expected to be obtained with a threshold of S/N = 3.08 (shown by the vertical red line in the plot).

5.3. Cleaning method

The algorithm starts by identifying the pixel with the largest S/N and stores a unique identification number (ID) for this detection, the pixel position (Xpix, Ypix), the sky position (RA and Dec) in degrees θ_d, the signal-to-noise ratio (S/N) and the associated value of the amplitude map (A). Then, it proceeds by subtracting the signal contribution of this detection from the amplitude map to facilitate the detection of other possible nearby contributors. This is done by removing the filtered lensing signal carried by each individual galaxy that we model using the same shear template adopted for the filter construction. Given a detection centred at the position θ_d, the expected shear modulus at the location θ is

$$|{\gamma }_d(\theta ;{\theta }_d)|=A·\tau (\theta ;{\theta }_d),$$(30)

where the modelled shear profile τ(θ) is scaled by the detection amplitude A associated with the detection. The resulting complex shear attributed to the detection at each position is thus

$${\gamma }_d={\gamma }_{d,1}+i{\gamma }_{d,2}=|{\gamma }_d|(cos2\phi +isin2\phi ),$$(31)

which is tangential with respect to the detection location by construction.

Based on that, the two components of the tangential shear are

$$\begin{array}{cc}\hfill {\gamma }_{d,1}(\theta ;\mathrm{\Delta }\theta )& =-|{\gamma }_d(\theta ;\mathrm{\Delta }\theta )|·cos\left(2arctan\frac{\mathrm{\Delta }{\theta }_1}{\mathrm{\Delta }{\theta }_2}\right)\hfill \\ \hfill {\gamma }_{d,2}(\theta ;\mathrm{\Delta }\theta )& =-|{\gamma }_d(\theta ;\mathrm{\Delta }\theta )|·sin\left(2arctan\frac{\mathrm{\Delta }{\theta }_1}{\mathrm{\Delta }{\theta }_2}\right),\hfill \end{array}$$(32)

where, in the flat sky approximation, we can write Δθ=θ−θ_d=Δθ₁+iΔθ₂, with Δθ₁ and Δθ₂ being the difference in position of the galaxy and detection with respect to the two Cartesian coordinates. With this as the new shear components of the galaxy, the code computes the lensing contribution and removes it from the pixels of the map in the region surrounding the detections.

After the removal of the detection contribution, that is, by applying the described procedure to all galaxies, the S/N is re-evaluated, and the algorithm continues by considering the pixel with the next largest S/N, reiterating the previous process. It may happen that, after cleaning, some peaks that initially had a S/N value below the detection threshold can be boosted to be detectable. Therefore, detections are not identified with a strict monotonic decrease in their value of S/N.

A visual impression of the cleaning procedure is shown in Fig. 7, where four peaks are progressively detected and removed with the discussed procedure. When the algorithm reaches the minimum S/N threshold imposed by the user (for instance, the one based on the method described in Sect. 5.2), the procedure ends, and the iteration stops. Together with the catalogue of detections, the code outputs the amplitude map where all the lensing signals from the dark matter halo candidates identified during the run have been removed. Even though the cleaning procedure improves the number of detections by only a few percentage points (2% in the COMPLETE run) at low S/N, it ensures a more precise identification of peak positions, particularly in cases of blending, leading to a more robust matching. Additionally, the cleaned map can be used to assess the efficiency of the chosen S/N threshold and, more generally, to evaluate the performance of the detection procedure by comparing it to the B-mode map.

Fig. 7. Example of the cleaning procedure on a 0.36 deg2 region containing 4 signal peaks, indicated as a, b, c and d. Each circle is centred on the pixel with the maximum amplitude. The corresponding S/N is also reported, as ‘SNR’. Top panel: Amplitude map before cleaning. Smaller bottom panels: Amplitude map after each of the sequential cleaning steps. The algorithm detects and removes each peak in descending S/N order: first peak a (top-left), then peak b (top-right), peak c (bottom-left), and finally peak d (bottom-right), leaving at the end a map without any significant signal.

6. Results and statistical analysis

We now describe the results of the application of the AMICO-WL algorithm with the four filters discussed in Sect. 4.2 to the simulated catalogue of galaxies’ ellipticities described in Sect. 3, with the four different redshift cut-offs. From a visual inspection, the resulting E-mode, B-mode, and variance maps of these four cases appear similar. Still, the differences cannot be neglected, especially in the low S/N tail, because the application of the filters shown in Fig. 4 produces maps with different values. Moreover, the normalisation of the maps changes as the effective number density of sources n_g is different for the four catalogues of galaxies: cutting out foreground galaxies leads to a decrease in the surface number density from n_g = 30.0 gal/arcmin² for COMPLETE to n_g = 20.7 gal/arcmin² for FR06. Last but not least, a significant variation in the variance map comes from the change of ${\sigma }_{\mathrm{LSS}}^2$ that depends on the specific filter Ψ, and on the cosmic shear power spectrum P_γ, which is related to the redshift distribution of sources (see Eq. (28)).

For the four applications, we report in the upper part of Table 3 the resulting S/N thresholds corresponding to an expected purity of 70% as obtained with the method described in Sect. 5.3. The table also shows the corresponding number of detections.

In Fig. 8 we show the S/N distribution of the detections. In the four cases, most of the detections have S/N<4.0, namely 76% for COMPLETE, 83% for FR02, 79% for FR04 and 67% for FR06. In the high-end S/N tail, we found the same peaks in all four realisations but with different values of S/N.

Fig. 8. Distribution of signal-to-noise ratio S/N in the catalogues of detections in ΔS/N = 2.0 bins: COMPLETE detections in black, FR02 in green, FR04 in light blue and FR06 in pink bars. The counts are on a logarithmic scale.

6.1. Lens plane matching procedure

To assess the performance of AMICO-WL and compare the results obtained with the four different implementations of the optimal filter, it is necessary to cross-match the detections with the haloes in the simulations and distinguish the matched detections from those caused by noise fluctuations and contribution due to the large-scale structures. The most common procedure to perform this matching is to use a fixed physical, comoving, or angular radius centred on the detection and associate it to the closest halo inside the radius (see e.g. Hennawi & Spergel 2005; Gavazzi & Soucail 2007; Miyazaki et al. 2018; Hamana et al. 2020; Oguri et al. 2021). Even using an adaptive matching radius, the method is afflicted by the lack of redshift information on the weak-gravitational lensing side, leading to a high probability of including random matches due to the large number of haloes displaced along the line of sight of a detection. Moreover, this approach does not allow us to identify the actual halo responsible for the largest lensing contribution (i.e. the one causing the detection) and does not allow us to distinguish the cases in which a detection is caused by a series of small haloes displaced along the line of sight, which would not be detected individually.

For these reasons, we followed a method similar to Pace et al. (2007) to gain some redshift information regarding the weak lensing detections when performing their cross-matching with the haloes in the simulations. This approach consisted of a ‘blinking’ procedure exploiting the fact that haloes are well-localised structures, i.e. their signal originates from a single thin lens plane. In contrast, observational noise, uncorrelated structures, and projection effects originate from all redshift planes. Our matching method relied on the application of AMICO-WL to 10 different realisations built from the same simulation, but in which we removed from the ray-tracing the contribution of all haloes contained in redshift slices of Δz = 0.1, one at a time (see Sect. 3 for a discussion on the effect on the convergence maps and Figs. 1 and 3 for a visual impression). The thickness of the redshift slices was chosen for practical reasons: they are sufficiently thin for the purpose but large enough to keep the computational costs small. In this way, the lensing signal of a detection caused by a halo contained in the removed redshift slice is completely removed. In contrast, false detections due to noise, uncorrelated structures, or projection effects are only marginally affected by the removal, as their contribution comes from more redshift planes. This allowed us to gain information about the redshift of detections, thus improving the reliability of the cross-matching with the halo catalogue.

Pace et al. (2007) first matched the detections with the catalogue of haloes (over the entire redshift range) with a simple positional match and then removed the redshift plane that contained the matched halo. If the S/N drops by more than 25%, the candidate was considered to be connected with the single dark matter halo in the plane.

In this work, we used a different criterion based on a two-step procedure. First, we identified in the original catalogue the detections that disappear from at least one of the 10 realisations where a redshift slice has been removed: we will refer to them as ‘vanished detections’ from now on. This was done using an angular match between the original catalogue of detections and the 10 realisations adopting a matching radius of 0.6 arcmin to find the ones that have no counterparts. The radius has been chosen, studying the typical extension of the peaks in the amplitude map, to be large enough to allow possible positional fluctuations of the maxima of the “lensing blob” but small enough to reduce most of the fortuitous matches between close or even blended detections. As an example, Fig. 9 shows a 0.5 deg² zoom-in of the E-mode map resulting from the analysis of the simulation containing all lensing planes (top panel) and the one where the z_l∈[0.2,0.3] slice has been excluded (bottom panel). The large peak visible in the top panel, highlighted by the white circle, disappears completely from the map in the bottom panel, where it is therefore not detected by AMICO-WL. In contrast, all other structures, due to noise and LSS, are basically unaffected. With this first step, we identified the peaks whose lensing signal is well localised, and thus it was possible to assign a redshift estimate to the vanished detections given by the redshift range of the plane, which, when removed, causes the signal to drop.

Fig. 9. E-mode lensing signal map for a 0.5 deg2 region, extracted from the COMPLETE run. Top panel: E-mode map including all lens planes. Bottom panel: E-mode map created by excluding the zl∈[0.2,0.3] lens plane. The map was created excluding the plane which contains the cluster responsible for the highest peak, whose position is identified by the white circle in both panels.

The second step consisted of matching the vanished detections in each realisation with the haloes present in the corresponding removed redshift slice. By doing so, the matching criterion was extended to three dimensions with a Δz = 0.1 resolution in redshift. For this second step, we chose an angular radius around the haloes of 3 arcmin. The choice of the radius was still based on the typical size expected for the weak lensing imprints in the data but, differently from the first step, we relaxed the size to account for the discretised positions of detections in the amplitude maps and for the presence of a possible shift of the position of the peaks coming from the addition of the weak lensing signal of the LSS present along the line of sight. Let us note that a larger matching radius would create an artificial increase in successful matching due to chance alignments. The matching algorithm also implemented an a priori sorting of the halo catalogue by the virial mass and of the detections by S/N: high-mass haloes were first matched with higher S/N peaks inside the matching radius. This step produced the ‘real detections’ catalogue.

When matching was performed, we considered haloes with a virial mass M₂₀₀>5×10¹³ h⁻¹ M_⊙ and redshift z<1.0. This selection was motivated by studies on the modelling of the expected S/N of weak lensing surveys, which showed that smaller and/or very distant haloes cannot be detected solely with weak gravitational lensing (see e.g. Maturi et al. 2005; Pace et al. 2007; Andreon & Bergé 2012). Furthermore, since the number of haloes rises steeply at lower masses and higher redshifts, the matching would be increasingly prone to random associations. Including haloes with lower masses or at higher redshifts would significantly increase the number of potential matches, but most of these would fall below the detection limit of our method, leading to a higher rate of spurious or random associations.

Inside the 5 × 5 deg field of view, the selected catalogues for the 10 different lens slices, from z_l∈[0.0,0.1] to z_l∈[0.9,1.0] contained 3, 15, 42, 117, 82, 105, 151, 150, 144, 154 dark matter haloes, respectively. This redshift distribution is also shown as a filled orange histogram in the top panel of Fig. 10, while the corresponding virial mass distribution is displayed in the lower panel.

Fig. 10. Redshift (top panel) and virial mass (bottom panel) distributions of the simulated matter haloes of the reference catalogue used for the matching procedure (orange histogram). The corresponding distributions of the matched haloes in the COMPLETE run are shown as blue solid line.

Note that, in some cases, we associated the same detection with more than one cluster when its S/N falls below the detection threshold for more than one removed lens slice. This can happen when (1) multiple haloes along the line of sight contribute cumulatively to a single lensing peak. Removing separately the planes that host these haloes reduces each time the peak signal to below the detection threshold, leading to multiple associations. Alternatively, (2) when the S/N of the peak is close to the threshold, then even a possible small variation in the signal given by noise fluctuations or by the removal of large-scale structures along the line of sight can bring the peak below the detection threshold several times so that it is associated with multiple lens planes. The opposite case in which none of the removed planes results in a non-detection, falls into the category of spurious detections and it is discussed in Sect. 7. To define unique associations, we associated each detection with the halo with the highest virial mass within the multiple associations to obtain the final catalogue of detected clusters. In the COMPLETE case, fewer than 7% of the detections have multiple associations.

Figure 11 illustrates how the ‘blinking’ procedure confirms the expectations of lensing theory. Specifically, we show, as a function of the lens plane redshift and for the COMPLETE run only, the fraction of primary detections which are vanished or matched (dashed and solid lines, respectively). These fractions are not constant but peak at z_l∈[0.3,0.4], consistent with lensing theory, which predicts a stronger lensing signal (and thus a higher detection rate) when the lens is located at an intermediate position between the observer and the sources.

Fig. 11. Fraction (in percentage) of detections in the primary catalogue which are vanished (dashed line) and matched (solid line) as a function of the removed lens-plane. The results refer to the COMPLETE run.

In Fig. 12, we show the distribution of the scatter in position between the detections and the corresponding matched haloes. Even if the matching radius was set to 3 arcmin, the distribution drops well below this value, demonstrating that the number of random matches is negligible, as it would otherwise have produced a produced a distribution increasing with the separation, since the probability of random matches is proportional to the annulus area, where P(random match)∝annulus area = 2πrdr. This further shows the validity of the adopted matching approach.

Fig. 12. Distribution of the scatter in position between the halo and the matched detection in the real detections catalogue. The red dashed vertical line represents the size of the matching radius of 3 arcmin used in the second step of the lens plane matching procedure.

The 963 simulated haloes with M₂₀₀>5×10¹³ h⁻¹ M_⊙ were matched by the following number of AMICO-WL detections:

• COMPLETE run: 131 dark matter haloes, 9 of which have a double or triple association;

• FR02 run: 132 dark matter haloes, 7 of which have a double association;

• FR04 run: 141 dark matter haloes, 3 of which have a double association;

• FR06 run: 129 dark matter haloes, 6 of which have a double association.

Detections matching two or three haloes along the line of sight usually occur for S/N≲4.0, falling into the category of peaks near the detection threshold. The distribution of the matched haloes with respect to the total catalogue is reported in Fig. 10 as a function of redshift (top panel) and virial mass (bottom panel) for the COMPLETE and FR02 runs. Note that the bulk of the dark matter haloes’ distribution is at low mass and high redshift values, whereas weak lensing detection has a low efficiency. This limitation emphasises the need for complementary data to enhance the sensitivity of weak lensing detections, especially in challenging low-mass, high-redshift environments.

6.2. Purity and completeness

In this section, we describe the results of the statistical analysis performed using the lens-plane-removed matching procedure to the applications of the four filters in AMICO-WL discussed in Sect. 6.1. In Table 3, we report, for the four total samples we produced, the values of the sample completeness, defined as

$$\text{COMP}=\frac{\#\hspace{0.5em}\mathrm{matched}\hspace{0.5em}\mathrm{haloes}}{\#\hspace{0.5em}\mathrm{selected}\hspace{0.5em}\mathrm{reference}\hspace{0.5em}\mathrm{haloes}},$$(33)

and the purity, defined as

$$\text{PUR}=\frac{\#\hspace{0.5em}\mathrm{real}\hspace{0.5em}\mathrm{detections}}{\#\hspace{0.5em}\mathrm{primary}\hspace{0.5em}\mathrm{detections}}·$$(34)

The FR04 and especially the FR02 samples display a higher purity with respect to the COMPLETE one, demonstrating the efficiency of the simple foreground removal we implemented. For the FR02 sample, we found that the purity aligns with the estimates provided in Sect. 5.2 for the S/N threshold, supporting the validity of the selection procedure.

We also list the S/N thresholds leading to a sample purity of 70% and 80% with the corresponding number of detections and completeness. Such a comparison at a fixed sample purity allowed for a more coherent comparison of the sample completeness achieved with the different applications. Considering a 70% purity, the FR02 run delivered the highest completeness (13% against 6.5% of the COMPLETE run), while at 80% purity, the performance of the four applications was almost equivalent. This is because the selected peaks were those with higher S/N and were related to the haloes with higher virial mass and redshift in the high lensing efficiency range.

In Fig. 13 we summarise with a ROC curve the performance of AMICO-WL in terms of purity and completeness as a function of S/N. As expected, at low S/N more objects were found (high completeness) but with a high probability of detecting spurious sources (low purity). Vice versa, at high S/N, only a few (low completeness) but very reliable detections were found (high purity). The runs with foreground removal had a higher completeness than the COMPLETE run for the same purity level. The best performances were delivered by FR02.

Fig. 13. Completeness as a function of purity for different S/N thresholds. Different lines refer to the COMPLETE run (black solid line), FR02 run (green dashed line), FR04 run (cyan dash-dotted line) and FR06 run (pink dotted line). In the case of COMPLETE and FR02, the S/N threshold at which the purity and completeness are computed is noted for each point.

The completeness of the sample was not uniquely defined, as it depends on the mass cut-off used to define the reference true cluster sample, which can be chosen arbitrarily. It was, therefore, of great interest to investigate the completeness as a function of both the true redshift and the virial mass M₂₀₀ of the matched clusters. This is shown in Fig. 14 for the COMPLETE and FR02 samples considering a 70% purity level. In the same figure, the side plots show the completeness as a function of redshift (top panels) and mass (right panels).

Fig. 14. The sample completeness in two-dimensional bins of redshift and virial mass of the matched haloes for the COMPLETE (left panel) and FR02 (right panel) runs. Both samples are cut at the S/N threshold corresponding to a 70% purity. The colour map and the values reported inside the two-dimensional bins both refer to completeness. The side plots show the completeness as a function of redshift (on the top) and mass (on the right).

In both cases, the completeness is maximum, that is 100%, in the top left part of the plots that refer to the haloes with high mass (M₂₀₀>10^14.4 h⁻¹ M_⊙) and redshift in the range z=[0.2,0.4]. This was expected because these are the haloes with the stronger lensing effect: large mass and redshift roughly in between the background source and the observer. The completeness of the FR02 run was higher in all bins, especially for those at lower redshift. At the lowest redshifts, z=[0.0,0.1], we typically found only one or two haloes per mass bin. Without foreground removal, these detections fell below the 70% purity threshold we set for the plot and are therefore excluded, resulting in a completeness of 0%. However, the increase in S/N provided by the foreground removal in the FR02 run put them over the purity threshold, then we observed the jump in completeness to 100%. We refer the reader to the following section, where we analyse the effect of foreground removal also on the spurious detections.

7. Analysis of the spurious detections

In this section, we investigate the nature of spurious detections to gain an in-depth understanding of the cause of their origin and, at the same time, to check for possible failures in the detection method or in the matching procedure. Here we focus on the results of the COMPLETE run to assess the spurious detections within the standard AMICO-WL run and to further check the efficiency of the foreground removal in avoiding them.

A spurious detection is a galaxy cluster candidate identified by AMICO-WL that did not find a counterpart in the mock catalogue of dark matter haloes through the lens-plane matching procedure. Here we study whether they are noise peaks mistakenly identified as real or false negatives due to failures in the matching procedure.

We identified two main categories of spurious detections: the ‘non-vanished’ ones, which were always found regardless of the removal of any of the 10 redshift slices of matter, and the ‘non-matched’ ones, which vanished when one or more lens planes were removed from the ray-tracing but that do not match any dark matter halo in those redshift slices. Of the 82 total spurious detections in the COMPLETE run, 17 were non-vanished and 65 non-matched. Their distribution with respect to S/N is shown in Fig. 15. Half of the spurious detections have small signal-to-noise ratios, S/N<4.0, but several, especially the non-matched ones, have relatively high S/N up to approximately S/N≈6.0.

Fig. 15. Fraction of spurious detections as a function of S/N, in the COMPLETE run. The blue dashed line refers to the fraction of the non-vanished spurious, while the blue solid one refers to the total spurious (non-vanished plus non-matched).

The non-vanished spurious were produced by non-compact physical matter distribution or noise fluctuations covering several redshift slices along the line of sight. The latter was the main cause of their origin, as most of them did not have a corresponding convergence peak found in the noiseless maps, confirming their ”non-physical” origin.

This was not the case for the two non-vanished spurious detections with the highest S/N values, S/N = 4.88 and S/N = 4.31. Their S/N value dropped significantly when removing certain redshift slices, but always remained above the detection threshold. We found massive dark matter structures along their line of sight in the redshift plane in which the S/N decreased. They were, therefore, caused by projection effects due to interlopers. To account for similar cases, we could relax the selection criterion of the detections performed in the first step of the matching procedure by labelling as also vanished the detections that have more than a 3σ variation of S/N evaluated with a K-sigma clipping excluding outliers within the 10 values of S/N. The result is that the vanished detections increased for every lens plane, thus reducing the number of non-vanished spurious. Still, the number of real detections did not increase significantly, and the increase in purity is of ∼2%. We conclude that the negligible increase in the performance does not justify the introduction of this specialised refinement for the identification of vanished detection. On the other hand, when cross-checking these two detections in the FR02 run, which is the most effective foreground removal application, we found that the first detection remains spurious, while the second one had a positive match. In total, for the FR02 run, the number of non-vanished spurious decreased from 17 to 6: this shows that reducing the foreground contribution was also effective in decreasing spurious detections.

To investigate the origin of the non-matched spurious detections, we studied the distribution of the total mass of the haloes with M₂₀₀>10¹² h⁻¹ M_⊙ laying along the line of sight and within the redshift slice associated with each spurious detection. We thus found that these detections were generated by “compact chains” of low-mass haloes aligned along the line of sight with a significant total mass. This analysis revealed additional types of spurious detections, primarily those associated with the redshift plane z_l∈[0.0,0.1]. This redshift plane was sparsely populated by isolated haloes, primarily of mass which is lower than those in the selected reference catalogue used for matching. Despite their low mass, these haloes could still generate S/N peaks comparable to those found in higher-redshift planes. The line of sight analysis was also complemented with the addition of the ‘sub-haloes’. The sub-haloes are bound structures identified by SUBFIND present within the virialised region of dark matter haloes (Springel et al. 2001; Dolag et al. 2009; Giocoli et al. 2010). We found that a few significant spurious detections arise from sub-haloes excluded from the reference halo catalogue. Including sub-haloes in the reference catalogue could marginally increase completeness by fewer than 1%. However, this inclusion would also significantly increase the complexity of the halo selection process.

To summarise, we identified four main typologies associated with non-matched spurious detections: (i) a substantial group without the presence of other structures along the line of sight, where S/N fluctuations appeared solely due to the extra contribution of observational noise and LSS; (ii) several detections associated with a collection of aligned haloes and/or sub-haloes of smaller masses, which, due to the nature of gravitational lensing, can emulate the effect of a single more massive structure, resulting in false positives; (iii) detections vanishing in the closest redshift plane, z_l∈[0.0,0.1], in some case with high S/N, which are produced by dark matter haloes with masses M₂₀₀<5×10¹³ M_⊙/h; (iv) a final group of spurious detections arising exclusively from the most massive sub-haloes, which were excluded from the main halo catalogue used for the matching and thus remain unmatched. The spurious detections in groups (i) and (ii) are related to the nature of the data and the superposition principle holding for weak lensing. In the FR02 run, the number of ‘non-matched’ spurious detections decreased from 65 to 42, mainly due to the noise reduction that acted on the group (i). In contrast, the spurious groups (iii) and (iv) are false negatives given by the selection criteria employed for the reference catalogue used for matching.

8. Conclusions

The advent of upcoming wide-field surveys, such as Euclid, will provide an unprecedented wealth of data for studying galaxy clusters through weak lensing. However, effectively extracting cluster information from these massive datasets requires robust and accurate tools for data analysis. In this work, we introduced AMICO-WL, an extension of the well-tested AMICO algorithm, that allows the detection of galaxy clusters through the weak-gravitational lensing effect rather than the properties of the cluster's galaxies exploited in the original version of the code. The core of the algorithm is the linear optimal matched filter introduced by Maturi et al. (2005), to which we have added several improvements.

We evaluated the performance of AMICO-WL on a 25 deg² field mimicking the expected distribution of galaxies in a Euclid-like survey with the simulated shear measurements produced with DUSTGRAIN-pathfinder. AMICO-WL was applied to the field in four different configurations with optimal filters tailored to different redshift cuts of the input catalogue to reduce foreground contamination: COMPLETE, FR02, FR04, FR06 where FR stands for “foreground removal” and the digits represents the redshift cuts applied (z_min = 0.2, z_min = 0.4 and z_min = 0.6).

We used an accurate matching procedure based on the ‘blinking’ of individual redshift planes of the matter simulation to gain redshift information (otherwise missed in weak-lensing analysis) and analysed the code performance in terms of purity and completeness. We then performed a detailed statistical analysis of the spurious detections that are not associated with real sources to better understand their origin.

These are the main results:

• AMICO-WL applies the robust and well-tested optimal matched filtering technique, originally developed for photometric galaxy cluster detection. Its efficiency is enhanced through the AMICO's computational infrastructure and parallelisation. AMICO-WL is designed to be simple and user-friendly, making it well-suited for application to large datasets while producing reliable results. Future improvements may include a dynamical filtering mechanism to vary the sensitivity of the filter for different cluster features (e.g. mass, redshift).

• The detection process includes a cleaning procedure that removes from the E-mode amplitude map the lensing signal associated with each detected candidate cluster. This is achieved by subtracting the filtered model signal of the detection at each galaxy position using the same template adopted for the filter construction. The cleaning ensures that overlapping or nearby peaks do not affect subsequent detections, improving the robustness of the identification process.

• The matching procedure identifies halo-detection associations by analysing their response to the removal of individual redshift slices, thus adding redshift information to the weak lensing matching, which allows a reduction of spurious associations. For the COMPLETE run, 963 simulated haloes with M₂₀₀>5×10¹³ h⁻¹ M_⊙ are used as a reference, of which 131 objects are matched. This includes 9 detections associated with multiple haloes located along the line of sight. The fraction of matched detections peaks at z_l∈[0.3,0.4], where the lensing efficiency is maximised. We obtained a purity of 61% for COMPLETE and 71%, 64% and 61% for FR02, FR04 and FR06, respectively.

• The foreground removal technique, based on redshift cuts (with z_min = 0.2,0.4,0.6), reduces the noise coming from low-redshift galaxies, which act as foreground or contribute with very little lensing signal. This foreground removal significantly enhances detection reliability, improving both purity and completeness. For instance:

  • Considering a 70% purity level, completeness increases from 6.5% for the COMPLETE run to 13.0% for the FR02 run.

  • Considering a more conservative 80% purity level, the completeness values for all runs converge since only high-S/N detections related to the most massive haloes remain in the sample. In particular, 5.3% for COMPLETE and 5.8% for FR02.

  The FR02 run demonstrates higher efficiency in detecting lower-mass haloes compared to the COMPLETE run at similar purity thresholds. Furthermore, the spurious detection rate in the FR02 run is reduced by 36% compared to the COMPLETE run. A dynamic redshift-dependent foreground removal strategy could further improve both purity and completeness.

• Spurious detections are categorised into ‘non-vanished’ detections, associated with noise or filamentary structures, and ‘non-matched’ detections, often arising from alignments of low-mass haloes along the line of sight. Of the 82 spurious detections in the COMPLETE run:

  • 17 are non-vanished and persist regardless of lens-plane removal.

  • 65 are non-matched, largely attributed to alignments of low-mass haloes emulating massive structures.

  From this analysis of the spurious detections, we quantified the power and limits of both the AMICO-WL detection method and the lens plane matching. We suggested some mitigation strategies to reduce the false negatives: refinement of the matching procedure and of the halo catalogue selection. The implementation of these intricate and specialised mitigation measures for false-negative spurious detections is not justified in the general analysis because of the negligible gain in purity they achieve.

In conclusion, AMICO-WL shows its effectiveness as a weak-lensing cluster detection tool, and the matching procedure that has been put in place has proved crucial in determining the characteristics of the structures that have been found. AMICO-WL is currently being evaluated in comparison with other algorithms in the Euclid Weak Lensing Cluster Challenge, demonstrating its versatility and competitiveness. We plan to improve the code by implementing an adaptive filter that would increase the flexibility of the algorithm and allow the detection of dark matter haloes in a broader range of masses and redshift. The ultimate goal is to build an optimised routine for weak-lensing detection of clusters complementary to those developed for photometric detection. The combination of the parent AMICO photometric algorithm with the new AMICO-WL represents an extremely versatile and powerful tool for detecting clusters that have different characteristics. The next step will be the application of this code to real observational data coming from existing or upcoming wide-field surveys, such as KiDS, CHFTLens, DES, Euclid, and LSST⁴ from the Vera C. Rubin Observatory.

Acknowledgments

We thank the anonymous referee for the careful reading and the useful comments which helped to improve the quality of the paper. LM and CG acknowledge the financial contribution from the PRIN-MUR 2022 20227RNLY3 grant ‘The concordance cosmological model: stress-tests with galaxy clusters’ supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 “Attivit scientifiche per la missione Euclid – fase E”. GC thanks the support from INAF theory Grant 2022: Illuminating Dark Matter using Weak Lensing by Cluster Satellites.

References

All Tables

All Figures

Fig. 1. Convergence map (5 deg on a side) for fixed source redshift zs = 1. In the left panel, we consider all the matter present along the line of sight from z = 0 to z = 1, while in the central and right panels, we remove the matter located in a redshift slice as indicated in the legend. The red and yellow squares highlight some particular high-density regions that disappear in the central and right panels when discarding the matter present in the considered redshift slices.

In the text

	Fig. 2. Redshift distribution of the weak lensing simulated sources. The black empty histogram refers to the COMPLETE catalogue, while the hatched histograms represent the foreground-removed catalogues FR02, FR04 and FR06 (described in Sect. 4.2), shown as green, light-blue and magenta, respectively.
In the text

	Fig. 3. Reconstructed convergence map from the source redshift distribution catalogue. The left panel refers to the full light-cone up to z = 4, while in the central and right panels, we show the maps where we exclude the matter located between z = 0.2−0.3 and z = 0.4−0.5, respectively.
In the text

	Fig. 4. Radial profiles (in real space) of the optimal filters employed for the four AMICO-WL runs: COMPLETE (black solid line), FR02 (green dashed line), FR04 (cyan dash-dotted line) and FR06 (pink dotted line).
In the text

	Fig. 5. The E-mode amplitude and its variance lensing signal maps are shown in the top-left and bottom-left panels, respectively, while the B-mode amplitude and its variance lensing signal maps are shown in the top-right and bottom-right panels. The maps cover a field of view 5° × 5° and refer to the COMPLETE run.
In the text

Fig. 6. S/N threshold determination analysis. Top panel: Cumulative S/N distributions of positive (ES/N) and absolute negative peaks (${\hat{E}}_{S/N}$) in the E-mode map detected with SExtractor. Bottom panel: Difference in cumulative counts of positive and negative peaks as a function of the S/N. The S/N threshold achieving a 70% purity is defined where the cumulative positive peak count reaches 70% of the total (horizontal line) and is marked across both panels (red vertical dashed line).

In the text

Fig. 7. Example of the cleaning procedure on a 0.36 deg2 region containing 4 signal peaks, indicated as a, b, c and d. Each circle is centred on the pixel with the maximum amplitude. The corresponding S/N is also reported, as ‘SNR’. Top panel: Amplitude map before cleaning. Smaller bottom panels: Amplitude map after each of the sequential cleaning steps. The algorithm detects and removes each peak in descending S/N order: first peak a (top-left), then peak b (top-right), peak c (bottom-left), and finally peak d (bottom-right), leaving at the end a map without any significant signal.

In the text

	Fig. 8. Distribution of signal-to-noise ratio S/N in the catalogues of detections in ΔS/N = 2.0 bins: COMPLETE detections in black, FR02 in green, FR04 in light blue and FR06 in pink bars. The counts are on a logarithmic scale.
In the text

	Fig. 9. E-mode lensing signal map for a 0.5 deg2 region, extracted from the COMPLETE run. Top panel: E-mode map including all lens planes. Bottom panel: E-mode map created by excluding the zl∈[0.2,0.3] lens plane. The map was created excluding the plane which contains the cluster responsible for the highest peak, whose position is identified by the white circle in both panels.
In the text

	Fig. 10. Redshift (top panel) and virial mass (bottom panel) distributions of the simulated matter haloes of the reference catalogue used for the matching procedure (orange histogram). The corresponding distributions of the matched haloes in the COMPLETE run are shown as blue solid line.
In the text

	Fig. 11. Fraction (in percentage) of detections in the primary catalogue which are vanished (dashed line) and matched (solid line) as a function of the removed lens-plane. The results refer to the COMPLETE run.
In the text

	Fig. 12. Distribution of the scatter in position between the halo and the matched detection in the real detections catalogue. The red dashed vertical line represents the size of the matching radius of 3 arcmin used in the second step of the lens plane matching procedure.
In the text

	Fig. 13. Completeness as a function of purity for different S/N thresholds. Different lines refer to the COMPLETE run (black solid line), FR02 run (green dashed line), FR04 run (cyan dash-dotted line) and FR06 run (pink dotted line). In the case of COMPLETE and FR02, the S/N threshold at which the purity and completeness are computed is noted for each point.
In the text

Fig. 14. The sample completeness in two-dimensional bins of redshift and virial mass of the matched haloes for the COMPLETE (left panel) and FR02 (right panel) runs. Both samples are cut at the S/N threshold corresponding to a 70% purity. The colour map and the values reported inside the two-dimensional bins both refer to completeness. The side plots show the completeness as a function of redshift (on the top) and mass (on the right).

In the text

	Fig. 15. Fraction of spurious detections as a function of S/N, in the COMPLETE run. The blue dashed line refers to the fraction of the non-vanished spurious, while the blue solid one refers to the total spurious (non-vanished plus non-matched).
In the text