© ESO 2022

1. Introduction

Statistics is an essential part of astronomy (Heck et al. 1985; Feigelson 1988, 2009; Feigelson & Babu 2013). The field relies on inferring physical properties, which cannot be determined directly, from observable quantities, which in turn need to be corrected for systematic effects as well as instrumental and observational biases. The key question that always needs to be answered when interpreting observations and results – before discussing how accurately these results can be constrained – is what one is actually seeing.

Weak gravitational lensing, caused by the deflection of light rays by density variations along the traveled path, has been a case in point for the last three decades. Gravitational lensing is a convex focusing effect that can magnify and shear affected background sources. The observed shapes and number counts can conversely yield information about these density variations but need to be disentangled statistically from the unknown intrinsic properties of background sources, such as distance, size (and luminosity), and shape.

The first detections of coherent alignments of galaxy shapes were observed in the background of clusters (Tyson et al. 1990), and subsequently in the emerging fields of galaxy-galaxy lensing (where the lensing “structure” is itself an ensemble of lenses; Brainerd et al. 1996) and cosmic shear (the weak lensing induced by large-scale structure; Wittman et al. 2000; Bacon et al. 2000; Kaiser et al. 2000; Van Waerbeke et al. 2000). Since then, techniques have progressed rapidly, and demands on accuracy have become increasingly stringent.

This is the second in a set of papers wherein we focus on the statistical aspects of inferring the lensing signal from the intrinsic shapes and the estimated lensing geometry, which depends on the distances between the observer, the moment of deflection, and the background sources. Assuming the cosmological principle, the intrinsic shapes of a sample¹ of background galaxies, including their orientation, are random, and the intrinsic galaxy shapes should average out from a sufficiently large sample, leaving the weak lensing signal as a net ellipticity. The common approach has been to take a weighted mean of galaxy ellipticities, which has computational and analytical advantages and, most importantly, is an unbiased estimator of the shear in the absence of pixel noise in the galaxy images (Seitz & Schneider 1997).

In practice, however, there are many sources of noise and the mean is known to be biased, underestimating the underlying shear signal (Melchior & Viola 2012; Viola et al. 2014; Sellentin et al. 2018; Mandelbaum 2018). The distribution of intrinsic galaxy shapes is well known to be non-Gaussian (Lambas et al. 1992; Rodríguez & Padilla 2013) and, in fact, centrally peaked. In Smit & Kuijken (2018, hereafter Paper I), we explored alternative estimators besides the mean that could potentially be better suited for such a cuspy distribution. It was found, using realistic simulated distributions and resampling of Canada-France-Hawaii Lensing Survey (CFHTLenS) shape measurements (Heymans et al. 2012b), that L¹ norm regression, also known as least absolute deviations (LAD), reduces bias from between ∼ − 4% and ∼ − 4.5% to between ∼ + 1% and ∼ − 3%, while at the same time reducing uncertainty by ∼9% to ∼23%.

In this paper we extend this study by applying these statistics to a weak lensing analysis of 6925 galaxy clusters in the Adaptive Matched Identifier of Clustered Objects (AMICO) cluster catalog (Bellagamba et al. 2011, 2018, 2019; Radovich et al. 2017; Maturi et al. 2019) of the third data release of the Kilo-Degree Survey (KiDS-450; de Jong et al. 2017). As opposed to Paper I, in this case the true lensing signal (here in the form of the excess surface density of the clusters) is unknown. We therefore study the relative biases and uncertainties between LAD and the mean, and we compare results to our findings in Paper I.

An important application is then to study the relation between the observable properties of clusters and groups and the physical quantities derived from the lensing signal (i.e., the matter distribution) to better our understanding of galaxy and cluster formation and cosmological models (e.g., Kautsch et al. 2008; Leauthaud et al. 2010; Lesci et al., in prep.). We calculate halo masses from the obtained lensing signals and derive a scaling relation between the observed r-band luminosity and the lensing mass, investigating the impact of estimator choice on the resulting constraints.

The order of magnitude of this estimated bias in the weak lensing results can be dominant compared to other sources of uncertainty in the process. Developments in the field have led to current constraints of the multiplicative bias in shape measurements on the order of ∼1% (Bernstein & Jarvis 2002; Hirata & Seljak 2003; Heymans et al. 2006, 2012a; Massey et al. 2007; Miller et al. 2007; Kitching et al. 2008; Bridle et al. 2010; Voigt & Bridle 2010; Bernstein 2010; Kitching et al. 2012; Kacprzak et al. 2012; Melchior & Viola 2012; Refregier et al. 2012; Mandelbaum et al. 2015; Viola et al. 2015; Fenech Conti et al. 2017). The uncertainty in the lensing geometry between the observer, lens, and background sources, introduced by the estimation of the photometric redshift probability distributions, can be a few percent (Hildebrandt et al. 2017; Bellagamba et al. 2019, and Appendix A.1). The broad category of selection biases, for example those introduced by intrinsic alignments, contamination of the background sample by cluster member galaxies, blending, detection, and subsequent selection effects, typically accumulate up to a few percent (Miyatake et al. 2015; van Uitert et al. 2017; Bellagamba et al. 2019) for cluster weak lensing. For instance, estimations on background selection yield a foreground contamination on the order of 2%, which can be partly corrected for, but does increase the uncertainty (Dvornik et al. 2017; Bellagamba et al. 2019, and Appendix A.2). In this study we investigate the usability of background sources to radii smaller than in Bellagamba et al. (2019).

These demands on accuracy and precision become higher as the data yield, and therefore the statistical power of surveys, increases dramatically (Mandelbaum 2018), as achieved by COSMOS² (Leauthaud et al. 2007), CFHTLenS³ (Heymans et al. 2012b), RCSLenS⁴ (Hildebrandt et al. 2016), KiDS⁵ (de Jong et al. 2013), and DES⁶ (Dark Energy Survey Collaboration 2016), and foreseen for future surveys such as LSST⁷ (Ivezić et al. 2019) and Euclid⁸ (Laureijs et al. 2011). While these two future surveys will require constraints on systematic uncertainty of order ≤2 × 10⁻³ (Mandelbaum 2018), we show that, even for weak lensing analyses in the last decade, the bias in shear inference can dominate other sources, such as the aforementioned multiplicative shape measurement bias that is commonly corrected for, as in Viola et al. (2015), Dvornik et al. (2017), and Bellagamba et al. (2019).

Several other approaches have been made to address this, including analytic modeling of the bias (e.g., Viola et al. 2014), weight corrections and priors (Bonnet & Mellier 1995; Van Waerbeke et al. 2000; Bernstein & Jarvis 2002), or nulling techniques (Herbonnet et al. 2017). The calculation of the main observable, the shapes of lensed background sources, itself relies on statistical methods. These are based mainly on surface brightness moments (Kaiser et al. 1995; Rhodes et al. 2000) or model fitting (Kuijken 1999, 2006; Bernstein & Jarvis 2002; Hirata & Seljak 2003; Refregier & Bacon 2003; Miller et al. 2007; Kitching et al. 2008). This means the most common approaches are corrections on a statistic that remains fundamentally skewed (Sellentin et al. 2018; Mandelbaum 2018).

Promising alternative approaches by Bernstein & Armstrong (2014) and Schneider et al. (2015) do not reproduce individual background shapes, but directly determine the underlying shear field from ensembles of background sources, reconsidering these steps in the chain of statistical inference. While future lensing surveys will require innovative improvements, these methods and their priors need to be gauged by deep observations of high signal-to-noise, and it is of fundamental importance that these calibrations are well constrained and do not suffer from even subtle systematic biases. In other words, the comparison of several perspectives is paramount in determining what we actually see.

The remainder of this paper is organized as follows. We introduce the definitions of galaxy shapes and the weak lensing formalism in Sect. 2 and relate these to our statistical approach. Data, analysis methods, and selection criteria are described in Sect. 3, while Sect. 4 states our results and analysis. Section 5 gives a summary of our conclusions.

Throughout this paper we assume a Planck (Planck Collaboration XVI 2014) cosmology with Ω_M = 0.315, Ω_Λ = 0.685, and H₀ = 100.0 h km s⁻¹ Mpc⁻¹. All measurements are in comoving units, unless specifically noted otherwise, such as in Sect. 2.

2. Weak gravitational lensing statistics

We briefly review the principles of weak gravitational lensing and relate the central concepts to our statistical approach, introducing the terminology and notation conventions used in this paper. We refer the reader to excellent reviews, such as Bartelmann & Schneider (2001), Schneider (2006), Hoekstra & Jain (2008), and Bartelmann & Maturi (2017), for more in-depth approaches.

2.1. Principles of weak lensing

Rays of light are deflected by the curvature or space-time due to mass inhomogeneities along their path. A mass overdensity acts as a convex lens on the light rays from distant sources behind that lens to an observer. In this section, we use D_l to denote angular-diameter distances from the observer to the lens, D_ls from the lens to the background source, and D_s from the observer to the background source (see Fig. 1), and in the remainder of this paper we translate quantities to comoving units where necessary.

Fig. 1. Representation of a gravitational lens system, showing the displacement of a source at position S(x, y) to an image at position I(x′,y′), where we take the origin of the source plane to be collinear with the position of the lens, L, and the observer, O.

For the purposes of this work, the extent of the lensing mass along the line of sight, compared to the distances from the observer to the lens and from the lens to the background source, can be considered negligible. In this so-called thin-lens approximation, the deflection of light rays by a deflection angle, $\widehat{\mathit{\alpha }}$, leads to an effective angular displacement (again, see Fig. 1),

$$\begin{array}{c}\hfill \mathit{\alpha }=-\frac{D_{\mathrm{ls}}}{D_{\mathrm{s}}}\widehat{\mathit{\alpha }},\end{array}$$(1)

also called the reduced deflection angle, which relates the observed position, θ, of a distant point source to its unlensed position, β, by the lens equation

$$\begin{array}{c}\hfill \mathit{\beta }=\mathit{\theta }-\mathit{\alpha }.\end{array}$$(2)

It can be shown through the relation between $\widehat{\mathit{\alpha }}$ and the three-dimensional gravitational potential, Φ, that this displacement is then described by α = ∇_θ ψ, where

$$\begin{array}{c}\hfill \psi =\frac{2}{c^2}\frac{D_{\mathrm{ls}}}{D_{\mathrm{l}}{D}_{\mathrm{s}}}{\displaystyle \int }\mathrm{\Phi }\mathrm{d}z\end{array}$$(3)

is called the (two-dimensional) lensing potential.

The differential effect of the deflection of light on the images, I(x, y), of extended background sources can to first order be described as a coordinate transformation by taking the derivatives in the lens equation (Eq. (2)) of the original angular position, β, with respect to the observed position, θ. Substituting ∇_θ ψ for α, we obtain the Jacobian matrix,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}1-{\psi }_{11}& -{\psi }_{12}\\ -{\psi }_{21}& 1-{\psi }_{22}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right),\end{array}$$(4)

with

$$\begin{array}{c}\hfill {\psi }_{\mathit{ij}}=\frac{{\partial }^2\psi }{\partial {\theta }_i\partial {\theta }_j},\end{array}$$(5)

resulting in the lensed image I(x′,y′), which is the key observable in our weak lensing study.

2.1.1. Critical surface mass density

To interpret the effect on the source image, we note that such a transformation can be decomposed into three parts, namely the identity (I), an isotropic part that describes a multiplication, and an anisotropic traceless part that describes a shearing of the image:

$$\begin{array}{c}\hfill \mathit{I}-\frac{1}{2}({\psi }_{11}+{\psi }_{22})\mathit{I}+\left(\begin{array}{cc}-\frac{1}{2}({\psi }_{11}-{\psi }_{22})& -{\psi }_{12}\\ -{\psi }_{21}& \frac{1}{2}({\psi }_{11}-{\psi }_{22})\end{array}\right).\end{array}$$(6)

To relate ψ_ij with the density of the lensing mass, we start with the isotropic term, which is half the Laplacian of the lensing potential: $\frac{1}{2}({\psi }_{11}+{\psi }_{22})=\frac{1}{2}{\mathrm{\nabla }}_{\theta }^2\psi$. From Eq. (3), we obtain

$$\begin{array}{c}\hfill \frac{1}{2}{\mathrm{\nabla }}_{\theta }^2\psi =\frac{1}{c^2}\frac{D_{\mathrm{l}}{D}_{\mathrm{ls}}}{D_{\mathrm{s}}}{\displaystyle \int }4\pi G\rho \mathrm{d}z,\end{array}$$(7)

which is a dimensionless quantity. Defining the surface mass density as

$$\begin{array}{c}\hfill \mathrm{\Sigma }\equiv {\displaystyle \int }\rho \mathrm{d}z\end{array}$$(8)

and gathering the rest of the right-hand side into

$$\begin{array}{c}\hfill \frac{4\pi G}{c^2}\frac{D_{\mathrm{l}}{D}_{\mathrm{ls}}}{D_{\mathrm{s}}}\equiv {\mathrm{\Sigma }}_{\mathrm{cr}}^{-1},\end{array}$$(9)

with Σ_cr being the critical surface mass density, we find that the isotropic term can be written as

$$\begin{array}{c}\hfill \kappa \equiv \frac{1}{2}{\mathrm{\nabla }}_{\theta }^2\psi =\frac{\mathrm{\Sigma }}{{\mathrm{\Sigma }}_{\mathrm{cr}}},\end{array}$$(10)

with κ a normalized dimensionless surface mass density. Recognizing that ${\nabla }_{\theta }^2\psi =\mathbf{\nabla }\cdot \mathit{\alpha }$ is the divergence of the deflection of the light rays (i.e., the manner in which those light rays converge due to the lensing effect), κ is simply the convergence.

2.1.2. Shear and intrinsic ellipticity

The shear matrix in Eq. (6) has two independent components, simply called the shear γ = γ₁ + iγ₂, with ${\gamma }_1=\frac{1}{2}({\psi }_{11}-{\psi }_{22})$ and γ₂ = ψ₁₂ = ψ₂₁. Equation (4) then becomes

$$\begin{array}{c}\hfill \left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2\\ -{\gamma }_2& 1-\kappa +{\gamma }_1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right).\end{array}$$(11)

This transformation leads to the magnification and distortion of the light distribution of background sources. In this work, we focus on the most commonly used net distortion or reduced shear g = g₁ + ig₂ ≡ (γ₁ + iγ₂)/(1 − κ),

$$\begin{array}{c}\hfill \left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=(1-\kappa )\left(\begin{array}{cc}1-g_1& -g_2\\ -g_2& 1+g_1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right),\end{array}$$(12)

where the transformation is written as a multiplication of (1 − κ) and a distortion matrix describing the alignment of lensed sources in the foreground potential.

The effect on a circular source is a shearing into an ellipse with axis ratio $q=\frac{b}{a}$ as

$$\begin{array}{c}\hfill q=\frac{1-|g|}{1+|g|}\iff |g|=\frac{1-q}{1+q}=\frac{a-b}{a+b}\end{array}$$(13)

and position angle φ via

$$\begin{array}{c}\hfill g=|g|(cos2\phi +isin2\phi ).\end{array}$$(14)

As mentioned before, we do not measure this gravitational distortion directly. Background sources have an intrinsic shape distribution, and we effectively measure the combined effect of their intrinsic shape and a weak lensing distortion. It is adequate to describe images by their quadrupole brightness moments or their ellipticities as well as the respective response to weak shear distortions. It is straightforward to use the common definition⁹ of ellipticity, defined as the reduced shear needed to create the intrinsic shape ϵ = ϵ₁ + iϵ₂ of a source from an image with circular isophotes (Bernstein & Jarvis 2002; Kuijken 2006). The resulting ellipticity, ϵ, after transforming an image with intrinsic¹⁰ ellipticity ϵ^I by a distortion, g, is then given by (Seitz & Schneider 1997)

$$\begin{array}{c}\hfill ϵ=\frac{{ϵ}^I+g}{1+g^{\ast }{ϵ}^I}\mathrm{for}|g|\le 1,\end{array}$$(15)

with g^* the complex conjugate of g.

The intrinsic shape distribution is called the shape noise and, assuming no preferred direction on the sky, should average to zero: ⟨ϵ^I⟩ = 0. This way, each background shape measurement, ϵ, is then an independent estimate of the underlying reduced shear, g.

In this paper we make use of the fact that the lensing signal is weak (i.e., κ ≪ 1) and assume g ≈ γ.

2.2. Estimation of the surface density profile

The shear induced by gravitational lensing is sensitive to the density contrast. For an axisymmetric lens, we can write $|\gamma |(R)=\overline{\kappa }(\le R)-\kappa (R)$, where $\overline{\kappa }$ is the average convergence within radius R. In fact, this relation holds for other mass distributions if we average azimuthally around the lens. In this work, we study the stacked signal of many lenses and assume a net axisymmetry (see, e.g., Evans & Bridle 2009; Oguri et al. 2010; Clampitt & Jain 2016; van Uitert et al. 2017, for weak lensing studies on elliptical lenses).

Since it can be seen from Fig. 1 that the gravitational shear acts in the radial direction, we define the tangential and cross components of the shear as

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\gamma }_{+}\\ {\gamma }_{\times }\end{array}\right)=\left(\begin{array}{cc}\hfill -cos(2\varphi )& \hfill -sin(2\varphi )\\ \hfill sin(2\varphi )& \hfill -cos(2\varphi )\end{array}\right)\left(\begin{array}{c}{\gamma }_1\\ {\gamma }_2\end{array}\right),\end{array}$$(16)

with ϕ the counterclockwise angle between the positive x axis¹¹ and the vector from lens to source. This gives

$$\begin{array}{c}\hfill \overline{\mathrm{\Sigma }}(\le R)-\mathrm{\Sigma }(R)\equiv \mathrm{\Delta }\mathrm{\Sigma }(R)={\gamma }_{+}(R){\mathrm{\Sigma }}_{\mathrm{cr}},\end{array}$$(17)

with ΔΣ(R) the excess surface density (ESD) at a radius R around the lensing mass. In axisymmetric lenses, the cross component of the shear cannot arise from gravitational lensing and should average to zero, if only produced by intrinsic source orientations, and can therefore be used as an indication of systematic effects, such as imperfect corrections for the point-spread function (PSF; Schneider 2003, and Appendix A.3).

The ESD is then estimated using the observed ellipticities of an ensemble of sources around the lens

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)=⟨{ϵ}_{+}{\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}⟩(R),\end{array}$$(18)

with each ϵ₊Σ_cr, ls an independent, albeit noisy, estimate of the ESD. Here, ⟨ ⋅ ⟩ denotes a weighted average, with weights to be specified.

The Σ_cr, ls behaves as a geometric scaling factor, indicating the lensing efficiency for each lens-source combination. Since the variance of the noise in ΔΣ is then affected by ${\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}^2$, the relative precision, or inverse variance, carried by each ϵ₊ scales as ${\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}^{-2}$.

In this paper we study the ESD profile in comoving radial bins, and we therefore use the comoving critical surface density

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{com}}={(1+z_{\mathrm{l}})}^2{\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{prop}}.\end{array}$$(19)

In practice, the distance to each background source is not known exactly and is estimated by its redshift probability distribution, p(z_s). Taking this into account, we estimate the comoving critical surface density via

$$\begin{array}{c}\hfill \langle {\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}^{-1}\rangle =\frac{4\pi G}{c^2}D(z_{\mathrm{l}}){(1+z_{\mathrm{l}})}^2{\displaystyle \int }\frac{D(z_{\mathrm{l}},z_{\mathrm{s}})}{D(z_{\mathrm{s}})}p(z_{\mathrm{s}})\mathrm{d}{z}_{\mathrm{s}}.\end{array}$$(20)

2.3. Statistical framework

In this section we discuss the estimation of ⟨ϵ₊ Σ_cr, ls⟩. We refer to Paper I for a complementary discussion.

Important aspects of a good estimator, $\widehat{ϵ}$, are: (i) minimal bias, defined as the difference between the expected value of the estimator, $\langle \widehat{ϵ}\rangle$, and the value of the quantity being estimated, for instance the shear (γ) or, in this case, ΔΣ; (ii) high efficiency, proportional to the inverse variance of the estimator, ${\sigma }_{\widehat{ϵ}}^{-2}$; and (iii) robustness, meaning the estimator retains these properties for a sufficient range of likely parameter distributions.

2.3.1. Bias

Even though the measured ellipticity, ϵ, is not a linear combination of the intrinsic shape, ϵ^I, and the shear, γ, it can be shown (Seitz & Schneider 1997) that, in the absence of further uncertainties, the mean μ(ϵ) is an unbiased estimator for the underlying shear and that this is independent of the intrinsic ellipticity distribution, P(ϵ^I). In the canonical approach, the ESD is therefore estimated as a weighted mean of an ensemble of lens-source combinations,

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)=\frac{{\sum }_{\mathrm{ls}}{w}_{\mathrm{ls}}{ϵ}_{+,\mathrm{ls}}{\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}}{{\sum }_{\mathrm{ls}}{w}_{\mathrm{ls}}},\end{array}$$(21)

where we use

$$\begin{array}{c}\hfill w_{\mathrm{ls}}=w_{\mathrm{s}}{\langle {\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}^{-1}\rangle }^2.\end{array}$$(22)

Here, the weight w_s is assigned to each measured ellipticity, scaled by the estimated lensing efficiency, as explained in the previous section (see, e.g., Viola et al. 2015; Dvornik et al. 2017; Bellagamba et al. 2019).

In practice, there are various sources of uncertainty at each step of the process, such as source selection bias, distortion by the PSF, and biases due to the measurement pipeline. These lead to convolutions of the ellipticity distribution, before and after the gravitational lensing effect. The result is a bias in the mean as an estimate of the ESD (Melchior & Viola 2012; Refregier et al. 2012; Kacprzak et al. 2012; Viola et al. 2014; Kacprzak et al. 2014). In this case, the intrinsic shape distribution will play a role.

The weighted mean, μ, is a statistic that, for a set of measurements ϵ_i with weights w_i, finds the estimate of γ that minimizes the loss function

$$\begin{array}{c}\hfill S_{\mu }={\displaystyle \sum _i{w}_i[{({ϵ}_{i,1}-{\gamma }_1)}^2+{({ϵ}_{i,2}-{\gamma }_2)}^2],}\end{array}$$(23)

that is, it is a least squares (LSQ) or L² norm regression method and arises naturally as the optimal estimator for Gaussian distributions.

Figure 2 shows that the measured ellipticity distribution, P(ϵ), displays crucial differences with a Gaussian distribution, showing a sharp peak and a slower decline, including a higher number of high ellipticities, |ϵ|. This central peak is an unbiased tracer of the underlying shear (Paper I). By Eq. (23), the mean is sensitive to outliers and therefore not robust when inferring the shear.

Fig. 2. Estimated ellipticity and axis ratio distributions of sources in the KiDS-450 catalog. Top: two-dimensional histogram of ellipticities. Middle: histogram of the absolute ellipticity, |ϵ|. Bottom: histogram of the ellipse axis ratio, q.

In contrast, LAD or L¹ norm regression minimizes the loss function,

$$\begin{array}{c}\hfill S_{\mathrm{LAD}}={\displaystyle \sum _i{w}_i\sqrt{{({ϵ}_{i,1}-{\gamma }_1)}^2+{({ϵ}_{i,2}-{\gamma }_2)}^2}.}\end{array}$$(24)

The LAD estimate is also known as the median in one dimension or the spatial median in higher dimensions. This estimator is more sensitive to the peak and less sensitive to high ellipticity outliers. Where the mean is expected to be biased low (Melchior & Viola 2012; Refregier et al. 2012; Kacprzak et al. 2012; Viola et al. 2014; Kacprzak et al. 2014), we expect this to be less so for the LAD (Paper I).

2.3.2. Efficiency

The formal definition of efficiency, $\stackrel{\sim }{\eta }$, relates the inherent (Fisher) information, ℐ, of a sample to the statistical variability around the expected value of the estimator, usually taken to be the variance, ${\sigma }_{\widehat{ϵ}}^2$, of the estimator:

$$\begin{array}{c}\hfill \stackrel{\sim }{\eta }=\frac{1}{\mathcal{I}·{\sigma }_{\widehat{ϵ}}^2}.\end{array}$$(25)

Since the variance of an unbiased estimator cannot be less than the reciprocal of the information, ${\mathcal{I}}^{-1}\le {\sigma }_{\widehat{ϵ}}^2$, we have $0\le \stackrel{\sim }{\eta }\le 1$ (Rao 1945; Cramer 1946).

As we are comparing two estimators with unknown bias, it is appropriate to use the relative efficiency,

$$\begin{array}{c}\hfill \eta =\frac{{\sigma }_{\mu }^2}{{\sigma }_{\mathrm{LAD}}^2},\end{array}$$(26)

with η < 1 indicating a higher efficiency for the mean, and vice versa.

Paper I showed that the LAD consistently performed better than the mean, with both higher efficiency and less bias, for various cusped intrinsic ellipticity distributions, including the shear catalog from CFHTLenS. In what follows, we take the CFHTLenS shear distribution shape to be representative of KiDS data as well since both surveys were processed with the THELI pipeline (Erben et al. 2013) and the shape measurement pipeline lensfit (Miller et al. 2007; Kitching et al. 2008; Miller et al. 2013).

2.4. Halo model

Studying the effects of estimator choice on the weak lensing signal forms the technical core of this paper. The scientific goal, however, is to assess the relevance on the inference of physical quantities, such as the derivation of a lensing halo mass from an ESD profile. Since we calculate the stacked signal for an ensemble of clusters with some common (observable) property (here a range in r-band luminosity), of interest is the scaling relation between the observable and derived lensing mass, M₂₀₀, where we use the definition with respect to the mean density of the universe.

To do so, we modeled the lens density profile the same way as Dvornik et al. (2017), using the halo model (Seljak 2000; Peacock & Smith 2000; Cooray & Sheth 2002; van den Bosch et al. 2013; Cacciato et al. 2013; Mead et al. 2015). The initial lens density profile is described by a Navarro-Frenk-White profile (Navarro et al. 1995). We used the mass-concentration relation given by Duffy et al. (2008) and allowed for a re-normalization factor, f_c (Viola et al. 2015).

A dominant source of systematic bias in stacked weak lensing analyses is a miscentering of the lenses, which can be due to an offset of the cluster halo with the visible distribution of galaxies (see, e.g., George et al. 2012) or the resolution of the cluster detection method (less than 0.1 Mpc h⁻¹ for AMICO; see Bellagamba et al. 2018). Following Johnston et al. (2007) as well as numerous subsequent works (e.g., Oguri et al. 2010; Viola et al. 2015; Dvornik et al. 2017; Bellagamba et al. 2019; Giocoli et al. 2021), we allowed a fraction, p_off, of clusters to be offset from the center of the galaxy distribution, effectively smoothing the central stacked ΔΣ profile with a characteristic radius, ℛ_off.

At large radii, typically beyond a few megaparsecs, the clustering of dark matter halos starts to dominate the signal. This “two-halo” term depends on the halo bias, b (Dvornik et al. 2017), and is modeled following Tinker et al. (2010). At small radii, the baryonic component of central galaxies can contribute to the signal, which is adequately described by a point mass, M_⋆, in the model (Viola et al. 2015; Dvornik et al. 2017).

In Table 1 we summarize these six free parameters for our halo model implementation, analogous to Dvornik et al. (2017).

In the AMICO cluster sample with 6925 lenses in 440 square degrees, many background sources are lensed by more than one cluster, contributing to the estimate of the ESD profile in various radial bins of different clusters. In the model fitting, we took the covariance between the ESD estimates into account, as described in Sect. 3.3.2 (Viola et al. 2015; Dvornik et al. 2017; Bellagamba et al. 2019).

3. Data and analysis

In this paper we use a lensing cluster catalog and a background source catalog from KiDS-450 (de Jong et al. 2017). KiDS is an optical wide-field imaging survey with OmegaCAM (Kuijken 2011) on the VLT Survey Telescope (VST; Capaccioli & Schipani 2011; de Jong et al. 2013). KiDS-450 consists of two patches, KiDS-N and KiDS-S (see Fig. 3), with 454 tiles of imaging data, for a total of 449.7 deg², in four optical filters, ugri. The survey was designed for lensing, ensuring a stable PSF, low seeing (< 0.96′, with an average of 0.66′ in r), and good photometric redshifts (photo-z; Hildebrandt et al. 2017).

Fig. 3. Overview of the KiDS-450 observations, with the KiDS-N (upper) and KiDS-S (lower) patches. The solid gray lines represent the planned KiDS survey area. Overplotted are the observed 1 square degree tiles, color coded with respect to their correspondence with the GAMA survey patches (G9 red, G12 yellow, G15 green, G23 blue, and GS purple; see Hildebrandt et al. 2017, for more details). The AMICO clusters analyzed in this work are represented by black dots.

The KiDS data were reduced with ASTRO-WISE (Valentijn et al. 2007; Verdoes Kleijn et al. 2012; Begeman et al. 2013; McFarland et al. 2013), as described in de Jong et al. (2015), Hildebrandt et al. (2017). Photometric redshifts, also termed z_B, were determined using a Bayesian photo-z estimation (BPZ; Benítez 2000; Coe et al. 2006) with PSF-matched photometry, as described in Hildebrandt et al. (2012, 2017), and Kuijken et al. (2015).

3.1. Lenses

We made use of the galaxy cluster catalog derived with AMICO (Bellagamba et al. 2011, 2018; Radovich et al. 2017), extracted from 440 tiles of KiDS-450 data and described in Maturi et al. (2019) and Bellagamba et al. (2019). For each cluster, the luminosity L₂₀₀ is defined¹² as the sum of r-band luminosities of bright candidate member galaxies, weighted by membership probability (see Maturi et al. 2019). We selected galaxies with k-corrected r-band magnitudes brighter than m^*(z_l)+1 within R₂₀₀(z_l), where z_l is the estimated cluster redshift and R₂₀₀ is derived from the adopted cluster model and is used in the construction of the cluster detection filter, as defined in Maturi et al. (2019). In this sense, L₂₀₀ is defined analogously to the apparent richness, λ^*, which is a sum of membership probabilities of galaxies with m < m^* + 1.5, within R₂₀₀.

We selected clusters in the range 0.1 ≤ z_l ≤ 0.6. We excluded clusters below z = 0.1 due to their unfavorable lensing geometry and above z = 0.6 due to the low density of background sources. For some clusters, no lens-source pairs were found, due to source selection criteria or masking. Our final selection comprises 6925 clusters, divided over the KiDS-450 survey area as shown in Fig. 3 and described in Table 2. The redshift distribution of these clusters is shown in Fig. 4, with a median redshift of z_l = 0.39.

Fig. 4. Redshift distribution of AMICO clusters (gray), with a median redshift of zl = 0.38 (dashed gray), and KiDS-450 background sources (purple), with a median redshift of z = 0.68 (dashed purple). In blue (DIR), we show the initial selection following Hildebrandt et al. (2017) and Dvornik et al. (2017). In red, we show the estimated redshift distribution of the gri color selection (COL), corresponding to the bottom panel of Fig. 5.

We divided the clusters into 13 bins of cluster L₂₀₀. The limits of these bins were chosen so that the signal-to-noise ratios of the ESD measurements were approximately the same in each bin. We give an overview of these bins, together with the estimated M₂₀₀, in Table 3.

3.2. Sources

We selected an initial sample of background sources using the same photometric redshift criteria as Hildebrandt et al. (2017), 0.1 < z_B ≤ 0.9, to reduce the outlier rate. We also applied the cut z_l + Δz < z_B, following Dvornik et al. (2017). Here, Δz = 0.2 is an offset between the redshift estimation, z_l, of the cluster by AMICO and the photometric redshift, z_B, of the source to sufficiently lessen the contamination of the background sample by cluster member galaxies (see also Appendix A.2).

Our selection of AMICO clusters is deeper than the lenses from the Galaxy And Mass Assembly (GAMA, Driver et al. 2011; Robotham et al. 2011) catalog used in Dvornik et al. (2017). As can be seen in Fig. 4, the redshift distributions of lenses and background sources significantly overlap, and the cut at Δz = 0.2 reduces the number density severely for clusters at higher redshift. Following Bellagamba et al. (2019), we also selected background sources using the color selection proposed by Oguri et al. (2012):

$$\begin{array}{c}\hfill g-r<0.3\vee r-i>1.3\vee r-i>g-r.\end{array}$$(27)

In Fig. 5 we show the photometric redshift distribution of this cut in the KiDS-450 catalog and compare it to the photometric and spectroscopic redshift distribution of the same cut in the spectroscopic redshift (spec-z) catalog used in Hildebrandt et al. (2017). Based on this analysis, we additionally required z_B ≥ 0.6 for this selection to reduce contamination by sources at low redshift and find that 98% of the galaxies in this color selection have z_spec > 0.6.

Fig. 5. Redshift distribution of background sources selected by color. Upper panel: the distribution of the photometric redshift, zB, of the sources in the KiDS-450 catalog that satisfy the color cut of Eq. (27), of which we select the sources with zB ≥ 0.6 (red) and discard those with zB < 0.6 (blue). Bottom two panels: the same selection applied to the spec-z catalog (Hildebrandt et al. 2017), plotted in terms of photometric redshift (middle) and spectroscopic redshift (lower). We find the contamination of sources with zB ≥ 0.6 and zspec < 0.6 is ∼2%.

3.2.1. Redshift distribution

To estimate the redshift distribution of background galaxies, we did not directly use the individual redshift probability distribution, p(z_s), per source galaxy. Instead, we applied a weighted direct calibration method (DIR), as motivated by Hildebrandt et al. (2017).

For each cluster, we used the spec-z catalog described in Hildebrandt et al. (2017) to select objects using the same selection criteria as described above. We then used the normalized spectroscopic redshift distribution, n(z_s), of this sample to calculate the comoving critical surface density analogous to Eq. (20):

$$\begin{array}{c}\hfill ⟨{\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{l}}^{-1}⟩=\frac{4\pi G}{c^2}D(z_{\mathrm{l}}){(1+z_{\mathrm{l}})}^2{\displaystyle \underset{z_{\mathrm{l}}+\mathrm{\Delta }z}{\overset{\infty }{\int }}\frac{D(z_{\mathrm{l}},z_{\mathrm{s}})}{D(z_{\mathrm{s}})}n(z_{\mathrm{s}})\mathrm{d}{z}_{\mathrm{s}}.}\end{array}$$(28)

The resulting redshift distribution for selected sources from the full KiDS-450 catalog is shown in Fig. 4.

3.2.2. Shape measurements

For shape measurements, the r-band data were reduced using the THELI pipeline, developed to meet the requirements for weak gravitational lensing analyses (Erben et al. 2005, 2009, 2013; Schirmer 2013). Galaxy shapes in the KiDS-450 catalog were then measured by lensfit (Miller et al. 2007, 2013; Kitching et al. 2008; Fenech Conti et al. 2017).

For each source, lensfit produces the ellipticity (ϵ₁,ϵ₂), an approximately inverse-variance weight w_s (see Miller et al. 2013), and a fitting quality parameter. We excluded sources with unreliable ellipticities from our source sample, using the same lensfit selection criteria as described in Hildebrandt et al. (2017).

In Paper I we compared the performance of estimators for ellipticity measurements in the CFHTLenS data with a subset of that catalog, selecting sources on the signal-to-noise ratio parameter ν_SN output by lensfit. We repeated that approach for a qualitative comparison here, using two subsets of the selected KiDS-450 sources. The first set selects sources with ν_SN ≥ 20, similar to Paper I, retaining ∼30% of the full background sample. The second set is a more stringent cut of the first set, additionally selecting objects with w_s ≥ 14.5, comprising ∼20% of the full sample.

3.2.3. Effective source density

The KiDS-450 catalog includes a filtering on general object detection and quality flags, for example, possibly blended sources or artifacts, as described by Kuijken et al. (2015) and Hildebrandt et al. (2017), and we discarded objects that lie in a mask. This removed approximately ∼12% of the sources. Our final selection comprises 14 124 197 sources, which translates to an effective number density of n_eff ≈ 8.23 arcmin⁻², as defined in Heymans et al. (2012b):

$$\begin{array}{c}\hfill n_{\mathrm{eff}}=\frac{1}{A}\frac{{\left({\sum }_i{w}_i\right)}^2}{{\sum }_i{w}_i^2},\end{array}$$(29)

with A the effective surface area, excluding masked regions.

3.3. Implementation

3.3.1. ESD estimation

Following Bellagamba et al. (2019), we measured the ESD in 14 logarithmic bins between 0.1 Mpc h⁻¹ and 3.16 Mpc h⁻¹. Not only does this make for an easy comparison of the results, but it has several other practical advantages.

We avoided radii smaller than the AMICO detection pixel size, which has a median size of 0.1 Mpc h⁻¹, to lessen the chance of a mismodeling the halo miscentering (Sect. 2.4). Here, the line of sight is also most contaminated by cluster members, which can lead to an overabundance by incorrectly including ellipticity measurements that carry no lensing signal, or by an obscuring and blending of background sources, which leads to an under-abundance of sources. While these effects may partially cancel out in the number counts, the effects on the ESD measurements do not cancel out, as the first leads to a diluted signal and the second to a very poor signal-to-noise ratio (see Appendix A.2 for an assessment of cluster member contamination).

At large radii, systematic additive biases can start to play a role (see, e.g., Dvornik et al. 2017, for this data set), which may differ for each KiDS survey patch (Fenech Conti et al. 2017; Hildebrandt et al. 2017). Another concern at larger separations is that the two-halo term becomes the dominant contribution to the ESD signal, which means we would need to properly constrain the halo bias, and we explain below how our approach does not fully take the clustering of dark matter halos into account.

The combination of background selection criteria from Dvornik et al. (2017) and Bellagamba et al. (2019) allows us to retain the three inner radial bins between 0.1 and 0.2 Mpc h⁻¹. We justify this inclusion in Appendix A.2, where we repeat the tests of Dvornik et al. (2017).

Each lens-source pair was then assigned a combined weight of

$$\begin{array}{c}\hfill w_{\mathrm{ls}}=w_{\mathrm{s}}{\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{l}}^{-2},\end{array}$$(30)

as motivated in Sect. 2.2. For LAD optimization, that is, the estimator that minimizes the L¹ norm (Eq. (24)), there exists no general analytic solution. The problem can, however, be formulated as a linear optimization, which can be solved iteratively (e.g., with simplex-based methods; Barrodale & Roberts 1973). In our weak lensing analyses, we find that convergence is robust.

To derive the covariance matrices for the ESD estimates using the mean and LAD in the same way, we can therefore also not employ the analytical prescription of Viola et al. (2015) used in earlier KiDS analyses (e.g., Sifón et al. 2015; van Uitert et al. 2016; Brouwer et al. 2016). Instead, we used a bootstrap approach.

Since the cluster bins of highest r-band luminosity, L₂₀₀, contain only a small number of clusters, covering only a small fraction of the KiDS-450 tiles, we cannot use the same bootstrap approach as Viola et al. (2015) and Dvornik et al. (2017) by bootstrapping 1 deg² tiles with replacement. Instead, we bootstrapped the source catalog, in accordance with Bellagamba et al. (2019).

This means that we are not sensitive to the clustering effect of dark matter halos, which justifies our choice of radial lens-source separation mentioned above. To assess the accuracy of these assumptions, we estimated the covariance matrix of the full 6925 cluster sample by bootstrapping the sources and by bootstrapping by KiDS-450 tiles in Appendix A.4. We conclude that our bootstrapping method yields a good estimate of the covariance matrix.

3.3.2. Halo model fitting

Having produced the LSQ and LAD shear profiles for the stacked clusters, we fit a halo model to the results. We used the fitting procedure described in Dvornik et al. (2017), producing the full posterior probabilities by a Bayesian inference technique, via a Monte Carlo Markov chain (MCMC) maximum likelihood approach. We assumed a Gaussian likelihood and made use of the full covariance between radial bins:

$$\begin{array}{c}\hfill \mathcal{L}\propto exp[-\frac{1}{2}{\mathit{R}}^{\mathrm{T}}{\mathit{C}}^{-1}\mathit{R}],\end{array}$$(31)

where the R are the residuals and C is the covariance matrix.

We used the emcee Python package (Foreman-Mackey et al. 2013) for the MCMC procedure, setting flat priors for all parameters. For the evaluation of the power spectrum and the halo mass function, we used the median redshift for each cluster luminosity bin.

4. Results

We present our results, starting with the derived ESD profiles obtained with the mean and LAD estimators, discussing potential biases and efficiency. Then, we show the results of the halo model fitting (i.e., M₂₀₀ for each luminosity bin) and conclude with the scaling relation between L₂₀₀ and M₂₀₀. We visualize the results for the case in which all 6925 clusters are stacked together, giving the numerical results of the 13 luminosity bins in Table 3.

4.1. ESD profiles

We calculated the ESD profiles using 10⁴ bootstraps with replacement. We estimated the ESD signal in the 14 radial bins, using both the mean and the LAD estimators, for each bootstrapped sample, preserving the bootstrap order of all 28 values throughout the whole process.

We find the estimator distribution to be almost perfectly normal, as expected from the central limit theorem. The correlation between the 14 bins of the full stack of clusters is shown in Fig. 6 and is given by

$$\begin{array}{c}\hfill {\rho }_{\mathit{ij}}\equiv \frac{{\mathrm{Cov}}_{\mathit{ij}}}{{\sigma }_i{\sigma }_j},\end{array}$$(32)

Fig. 6. Correlation matrix between the stacked ESD signals in different radial bins, using the full AMICO cluster catalog and 104 bootstraps. The upper-left triangle shows the correlation, ρLAD, between the LAD estimations. The lower-right triangle shows the correlation, ρmean, between the ESD estimates using a weighted mean. We note the general similarities in the two patterns.

where i and j denote the radial bin subscripts.

The upper-left part of the matrix shows the correlation between the LAD estimates of the radial bins, and the lower-right part the shows the mean results. Although the correlation between bins is very low, it is clear that the overall trends are the same for the two estimators.

The signal-to-noise ratio of the recovered ESD profile of the full stack, which is shown in Fig. 7, is high enough to allow us to notice the difference between the estimators, which indicates a small relative bias. The blue points show the LAD estimates, and the red points represent the mean estimates, with error bars in both cases defined as the square root of the diagonal elements of the covariance matrices (i.e., the classical standard deviation).

Fig. 7. Estimated ESD profile from the full AMICO cluster catalog, using the LAD estimator (blue) and a weighted mean (red). The error bars are the square roots of the diagonal values of the respective covariance matrices. The solid lines represent the best fitting halo model obtained by the MCMC fit. The shaded regions show the 68.3% confidence bands, estimated using the 15.9th and 84.1th percentiles of the MCMC realizations.

Tests for systematic effects, such as the cross signal, and a test for systematic additive noise around random points were already conducted by Dvornik et al. (2017) and Bellagamba et al. (2019). In Appendix A we repeat these tests for completeness since we use the KiDS-S field and an extended source selection with respect to Dvornik et al. (2017) and use a different source selection and three smaller radial bins with respect to Bellagamba et al. (2019). Our results show no residual systematic effects, in accordance with these papers.

4.2. Bias and efficiency

A possible bias is expected to depend on the strength of the underlying shear field since a zero lensing signal would imply no bias. In that case, the expected relevant distributions, tangential ellipticities or noise, are symmetric around zero ellipticity.

To quantify the difference between the ESD estimates, which we call the relative bias, ΔΣ_Mean − ΔΣ_LAD, we assumed¹³ to first order

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{Mean}}-\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{LAD}}=m·\mathrm{\Delta }\mathrm{\Sigma },\end{array}$$(33)

where we arbitrarily¹⁴ use Σ_LAD for Σ.

We used the full stack for its high signal-to-noise ratio, using the 15.9th and 84.1th percentiles of the differences in all bootstrap results to calculate uncertainties for each bin. We find m = −0.088 ± 0.020. In Fig. 8 we show this relative bias, plotting for visualization purposes

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{Mean}}-\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{LAD}}}{\mathrm{\Delta }\mathrm{\Sigma }}\end{array}$$(34)

Fig. 8. Difference between the recovered ESD signals in the radial bins, showing (ΔΣmean−ΔΣLAD)/ΔΣ. The solid white line represents the average difference, with the shaded region showing the formal 1σ error. As a possible bias in the recovered values is expected to increase with increasing shear, the differences are plotted versus the full ESD signal in each bin where we use ΔΣLAD, but we note that the small variations in the individual points, when plotting against ΔΣmean instead, give the same result, within statistical significance.

and a horizontal line at m = −0.88 to give a more intuitive impression of the relative error bars.

We reiterate that it is impossible to determine the absolute bias of each estimator as we did in Paper I, as we have no knowledge of the true ESD. However, the overall trend between the mean and LAD is similar in sign and order of magnitude, as we found for the CFHTLenS data in Paper I.

In Fig. 9 we show the derived relative efficiency (Eq. (26)) η = 1.047 ± 0.006, which is in accordance with the findings for CFHTLenS data in Paper I. The measured¹⁵ ellipticity distribution is expected to differ for shape measurements with a higher signal-to-noise ratio. For example, the cuspiness of the distribution shown in Fig. 2 can be smoothed out by noise convolutions. As the LAD estimation is more sensitive to the central peak, this will affect its precision.

Fig. 9. Relative efficiencies, $\eta \equiv {\sigma }_{\text{Mean}}^2/{\sigma }_{\text{LAD}}^2$, defined as the ratio of the diagonal elements of the covariance matrices. The relative efficiencies using the full KiDS-450 source catalog are shown in red. The purple and blue represent the higher S/N and higher lensfit weight selections, respectively. The solid lines represent the average η for each selection, with the shaded regions showing the 1σ errors.

In Paper I this was confirmed in the comparison of simulated data with and without noise, as well as in the results of the CFHTLenS sample with a stringent signal-to-noise selection compared with the full sample. As in Paper I, we compared the relative efficiencies for our selections with ν_SN ≥ 20 and with w_i ≥ 14.5, finding indeed a higher efficiency for less noisy shapes, namely η = 1.240 ± 0.010 and η = 1.386 ± 0.018, respectively.

4.3. Halo masses

We ran MCMC chains of 120 000 samples, using 120 walkers with 1000 steps each. The resulting chains were fully converged after the first 200 steps, so we discarded the first 24 000 samples.

We summarize the M₂₀₀ derived from the ESD estimation of the 13 luminosity bins in Table 3. Reduced χ², estimated between 0.730 and 2.528, are fairly consistent between derived results for mean and LAD.

For the full stack of clusters, we derived

$$\begin{array}{ccc}\hfill M_{200}=(0.{453}_{-0.030}^{+0.030})\times {10}^{14}{h}^{-1}{M}_{\odot },& {\chi }_{\nu }^2=1.25\hfill & \hfill (\mathrm{Mean})\end{array}$$(35a)

$$\begin{array}{ccc}\hfill M_{200}=(0.{487}_{-0.036}^{+0.033})\times {10}^{14}{h}^{-1}{M}_{\odot },& {\chi }_{\nu }^2=1.37\hfill & \hfill (\mathrm{LAD}).\end{array}$$(35b)

The confidence intervals are derived from the 15.9th and 84.1th percentiles of the posterior distributions. The best fitting ESD models are shown in Fig. 7. The 68.3% confidence bands overlap at some radii and are in tension at other radii. While the difference in ESD is significant, the 68.3% confidence intervals for M₂₀₀ just touch.

4.4. L₂₀₀ − M₂₀₀ scaling relation

We assumed a power-law relation between the derived halo masses and the median r-band luminosity of each cluster bin. We fit this relation in the form

$$\begin{array}{c}\hfill log\left(\frac{M_{200}}{M_{\mathrm{piv}}}\right)=a+blog\left(\frac{L_{200}}{L_{\mathrm{piv}}}\right),\end{array}$$(36)

with a the intercept and b the slope, where M_piv ≈ 10^14.1h⁻¹ M_⊙ and L_piv ≈ 10^11.8h⁻² L_⊙ are typical pivotal values of the halo mass and luminosity, derived from the fit itself. The fit was done in log basis as the derived posterior distributions of the halo mass are log-normal.

We did not take a redshift dependence into account, as Bellagamba et al. (2019) showed only a marginal and not very steep dependence of the halo mass on redshift.

We obtained the scaling relations

$$\begin{array}{cc}\hfill \frac{M_{200}}{{10}^{14.1}{h}^{-1}{M}_{\odot }}=(0.97\pm 0.06){\left(\frac{L_{200}}{{10}^{11.8}{h}^{-2}{L}_{\odot }}\right)}^{(1.24\pm 0.08)}& (\mathrm{Mean})\hfill \end{array}$$(37a)

$$\begin{array}{cc}\hfill \frac{M_{200}}{{10}^{14.1}{h}^{-1}{M}_{\odot }}=(1.03\pm 0.05){\left(\frac{L_{200}}{{10}^{11.8}{h}^{-2}{L}_{\odot }}\right)}^{(1.24\pm 0.08)}& (\mathrm{LAD})\hfill \end{array}$$(37b)

and plot the results in Fig. 10. As with the derived ESD profiles and halo masses for the full stack of clusters, the 68.3% confidence bands just touch at the pivot point

$$\begin{array}{c}\hfill (L_{200},M_{200})=({10}^{11.8}{h}^{-2}{L}_{\odot },{10}^{14.1}{h}^{-1}{M}_{\odot }),\end{array}$$(38)

Fig. 10. r-band luminosity-halo mass scaling relations, derived from the ESD profiles estimated using the weighted mean (red) and LAD (blue). At the pivot point, recognizable as the narrowest parts of the confidence bands, the relations just touch at the 68.3% confidence level.

recognizable as the narrowest parts of the confidence bands.

5. Summary and conclusions

We conducted a weak shear analysis of 6925 AMICO clusters in the KiDS-450 data We derived a tightly constrained scaling relation between r-band luminosity, L₂₀₀, and average lensing masses, M₂₀₀, in concordance with earlier results in the literature.

We investigated the impact of estimator choice for inferring the central moment of the cusped and skewed ellipticity distribution of background galaxies, finding a relative bias on the order of a few percent, as predicted in Paper I. We find that the constraints obtained via LAD regression are tighter than those obtained via LSQ regression, and they significantly improved as the signal-to-noise ratio of the shape measurements of the background galaxies increased. Complemented by simulations from Paper I, we give an alternative perspective on the problem of inferring the central shear value from the skewed distribution of background galaxy shapes, at the minor cost of increased, but still feasible, computation times for numerical iterativeregression.

5.1. L₂₀₀ − M₂₀₀ relation

Since the relative bias we found in both this research and in Paper I is approximately proportional to the ESD signal, it is expected that the LAD estimator will mainly have an effect on the intercept of the L₂₀₀ − M₂₀₀ scaling relation. This was confirmed by our results in Eq. (37).

The power-law index of the L₂₀₀ − M₂₀₀ scaling relation was constrained to 1.24 ± 0.08 (Eq. (37)), independent of estimator choice. This is in agreement with earlier work in the literature, such as Viola et al. (2015, and references therein), who cite 1.16 ± 0.13. This agreement is noteworthy since the AMICO clusters are derived from photometric redshifts, as opposed to the spectroscopically derived groups from GAMA (Driver et al. 2011; Robotham et al. 2011). The difference in confidence is explained by the increased number of lenses, which is not surprising given the large overlap in mass range: Viola et al. (2015) analyzed 1413 galaxy groups between z = 0.03 and z = 0.33, with r-band luminosity bin limits between 2.5 × 10⁹ and 5.0 × 10¹²h⁻² L_⊙, deriving halo masses between 1.4 × 10¹³ and 3.5 × 10¹⁴h⁻¹ M_⊙.

The choice of estimator produces a difference in intercept at just the 1σ level. Using a weighted mean to derive the ESD leads to an intercept of 0.97 ± 0.06 in the scaling relation, while LAD gives an intercept of 1.03 ± 0.05. This is to be expected as the relative bias seems roughly constant, when normalized by the ESD (see Fig. 8), and in good agreement with the bias of ∼5% found using simulations in Paper I.

As in Bellagamba et al. (2019), we estimate that systematic effects mostly affect the intercept. While the derived intercept is in agreement with the aforementioned papers, we note that the chosen definition¹⁶ of L₂₀₀, the difference in redshift range and definitions, and the completeness of group and cluster membershipcan account for possible differences on the same order of magnitude. This would not affect our conclusions on methodology, as the results from both estimators would be similarly affected.

In a further comparison with the scaling relation between richness and mass, cited in Bellagamba et al. (2019), we find similar significance in constraints on the slope. We define bins in luminosity instead of AMICO detection amplitude (Radovich et al. 2017) or richness, but since the cluster luminosity is tightly correlated with the richness, λ^*, this confirms our findings.

There are a few differences to consider. We chose not to account for a possible redshift dependence. This is motivated by Bellagamba et al. (2019) finding only a shallow dependence, which they point out may be driven mainly by the highest redshift bin. Our lensing analysis employs a slightly different background selection for an increased source density, combined with a different derivation of the associated redshift distribution, n(z). Another difference is the inclusion of radial bins at 0.1 ≤ R < 0.2 Mpc h⁻¹. This is expected to only have a minor effect as the contribution of the stellar mass is an order of magnitude lower than the halo term at these radii, while the contribution of miscentered halos only starts to become significant at larger radii (see also Rykoff et al. 2016; Oguri et al. 2018). In this sense, these findings are a confirmation of the robustness of the results across these papers.

5.2. Optimal estimators

Our results are in good agreement with Paper I, with a relative bias between the two estimators that shows the recovered lensing signal is higher with LAD, suggestive of a lower absolute bias. At the same time, LAD regression gives a small (albeit significant) gain in efficiency, giving a reduction in error bars of a few percent, and potentially up to 11%−18% for shape measurements of a higher signal-to-noise ratio. Least absolute deviation regression comes at the cost of a higher computation time, but at a step in the analysis process that does not dominate the total computational cost.

Both simulations (Paper I) and analyses on real data (this article) cite quantitative results of significance while at the same time showing similar trends between estimators on a qualitative level. We have conducted a cautious and thorough investigation but can never exclude the unknown: biases arising due to assumptions in the simulations of Paper I or uncorrected systematic effects in this research, or, most likely, both. However, given the range and realism in simulated distributions and the similarities in findings among those simulations, this research, and other work in the literature, we are confident that the recovered differences in results between the two estimators are real.

We note some differences between the two analyses. In Paper I we analyzed the regression of a sample of ellipticities with a single underlying value of the shear and, for each type of simulation, a single intrinsic ellipticity distribution. In this research, the situation is more complex. We studied the stacked signal around samples of lenses and of samples of background sources at a range of redshifts. This means also stacking noise that has been scaled by a range in lensing geometries, quantified by ${\mathrm{\Sigma }}_{\mathrm{cr}}^{-2}$. Furthermore, in each radial bin of each luminosity bin, we assume: (i) a constant lensing effect, which is in reality the stacked average of a range in L₂₀₀, and therefore a range in M₂₀₀, confounded by intrinsic scatter between these two quantities, and (ii) radial distance R from the lens, combined with a miscentering of halos.

All these effects tend to convolve the intrinsic galaxy shape distribution, which makes the level of agreement and significance between the two papers in fact remarkably robust. In conclusion:

• The combination of Paper I and this research shows that LAD regression is more naturally suited to the cusped intrinsic ellipticity distribution of background galaxies.

• Our simulations in Paper I showed a lower bias for LAD regression than for LSQ regression in the presence of noise in the background source shape measurements, while this research confirmed the same relative bias between the two estimators.

• Constraints obtained via LAD regression are comparable with or tighter than constraints obtained via LSQ regression.

An optimal estimator is, from a principled point of view, more objective and better suited than corrections to an approach, which is known to mismatch the sample distribution. More practically, LAD regression provides a robust consistency check for shear inference, which has been and still remains a major investment in the field of weak lensing. Keeping different perspectives, such as exploring these alternative statistical approaches, is fundamental for determining the way forward.

Acknowledgments

MSm acknowledges support from the Netherlands Organization for Scientific Research (NWO). AD acknowledges ERC Consolidator Grant (No. 770935) LM acknowledges the grants ASI-INAF n. 2018-23-HH.0 and PRIN-MIUR 2017 WSCC32 “Zooming into dark matter and proto-galaxies with massive lensing clusters”. MSe acknowledges financial contribution from contract ASI-INAF n.2017-14-H.0 and contract INAF mainstream project 1.05.01.86.10. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017 and 177.A-3018, and on data products produced by Target/OmegaCEN, INAF-OACN, INAF-OAPD and the KiDS production team, on behalf of the KiDS consortium. OmegaCEN and the KiDS production team acknowledge support by NOVA and NWO-M grants. Members of INAF-OAPD and INAF-OACN also acknowledge the support from the Department of Physics & Astronomy of the University of Padova, and of the Department of Physics of Univ. Federico II (Naples).

References

Appendix A: Tests for systematics

For completeness, we repeated some of the tests for systematics that were already carried out in Dvornik et al. (2017) and Bellagamba et al. (2019) because of our difference in sky coverage, background selection, and estimated redshift distribution compared with those two studies.

A.1. Photometric redshift

We used the same method as Dvornik et al. (2017) to determine the comoving critical density. There are two important differences that could affect the uncertainty in Σ_cr: We selected lenses at a significantly higher redshift, and we complemented our background source selection with the color selection described in Sect. 3.2.

We assessed the relative errors in Σ_cr by performing 10⁴ bootstraps of the spectroscopic catalog of Hildebrandt et al. (2017). We find the median error on Σ_cr to be ∼0.5%, as shown in Fig. A.1.

Fig. A.1. Relative errors in Σcr, estimated using 104 bootstraps of the spectroscopic catalog of Hildebrandt et al. (2017). The relative error on the ESD is negligible.

A.2. Contamination of the background sample by cluster galaxies

Dvornik et al. (2017) showed that an offset of Δz = 0.2 is enough to avoid a significant contamination of the background sources by unidentified GAMA group members. For lenses at a higher redshift, this contamination increases, while at the same time the density of available background sources decreases due to the observed depth of KiDS-450.

We used the same test as Dvornik et al. (2017) to assess the source density around AMICO clusters in order to determine the necessary Δz offset between the lens and the sources. We find that Δz = 0.2 is appropriate for our cluster selection (Fig. A.2).

Fig. A.2. Relative source densities around AMICO clusters as a function of radius, R, and the photometric redshift offset, zB ≥ zl + Δz, between the lens and the source. We note that some small under-densities around R = 0.2 h−1 Mpc may be due to the relative normalization.

A.3. Individual bin ESD profiles and cross signals

In Fig. A.3 we show the ESD profiles for the 13 cluster bins, including the cross signal, which is consistent with zero. We also show the derived halo model fits and their confidence intervals, comparing the fits using the full AMICO cluster catalog from Fig. 7.

Fig. A.3. Estimated ESD profiles and cross signal of the 13 cluster bins, using the LAD estimator (blue) and a weighted mean (red). The error bars are the square roots of the diagonal values of the respective covariance matrices. The solid lines represent the best fitting halo model obtained by the MCMC fit. The shaded regions show the 68.3% confidence bands, estimated using the 15.9th and 84.1th percentiles of the MCMC realizations. Individual bins are indicated by the range in normalized luminosity, 𝔏 ≡ L200 /(1010h−2 L⊙). The lower-right plot shows the ESD profile estimated from the full AMICO cluster catalog, also shown in Fig. 7. The average of the two best fitting halo models from Fig. 7 are shown in each panel as a dotted line for easy comparison.

A.4. Tile bootstrap

As described in Sect. 3.3, we could not use the same bootstrap approach as Dvornik et al. (2017), due to the sparsity of lenses in the highest lens luminosity bins. Since our bootstrap approach described in Sect. 3.3 does not account for cosmic variance and is not sensitive to the clustering effect of dark matter halos, we compare the covariances derived by the two bootstrap methods for the ESD of the whole lens selection in Fig. A.4 and find no significant differences or pattern beyond what is expected from statistical noise. Since we expect the contribution from cosmic variance to be even lower for subsets of lenses, we conclude that our bootstrap approach yields a good estimate of the covariance.

Fig. A.4. Left: Correlation derived from bootstrapping the signal by individual sources; same as Fig. 6, but with the color stretch adjusted to the middle plot. Middle: Correlation derived from bootstrapping the signal in 1 deg2 tiles. The upper-left corners show the correlations from LAD regression. The lower-right corners show the correlations from using the weighted mean. Right: Comparison of the errors obtained from bootstrapping sources (LAD: solid blue; mean: solid red) and bootstrapping 1 deg2 tiles (LAD: dashed blue; mean: dashed red).

Appendix B: Analysis of dependence on outer data points

In Fig. 10 it can be seen that the distribution of clusters in the two outermost luminosity bins is not symmetric. At the lower end, this is due to the selection criterion of λ^* in the AMICO catalog. At the higher end, we have only a few clusters.

We assessed the effect these two points have on the L₂₀₀ − M₂₀₀ scaling relation by repeating the fit without these bins. We find no difference within the statistical uncertainties, as given in Eq. B:

$$\begin{array}{cc}\hfill \mathrm{Mean}& \frac{M_{200}}{{10}^{14.1}{h}^{-1}{M}_{\odot }}=(0.98\pm 0.06){\left(\frac{L_{200}}{{10}^{11.8}{h}^{-2}{L}_{\odot }}\right)}^{(1.25\pm 0.10)},\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill \mathrm{LAD}& \frac{M_{200}}{{10}^{14.1}{h}^{-1}{M}_{\odot }}=(1.02\pm 0.06){\left(\frac{L_{200}}{{10}^{11.8}{h}^{-2}{L}_{\odot }}\right)}^{(1.23\pm 0.09)}.\hfill \end{array}$$(B.2)

All Tables

All Figures

	Fig. 1. Representation of a gravitational lens system, showing the displacement of a source at position S(x, y) to an image at position I(x′,y′), where we take the origin of the source plane to be collinear with the position of the lens, L, and the observer, O.
In the text

	Fig. 2. Estimated ellipticity and axis ratio distributions of sources in the KiDS-450 catalog. Top: two-dimensional histogram of ellipticities. Middle: histogram of the absolute ellipticity, |ϵ|. Bottom: histogram of the ellipse axis ratio, q.
In the text

Fig. 3. Overview of the KiDS-450 observations, with the KiDS-N (upper) and KiDS-S (lower) patches. The solid gray lines represent the planned KiDS survey area. Overplotted are the observed 1 square degree tiles, color coded with respect to their correspondence with the GAMA survey patches (G9 red, G12 yellow, G15 green, G23 blue, and GS purple; see Hildebrandt et al. 2017, for more details). The AMICO clusters analyzed in this work are represented by black dots.

In the text

Fig. 4. Redshift distribution of AMICO clusters (gray), with a median redshift of zl = 0.38 (dashed gray), and KiDS-450 background sources (purple), with a median redshift of z = 0.68 (dashed purple). In blue (DIR), we show the initial selection following Hildebrandt et al. (2017) and Dvornik et al. (2017). In red, we show the estimated redshift distribution of the gri color selection (COL), corresponding to the bottom panel of Fig. 5.

In the text

Fig. 5. Redshift distribution of background sources selected by color. Upper panel: the distribution of the photometric redshift, zB, of the sources in the KiDS-450 catalog that satisfy the color cut of Eq. (27), of which we select the sources with zB ≥ 0.6 (red) and discard those with zB < 0.6 (blue). Bottom two panels: the same selection applied to the spec-z catalog (Hildebrandt et al. 2017), plotted in terms of photometric redshift (middle) and spectroscopic redshift (lower). We find the contamination of sources with zB ≥ 0.6 and zspec < 0.6 is ∼2%.

In the text

	Fig. 6. Correlation matrix between the stacked ESD signals in different radial bins, using the full AMICO cluster catalog and 104 bootstraps. The upper-left triangle shows the correlation, ρLAD, between the LAD estimations. The lower-right triangle shows the correlation, ρmean, between the ESD estimates using a weighted mean. We note the general similarities in the two patterns.
In the text

Fig. 7. Estimated ESD profile from the full AMICO cluster catalog, using the LAD estimator (blue) and a weighted mean (red). The error bars are the square roots of the diagonal values of the respective covariance matrices. The solid lines represent the best fitting halo model obtained by the MCMC fit. The shaded regions show the 68.3% confidence bands, estimated using the 15.9th and 84.1th percentiles of the MCMC realizations.

In the text

Fig. 8. Difference between the recovered ESD signals in the radial bins, showing (ΔΣmean−ΔΣLAD)/ΔΣ. The solid white line represents the average difference, with the shaded region showing the formal 1σ error. As a possible bias in the recovered values is expected to increase with increasing shear, the differences are plotted versus the full ESD signal in each bin where we use ΔΣLAD, but we note that the small variations in the individual points, when plotting against ΔΣmean instead, give the same result, within statistical significance.

In the text

Fig. 9. Relative efficiencies, $\eta \equiv {\sigma }_{\text{Mean}}^2/{\sigma }_{\text{LAD}}^2$, defined as the ratio of the diagonal elements of the covariance matrices. The relative efficiencies using the full KiDS-450 source catalog are shown in red. The purple and blue represent the higher S/N and higher lensfit weight selections, respectively. The solid lines represent the average η for each selection, with the shaded regions showing the 1σ errors.

In the text

	Fig. 10. r-band luminosity-halo mass scaling relations, derived from the ESD profiles estimated using the weighted mean (red) and LAD (blue). At the pivot point, recognizable as the narrowest parts of the confidence bands, the relations just touch at the 68.3% confidence level.
In the text

	Fig. A.1. Relative errors in Σcr, estimated using 104 bootstraps of the spectroscopic catalog of Hildebrandt et al. (2017). The relative error on the ESD is negligible.
In the text

	Fig. A.2. Relative source densities around AMICO clusters as a function of radius, R, and the photometric redshift offset, zB ≥ zl + Δz, between the lens and the source. We note that some small under-densities around R = 0.2 h−1 Mpc may be due to the relative normalization.
In the text

Fig. A.3. Estimated ESD profiles and cross signal of the 13 cluster bins, using the LAD estimator (blue) and a weighted mean (red). The error bars are the square roots of the diagonal values of the respective covariance matrices. The solid lines represent the best fitting halo model obtained by the MCMC fit. The shaded regions show the 68.3% confidence bands, estimated using the 15.9th and 84.1th percentiles of the MCMC realizations. Individual bins are indicated by the range in normalized luminosity, 𝔏 ≡ L200 /(1010h−2 L⊙). The lower-right plot shows the ESD profile estimated from the full AMICO cluster catalog, also shown in Fig. 7. The average of the two best fitting halo models from Fig. 7 are shown in each panel as a dotted line for easy comparison.

In the text

Fig. A.4. Left: Correlation derived from bootstrapping the signal by individual sources; same as Fig. 6, but with the color stretch adjusted to the middle plot. Middle: Correlation derived from bootstrapping the signal in 1 deg2 tiles. The upper-left corners show the correlations from LAD regression. The lower-right corners show the correlations from using the weighted mean. Right: Comparison of the errors obtained from bootstrapping sources (LAD: solid blue; mean: solid red) and bootstrapping 1 deg2 tiles (LAD: dashed blue; mean: dashed red).

In the text