© ESO 2019

1. Introduction

According to the standard paradigm of galaxy physics, there must exist a strong correlation between the distributions of dark matter and galaxies (Mo et al. 2010). Semi-analytic models (SAMs) of galaxies are one method utilised to simulate the various complexities of galaxy physics by combining analytic prescriptions, where possible, with results from numerical simulations of the cosmological dark-matter density field (e.g. Springel et al. 2005 and Henriques et al. 2015, H15 hereafter). However, while the main physical mechanisms of galaxy evolution seem to be identified by observations, specific details concerning feedback and baryonic physics (of, in particular, star formation) are somewhat unclear. As a result, prescriptions that are encoded in SAMs are frequently phenomenological and employ free parameters calibrated to observables. Typical examples include the stellar-mass function of galaxies or, as for H15, the fraction of quiescent galaxies as function of redshift.

The galaxy-matter correlations as function of galaxy type, cosmic time (redshift), spatial scale, and environment encode useful statistical information about galaxy physics, and they can be used to test SAMs. An important tool for gathering this information is that of weak gravitational lensing, an effect which shears the shapes of distant galaxies (“sources”) through differential light deflection in the presence of a tidal gravitational field between the source and the observer (for a review, see Schneider et al. 2006). Importantly, the shear distortion, as described by the theory of general relativity, does not depend on the nature of the gravitating matter, making lensing an ideal probe for dark-matter physics.

Galaxy-galaxy lensing is one application in weak lensing that considers the correlation between positions of galaxies (lenses) and the shear signal of source galaxies in the background. This is a probe of the radial profile of the surface mass-density of matter around an average lens (e.g. Clampitt et al. 2017; Viola et al. 2015; Choi et al. 2012; Mandelbaum et al. 2006). Galaxy-galaxy-galaxy lensing is a recent extension of the second-order galaxy-galaxy lensing that probes the third-order correlations between the projected matter density and the galaxy number-density with two correlation functions (Simon et al. 2013, 2008; Watts et al. 2005; Schneider & Watts 2005). For the scope of this study, we only work with the lens-lens-shear correlation function and use this correlator synonymously with galaxy-galaxy-galaxy lensing. This kind of statistic is similar to galaxy-galaxy lensing in the sense that it measures the mean tangential shear around pairs of lenses or, after application of a lensing mass-reconstruction, the lensing convergence that is correlated with lens pairs (Simon et al. 2012). Because the lensing convergence is essentially the projected matter density on the sky, the map produced from a galaxy-galaxy-galaxy lensing visualises the typical matter environment of galaxy pairs in projection. With galaxy-galaxy-galaxy lensing being a connected three-point correlation function (by definition), the map shows the convergence in excess of the convergence around two individual galaxies. We therefore refer to this map as “excess mass map”. Since the introduction of galaxy-galaxy-galaxy lensing, alternative lensing measures of mass around average lens pairs have also been proposed and obtained from data, partly to probe the filamentary structure of the cosmic web (Epps & Hudson 2017; Clampitt et al. 2016; Johnston 2006).

As discussed in Saghiha et al. (2012) and recently shown by Saghiha et al. (2017, hereafter S17) galaxy-galaxy-galaxy lensing can test SAMs by comparing model predictions for the average matter density around galaxy pairs to measurements. The analysis in S17 is based on the CFHTLenS¹ measurements in Simon et al. (2013, herafter S13), and uses lens galaxies with stellar masses between ∼10⁹ − 10¹¹ M_⊙ and redshifts below z ≲ 0.6. S17 find the H15 model to be in good agreement with the CFHTLenS observations, while other models strongly disagree with the data. In addition, S17 argue that, for the same data, galaxy-galaxy-galaxy lensing has more discriminating power in this test than second-order galaxy-galaxy lensing. In particular, the strong dependence on galaxy morphology and galaxy colour makes galaxy-galaxy-galaxy lensing a powerful test for galaxy models.

We revisit the CFHTLenS data and the most promising galaxy model in S17 (i.e. the H15 model) for this paper, in an effort to gain more insight into the matter-galaxy relation on spatial scales of a few 100 h⁻¹ kpc. In particular, we create mass maps such that we are able to probe the matter environment of galaxy pairs on smaller physical scales than the related aperture statistics that are utilised in S17. These maps offer a better intuitive interpretation of the signal, by directly mapping out the average surface-matter density that is correlated with lens pairs for a fixed separation of the pair. Conversely the aperture statistics are useful for quantitative measurements because they are an average of the noisy galaxy-galaxy-galaxy-lensing correlation function for a broad range of separations, closely connected to the angular galaxy-matter bispectrum on the sky. We perform measurements of the excess mass around CFHTLenS galaxy pairs and, in a first study of this kind, compare these maps to the H15 model predictions. Moreover, for future studies, we introduce and investigate a promising new variant of the excess mass which we designate the “pair convergence”. In tests with synthetic data and from theoretical arguments, we find that it is less affected by nonphysical lens pairs that are merely close on the sky in projection (referred to as “chance pairs” in the following).

The structure of the paper is as follows. In Sect. 2, we introduce our notation as well as the definitions of all correlation functions relevant for second-order and third-order galaxy-galaxy lensing. Section 3 introduces the excess mass map and establishes its relation to the correlation function of galaxy-galaxy-galaxy lensing. Importantly, we discuss the effect of chance pairs on the excess map, and we introduce the pair convergence, a variation of the excess mass that is less affected by chance pairs. Our two investigated data sets, the simulated mocks and CFHTLenS data, are briefly described in Sect. 4. In Sect. 5 we outline the two estimators for the excess mass (or the pair convergence). We apply these estimators to CFHTLenS data and simulated H15 data in Sect. 6, verify our computer implementation of the CFHTLenS mapping code, and quantify its accuracy. Finally, in Sect. 7 we discuss our conclusions for the excess mass around CFHTLenS pairs and how they qualitatively compare to the H15 predictions.

2. Formalism

2.1. Lensing notation

Let δ_m(x) be the fractional density contrast $\delta {\rho }_{\mathrm{m}}/{\overline{\rho }}_{\mathrm{m}}$ in the matter-density field ρ_m(x) relative to the mean density ${\overline{\rho }}_{\mathrm{m}}$ at a comoving position x. We define positions x with respect to a fiducial light ray with the observer at the origin (Bartelmann & Schneider 2001). For this, let f_K( χ) θ at comoving distance χ be the transverse, comoving separation vector of a neighbour light ray from the fiducial ray, where f_K( χ) denotes the comoving angular diameter distance and θ the angular separation of the light ray on the sky. We assume a flat sky and denote angular positions by Cartesian vectors θ = θ₁ + iθ₂ in a complex notation with origin θ = 0 in the direction of the fiducial ray. For sources distributed along radial distance χ according to the probability density function (PDF) p_s( χ), the effective convergence at θ is a projection of δ_m(x) onto the sky:

$$\begin{array}{cc}\hfill \kappa (\mathit{\theta })& ={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi g(\chi ){\delta }_{\mathrm{m}}(f_K(\chi )\mathit{\theta },\chi );}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill g(\chi )& =\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{f_K(\chi )}{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{h}}}\mathrm{d}{\chi }^{\prime }{p}_{\mathrm{s}}({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })},}\hfill \end{array}$$(2)

where c is the vacuum speed of light, and H₀ is the Hubble constant; a( χ) is the scale factor with a( χ)=1 for χ = 0; and χ_h is the radius of the observable Universe (Schneider et al. 2006). The convergence is related to the (Cartesian) shear field γ_c(θ) up to a constant κ₀ through the convolution integral

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })-{\kappa }_0=\frac{1}{\pi }{\displaystyle \int {\mathrm{d}}^2{\theta }^{\prime }{\mathcal{D}}^{\ast }(\mathit{\theta }-{\mathit{\theta }}^{\prime }){\gamma }_{\mathrm{c}}({\mathit{\theta }}^{\prime })=:{\kappa }_{\mathrm{E}}(\mathit{\theta })+\mathrm{i}{\kappa }_{\mathrm{B}}(\mathit{\theta }),}\end{array}$$(3)

with the kernel

$$\begin{array}{c}\hfill {\mathcal{D}}^{\ast }(\mathit{\theta })=\frac{{\theta }_2^2-{\theta }_1^2+2\mathrm{i}{\theta }_1{\theta }_2}{{|\mathit{\theta }|}^4},\end{array}$$(4)

(Kaiser & Squires 1993). We call the real part κ_E(θ) of the convergence the E-mode of the convergence field and the imaginary part κ_B(θ) its B-mode. In the weak-lensing regime, lens-lens couplings are negligible (Hilbert et al. 2009) so that B-modes serve as indicator of systematic errors for our lensing analysis. In practical studies, we exploit that measurements of the image ellipticity of source galaxies can be converted into unbiased estimates of γ_c at the positions θ of the sources if |κ| ≪ 1 (e.g. Simon & Schneider 2017). In the following, we therefore denote by ϵ_i an unbiased estimator of γ_c(θ_i) at the position θ_i of a source galaxy.

We consider the tangential- and cross-shear components γ_t and γ_× of γ_c(θ) relative to an orientation angle φ, namely we define the φ-rotated shear,

$$\begin{array}{c}\hfill \gamma (\mathit{\theta };\phi )=-{\mathrm{e}}^{-2\mathrm{i}\phi }{\gamma }_{\mathrm{c}}(\mathit{\theta }),\end{array}$$(5)

and its decomposition

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{t}}(\mathit{\theta };\phi )+\mathrm{i}{\gamma }_{\times }(\mathit{\theta };\phi ):=\gamma (\mathit{\theta };\phi ).\end{array}$$(6)

Also, for mathematical convenience, we denote differences between two vectors θ_i and θ_j on the flat sky by

$$\begin{array}{c}\hfill {\mathit{\theta }}_{\mathit{ij}}:={\mathit{\theta }}_i-{\mathit{\theta }}_j={\theta }_{\mathit{ij}}{\mathrm{e}}^{\mathrm{i}{\phi }_{\mathit{ij}}},\end{array}$$(7)

with the polar coordinates θ_ij and φ_ij.

2.2. Galaxy clustering and galaxy-galaxy lensing

In galaxy-galaxy lensing, we correlate positions of lens galaxies with the tangential shear around these galaxies. For this, consider the number density N_g(θ) of lens galaxies on the sky and their density contrast

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{g}}(\mathit{\theta })=\frac{N_{\mathrm{g}}(\mathit{\theta })-{\overline{N}}_{\mathrm{g}}}{{\overline{N}}_{\mathrm{g}}},\end{array}$$(8)

relative to the mean number density N̅_g. Similar to κ(θ), the density contrast κ_g(θ) constitutes a projection along the line-of-sight θ. Specifically, let δ_g(x) be the relative fluctuations δn_g/n̅_g in the three dimensional number density of galaxies at a position relative to the fiducial light ray and p_d( χ) the PDF of galaxy distances χ inside the observed light cone. Then the density contrast of lenses on the sky is

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{g}}(\mathit{\theta })={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi p_{\mathrm{d}}(\chi ){\delta }_{\mathrm{g}}(f_K(\chi )\mathit{\theta },\chi )}\end{array}$$(9)

(e.g. Hoekstra et al. 2002).

Following Peebles (1980), we quantify the second-order angular clustering of lenses by the correlation of two density contrasts with separation ϑ,

$$\begin{array}{c}\hfill \omega (\vartheta )=\langle {\kappa }_{\mathrm{g}}(\mathit{\theta }){\kappa }_{\mathrm{g}}(\mathit{\theta }+\mathit{\vartheta })\rangle .\end{array}$$(10)

Owing to statistical isotropy and homogeneity of the random field δ_g, the correlation function ω(ϑ) depends only on the separation ϑ of the two points. Likewise, we assume isotropy and homogeneity for the matter-density fluctuations δ_m so that all following functions that correlate two or more points depend only on the mutual separation of points.

For a cross-correlation of lens positions with the projected matter density, we define the correlation between the number density of lenses and the mean tangential shear

$$\begin{array}{c}\hfill {\overline{\gamma }}_{\mathrm{t}}(\vartheta )=\frac{1}{{\overline{N}}_{\mathrm{g}}}\langle N_{\mathrm{g}}(\theta )\gamma (\mathit{\theta }+\mathit{\vartheta };\phi )\rangle =\langle {\kappa }_{\mathrm{g}}(\mathit{\theta })\gamma (\mathit{\theta }+\mathit{\vartheta };\phi )\rangle \end{array}$$(11)

with φ being the polar angle of ϑ = ϑ e^iφ. The imaginary part or cross-component ${\overline{\gamma }}_{\times }(\vartheta )$ of this correlator vanishes in the statistical average because of a parity invariance of the random fields (Schneider 2003). Physically, galaxy-galaxy lensing probes the (axis-symmetric) profile of the stacked (i.e. ensemble-averaged) convergence $\overline{\kappa }(\vartheta )=\langle {\kappa }_{\mathrm{g}}(\mathit{\theta })\kappa (\mathit{\theta }+\mathit{\vartheta })\rangle$ around lenses at separation ϑ:

$$\begin{array}{c}\hfill {\overline{\gamma }}_{\mathrm{t}}(\vartheta )=\left(\frac{2}{{\vartheta }^2}{\displaystyle {\int }_0^{\vartheta }\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }\overline{\kappa }({\vartheta }^{\prime })}\right)-\overline{\kappa }(\vartheta )\end{array}$$(12)

(Kaiser 1995).

2.3. Galaxy-galaxy-galaxy lensing

A correlation function similar to galaxy-galaxy lensing can be defined by measuring the average tangential shear at θ₃ around a pair of lenses at θ₁ and θ₂. This is introduced as one of two correlation functions of galaxy-galaxy-galaxy lensing in Schneider & Watts (2005) by

$$\begin{array}{cc}\hfill \mathcal{G}({\vartheta }_1,{\vartheta }_2,{\varphi }_3)=& \langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_1){\kappa }_{\mathrm{g}}({\mathit{\theta }}_2)\gamma \left({\mathit{\theta }}_3;\frac{{\phi }_1+{\phi }_2}{2}\right)\rangle \hfill \\ \hfill & =-{\mathrm{e}}^{-\mathrm{i}({\phi }_1+{\phi }_2)}\langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_1){\kappa }_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle ,\hfill \end{array}$$(13)

where ϕ₃ = φ₂ − φ₁ denotes the angle spanned by the two separation vectors θ₁ − θ₃ and θ₂ − θ₃ with polar angles φ₁ and φ₂ respectively. Figure 1 sketches the geometry. The tangential shear is thus defined relative to the line that bisects the angle between the two lenses.

Fig. 1. Geometry in the galaxy-galaxy-galaxy lensing correlation function 𝒢 for lenses at θ1 and θ2, and a source at θ3. The separations between source and lenses are ϑ1 and ϑ2. The separations of the lenses is θ12 = |θ1 − θ2|. This figure is copied from Schneider & Watts (2005).

From a mathematical point of view, 𝒢 is a correlator between the random field κ_g(θ) at θ₁ and θ₂, and the random field γ_c(θ) at θ₃. Since both fields vanish on average, i.e., ⟨κ_g(θ)⟩ = ⟨γ_c(θ)⟩ = 0 for all θ, the correlation function 𝒢 contains no additive contributions from correlators smaller than third order (the so-called unconnected terms) and thus vanishes for purely Gaussian random fields (e.g. Mo et al. 2010).

3. Mapping of mass correlated with lens pairs

3.1. Definitions of the excess mass

In the case of galaxy-galaxy lensing, the correlator ${\overline{\gamma }}_{\mathrm{t}}$ can be related to the stacked convergence around individual lenses. This is modified in the case of galaxy-galaxy-galaxy lensing, whereby the correlator 𝒢 can now be related to the stacked convergence around individual pairs of lenses (Simon et al. 2008, 2013). We show this by considering first the shear pattern around an average lens pair. Stacking the shear field around two lenses at θ₁ and θ₂ results in the average

$$\begin{array}{cc}& \frac{\langle N_{\mathrm{g}}({\mathit{\theta }}_1)N_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle }{\langle N_{\mathrm{g}}({\mathit{\theta }}_1)N_{\mathrm{g}}({\mathit{\theta }}_2)\rangle }=\frac{\langle [1+{\kappa }_{\mathrm{g}}({\mathit{\theta }}_1)][1+{\kappa }_{\mathrm{g}}({\mathit{\theta }}_2)]{\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle }{1+\omega ({\theta }_{12})}\hfill \\ \hfill & =\frac{\langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_1){\kappa }_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle +\langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_1){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle +\langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle }{1+\omega ({\theta }_{12})},\hfill \end{array}$$(14)

where we have used $\langle N_{\mathrm{g}}(\mathit{\theta })N_{\mathrm{g}}(\mathit{\theta }+\mathit{\vartheta })\rangle ={\overline{N}}_{\mathrm{g}}^2(1+\omega (\vartheta ))$, and define θ₁₂ as the separation of the lenses. In order to understand why the left-hand side of Eq. (14) equals the average shear around two lens galaxies, we consider the number-density field of lenses projected on a regular grid with a large number of micro-cells, each with a solid angle σ. We define the micro-cells such that they are sufficiently fine as to contain at most one lens, such that ∫_σdσ N_g(θ)=0 or 1. We then have a contribution to ⟨N_g(θ₁)N_g(θ₂)γ_c(θ₃)⟩ only for ∫_σdσ N_g(θ₁)=∫_σdσ N_g(θ₂)=1, while ⟨N_g(θ₁)N_g(θ₂)⟩ in the denominator is the probability to have a pair of galaxies at θ₁ and θ₂ at the same time (and is therefore a normalisation factor).

Now, through the definitions Eq. (5) and (11) for galaxy-galaxy lensing we additionally have $\langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_j){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle =-{\text{e}}^{2\text{i}{\phi }_{j3}}{\overline{\gamma }}_{\mathrm{t}}({\vartheta }_j)$ and can therefore cast Eq. (14) into

$$\begin{array}{cc}\hfill \langle {\kappa }_{\mathrm{g}}({\mathit{\theta }}_1){\kappa }_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle & =(1+\omega ({\theta }_{12}))\frac{\langle N_{\mathrm{g}}({\mathit{\theta }}_1)N_{\mathrm{g}}({\mathit{\theta }}_2){\gamma }_{\mathrm{c}}({\mathit{\theta }}_3)\rangle }{\langle N_{\mathrm{g}}({\mathit{\theta }}_1)N_{\mathrm{g}}({\mathit{\theta }}_2)\rangle }\hfill \\ \hfill & +{\mathrm{e}}^{2\mathrm{i}{\phi }_{13}}{\overline{\gamma }}_{\mathrm{t}}({\vartheta }_1)+{\mathrm{e}}^{2\mathrm{i}{\phi }_{23}}{\overline{\gamma }}_{\mathrm{t}}({\vartheta }_2).\hfill \end{array}$$(15)

This shows that 𝒢 is, apart from a phase factor, indeed related to the average shear around lens pairs given by the first term on the right-hand side – but rescaled with 1 + ω(θ₁₂) and in excess of the mean shear around individual lenses as given by the two terms that involve ${\overline{\gamma }}_{\mathrm{t}}$.

For the excess mass map, we consider a two-dimensional convergence map ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})$ that corresponds to the excess shear in Eq. (15) around lenses at given separation θ₁₂. We construct this map in a specific coordinate frame, for which ϑ is the relative separation from the map centre 𝒪. See the sketch in Fig. 2. The lenses are located at θ₁ = −θ₁₂ e_x/2 and θ₂ = +θ₁₂ e_x/2 with e_x being a unit vector in the x-direction. Applying the linear Kaiser–Squires transformation in Eq. (3) to Eq. (15) in this coordinate frame yields, up to a constant κ₀, a convergence map that can be expanded according to

$$\begin{array}{cc}& {\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})\hfill \\ \hfill & =(1+\omega ({\theta }_{12})){\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)+{\kappa }_0.\hfill \end{array}$$(16)

Fig. 2. Cartesian coordinate frame for a map of excess mass or the pair convergence. Lenses with separation θ12 are at the positions as indicated by the solid points. A location inside the map is identified by the two-dimensional vector ϑ.

Here the average convergence ${\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})$ around lens pairs corresponds to the shear stack ⟨N_gN_gγ_c⟩/⟨N_gN_g⟩ in Eq. (15), and the average convergence ${\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }|)=\langle N_{\mathrm{g}}(\mathit{\theta })\kappa (\mathit{\theta }+\mathit{\vartheta })\rangle /{\overline{N}}_{\mathrm{g}}$ around individual lenses corresponds to the average shear $-{\text{e}}^{2\text{i}\phi }{\overline{\gamma }}_{\mathrm{t}}(|\mathit{\vartheta }|)$, centred on the location of each lens. We emphasise that 1 + ω(θ₁₂) is a constant in this map, and that ${\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }|)$ is, by definition of 𝒢, the average convergence around all galaxies in the sample – including those that do not have a partner at separation θ₁₂.

Since the relation between κ(θ) and γ_c(θ) is only defined up to a constant κ₀, we cannot uniquely determine the excess mass map from the excess shear (see Eq. (3)). It is, however, reasonable to assume that ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})$, being the three-point correlation function ⟨κ_gκ_gκ⟩, quickly approaches zero for large ϑ which might be used to define κ₀. Alternatively, for the maps presented here, we fix κ₀ by asserting that ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}$ vanishes when averaged over the entire map area. We will neglect κ₀ in the following equations for convenience.

For what follows we also consider the pair convergence of lens pairs as a variant of the excess mass map, which is the straightforward difference signal

$$\begin{array}{cc}\hfill \overline{\mathrm{\Delta }\kappa }(\mathit{\vartheta };{\theta }_{12}):=& {\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)\hfill \\ \hfill & ={\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})-\omega ({\theta }_{12}){\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})\hfill \end{array}$$(17)

between the stacked convergence around lens pairs and the stacked convergence around two individual lenses. Since the excess mass map, originating from the connected correlation function 𝒢, is free of unconnected correlations (by definition), our interpretation is that the excess mass is the connected part of the pair convergence, and the extra term $-\omega ({\theta }_{12}){\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})$ is the unconnected part of the pair convergence.

While the excess mass exactly vanishes for Gaussian random fields, the pair convergence generally does not; although it is entirely determined by second-order correlations in this case. This can be seen from the definition Eq. (17) of $\overline{\mathrm{\Delta }\kappa }$ and Eq. (16) with ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}=0$, giving

$$\begin{array}{c}\hfill \overline{\mathrm{\Delta }\kappa }(\mathit{\vartheta };{\theta }_{12})=-\frac{\omega ({\theta }_{12})}{1+\omega ({\theta }_{12})}({\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)+{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)).\end{array}$$(18)

We visualise the excess mass as a two-dimensional map by plotting either ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}$, or $\overline{\mathrm{\Delta }\kappa }$ for the pair convergence, as function of ϑ for a fixed lens-lens separation θ₁₂ and orientation. The resulting maps have two known symmetries (Simon et al. 2008). First, there is a parity symmetry: correlation functions are unchanged under a reflection of shear and the lens density across an axis owing to the parity invariance of cosmological fields (Schneider 2003). As a consequence, quadrants in the maps are statistically consistent when mirrored across the line connecting two lenses. Second, another symmetry is present because we correlate density fluctuations κ_g at θ₁ and θ₂ from the same galaxy sample: a permutation of lens indices results in the same correlation function. These symmetries combine to produce, in the absence of noise, an exact reflection symmetry of the map with respect to both the x- and y-axes. We exploit this symmetry to enhance the S/N in the maps by averaging the quadrants inside each map.

Physically, the dimensionless quantities ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}$ and $\overline{\mathrm{\Delta }\kappa }$ are surface-mass densities in units of the (average) critical density ${\overline{\mathrm{\Sigma }}}_{\mathrm{crit}}:=1/\overline{{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}}$, defined by

$$\begin{array}{cc}& \overline{{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}}:=\hfill \\ \hfill & \frac{4\pi G_{\mathrm{N}}}{c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}{\chi }_{\mathrm{s}}\mathrm{d}{\chi }_{\mathrm{d}}{p}_{\mathrm{d}}({\chi }_{\mathrm{d}})p_{\mathrm{s}}({\chi }_{\mathrm{s}})\frac{a({\chi }_{\mathrm{d}})f_K({\chi }_{\mathrm{d}})f_K({\chi }_{\mathrm{s}}-{\chi }_{\mathrm{d}})}{f_K({\chi }_{\mathrm{s}})},}\hfill \end{array}$$(19)

where G_N is Newton’s gravitational constant. In our analysis, lenses have a typical distance of z_d ≈ 0.4 and sources of z_s ≈ 0.93 so that we estimate ${\overline{\mathrm{\Sigma }}}_{\mathrm{crit}}\approx 4.25\times {10}^{3}h{M}_{\odot }{\mathrm{pc}}^{-2}$ as the fiducial value for the critical surface mass density in our analysis below.

3.2. Impact of chance galaxy pairs

The following is a simple discussion exploring the impact of uncorrelated lens pairs on excess mass or pair convergence maps. In the construction of our maps, we select galaxy pairs within a sample by their angular separation θ₁₂ on the sky. Therefore, there will be pairs that are well separated in radial distance from each other, such that third-order correlations involving these are negligible for practical purposes. On the other hand, we will have pairs which have non-negligible third-order correlations with the lensing convergence as they are at similar distance. Naturally, however, making a distinction between what is negligible or otherwise is somewhat arbitrary. Nonetheless we could define a reasonable threshold for the correlation amplitude or a maximum radial separation of galaxies in a pair and use this to restrict our sample to those pairs with higher expected S/N. For the purpose of this simple discussion, however, we assume that we have a sample of lenses with a fraction 1 − p_tp of pairs for which a stack of convergence shall be free of any third-order correlations; thus giving the expectation value

$$\begin{array}{c}\hfill {\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})={\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)+{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)=:{\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{cp}}.\end{array}$$(20)

These are the “chance pairs”; sources which are pairs only in projection on-sky. The remaining fraction p_tp of (physically connected) pairs, to which we refer as “true pairs”, shall have a different yet unspecified stack ${\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})={\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}$ that depends in detail on the average surface mass-density around those pairs and the critical density ${\overline{\mathrm{\Sigma }}}_{\mathrm{crit}}$ at the distance of the pair. The stack around true pairs carries the interesting physical information so that, ideally, we would like to define an excess mass that is independent from chance pairs. This is neither true for an excess mass map nor for a pair convergence map.

We start by exploring the behaviour of pure samples of chance or true pairs. For pure samples of chance pairs (i.e. p_tp = 0) the excess mass and pair convergence vanish because mass cannot be correlated with two statistically independent lenses. Indeed, we will find ${\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})={\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{cp}}$, ω(θ₁₂)=0, and consequently ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}=\overline{\mathrm{\Delta }\kappa }=0$. On the other hand, if we have a pure sample of true pairs (p_tp = 1) with a clustering amplitude ω_tp (at separation θ₁₂), we will find

$$\begin{array}{c}\hfill {\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})=\\ \hfill & (1+{\omega }_{\mathrm{tp}}){\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)\hfill \\ \hfill & =:{\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}\hfill \end{array}$$(21)

for the excess mass. We obtain the equation for $\overline{\mathrm{\Delta }\kappa }(\mathit{\vartheta };{\theta }_{12})$ by setting ω_tp to zero in Eq. (21).

Usually we have a mixture of chance pairs and true pairs, and the impact of chance pairs on the excess mass is not immediately obvious. To explore this case, let us now assume that 0 <  p_tp <  1 and that the mean tangential shear ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$ is as unchanged from the pure case. In this mixture, the clustering amplitude is reduced to ω = p_tp ω_tp, and the convergence stack around all pairs,

$$\begin{array}{c}\hfill {\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})=p_{\mathrm{tp}}{\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){{|}_{\mathrm{tp}}+(1-p_{\mathrm{tp}}){\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12})|}_{\mathrm{cp}},\end{array}$$(22)

is the weighted average of the stacks in the pure samples. Then using Eqs. (20)–(22) in Eq. (16) results, after some algebra, in the excess mass of a mixed sample:

$$\begin{array}{c}\hfill {\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12})=\\ \hfill & p_{\mathrm{tp}}\frac{1+p_{\mathrm{tp}}{\omega }_{\mathrm{tp}}}{1+{\omega }_{\mathrm{tp}}}{\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}\hfill \\ \hfill & +p_{\mathrm{tp}}(1-p_{\mathrm{tp}})\frac{{\omega }_{\mathrm{tp}}^2}{1+{\omega }_{\mathrm{tp}}}({\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)+{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|)).\hfill \end{array}$$(23)

In conclusion, while the presence of chance pairs merely diminishes the overall amplitude of the true pair excess-mass ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}$ in the mixture, there is also a second term in the last line of Eq. (23) that changes the overall appearance of the map by adding an extra signal, mainly at the lens positions, that is proportional to the convergence stack around chance pairs.

This extra signal can be avoided in the pair convergence. Namely, by plugging Eqs. (20) and (22) into Eq. (17) we get

$$\begin{array}{c}\hfill \overline{\mathrm{\Delta }\kappa }(\mathit{\vartheta };{\theta }_{12})=p_{\mathrm{tp}}({\overline{\kappa }}_{\mathrm{pair}}(\mathit{\vartheta };{\theta }_{12}){|}_{\mathrm{tp}}-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_1|)-{\overline{\kappa }}_{\mathrm{ind}}(|\mathit{\vartheta }-{\mathit{\theta }}_2|))\end{array}$$(24)

for a mixture sample. The presence of chance pairs does at most change the overall amplitude in the pair convergence; each value in the map gives a lower limit to the pair convergence around true pairs inside the brackets of Eq. (24).

4. Data

4.1. The Canada-France-Hawaii Telescope Lensing Survey

The CFHTLenS is a multi-colour, wide-field lensing survey with measurements of galaxy photometry in the five bands u^*g′r′i′z′ (AB system), observed as part of the CFHT Legacy Survey Wide (Heymans et al. 2012). The survey covers 154 deg² of the sky, consisting of four contiguous fields: W1 (≈64 deg²), W2 (≈23 deg²), W3 (≈44 deg²), and W4 (≈23 deg²). The seeing is optimised for lensing measurements in the r′-band and lies typically between 0.​^″66 − 0.​^″82; the camera resolution is 0.​^″187 per CCD pixel. Each field is a mosaic of a set of Megacam pointings with 1 × 1 deg² field-of-view each. The data reduction uses the processing pipeline THELI (Erben et al. 2013), shear measurements of source galaxies are made using lensfit (Miller et al. 2013), and estimates for galaxy photometric redshifts and stellar masses are obtained using PSF-matched photometry and the computer code BPZ (Hildebrandt et al. 2012, Benitez 2000). The photometric estimator z_ph for the redshift is the maximum in the posterior redshift distribution of a galaxy. The final galaxy catalogue comprises these physical parameters of 7 × 10⁶ objects². For our analysis, we use only Megacam pointings that are flagged as “good” for lensing studies. This amounts to roughly 75 per cent of the pointings with an effective area of A ≈ 95 deg², which we determined by counting the number of unmasked pixels in the mask files of the 129 good pointings. For the selection of lens and source catalogues we follow the criteria in S13, where additionally several tests for systematic errors in the galaxy-galaxy-galaxy lensing measurements (for our samples) have already been performed. S13 also shows plots of the redshift distributions of galaxy samples.

Our source samples are galaxies with i′< 24.7, photometric redshifts of 0.65 ≤ z_ph <  1.2, and non-vanishing statistical weights w_i according to lensfit. This gives roughly 2.2 × 10⁶ sources inside the good pointings, and a mean w-weighted redshift of z̅_s ≈ 0.93. The lower limit in z_ph is chosen to reduce the overlap in redshift between lens and source samples in order to suppress systematic errors in the correlation functions. The effective number density of our sources is ${\overline{N}}_{\mathrm{s}}={(\sum _i{\mathit{w}}_i)}^2/(A\sum _i{\mathit{w}}_i^2)\approx 5.5{\mathrm{arcmin}}^{-2}$. To correct for the additive and multiplicative bias m_i in the estimators ϵ_i of the shear of the ith source, we follow the instructions in Miller et al. (2013). For this, the correction for m_i is easily included in our shear-related estimators by replacing ϵ_i ↦ ϵ_i (1 + m_i)⁻¹ and w_i ↦ w_i (1 + m_i). We note that we always have m_i >  −1 (see Appendix A in S13).

For the lens samples, we select galaxies with flux limit i′≤22.5 from two photo-z bins: 0.2 ≤ z_ph <  0.44 (“low-z”) and 0.44 ≤ z_ph <  0.6 (“high-z”). In addition, we select galaxies only from the stellar-mass range between 5 × 10⁹ ≤ M_sm/M_⊙ <  3.2 × 10¹¹ which combines all stellar-mass samples sm1 to sm6 in S13. This selection picks galaxies around the characteristic mass M_* ≈ 5 × 10¹⁰ M_⊙ of the stellar-mass function (Wright et al. 2017, and references therein). The estimates for stellar masses assume an initial-mass function according to Chabrier (2003) and have a typical RMS error of 0.3 dex (Velander et al. 2014). Counting only good pointings, we have in total 1.8 × 10⁵ galaxies in the low-z sample, yielding N̄_g ≈ 0.5 arcmin⁻², and 2.5 × 10⁵ galaxies in the high-z sample, N̅_g ≈ 0.7 arcmin⁻². The mean redshifts are z̅_d ≈ 0.35, 0.51 for low-z and high-z respectively. The RMS error of the redshift estimates is σ(z)≈0.04 (1 + z) with an outlier rate of roughly 3%.

The lower limit of allowed stellar masses in addition to the i′≤22.5 flux limit makes our lens sample approximately volume-limited inside the redshift intervals, as can be seen from Fig. A.1. The figure shows absolute rest-frame magnitudes and colours versus redshift for the CFHTLenS galaxies: the blue dots are galaxies with stellar-mass selection, and the orange dots are galaxies without stellar-mass selection. Only the few galaxies that are fainter than M_u ≳ −18.5 or M_g ≳ −19.5 are missing for z_ph ≳ 0.4. By comparing the distribution of blue and orange points in the colour plots we also see that the stellar-mass selection rejects galaxies bluer than M_g − M_r ≲ 0.3 or M_u − M_g ≲ 0.6 at all redshifts (left column in figure).

4.2. Synthetic lensing data and mock galaxies

Our mock lensing-data are generated by tracing the distortion of light bundles that traverse 64 independent light cones of the N-body Millennium Simulation (Springel et al. 2005; Hilbert et al. 2009). The Millennium Simulation is a purely dark-matter simulation with a spatial (comoving) resolution of 5 h⁻¹ kpc that is sampled by ∼10¹⁰ mass particles, populating a cubic region of comoving side length of 500 h⁻¹ Mpc. The fiducial cosmology of the Millenium Simulation has the following parameters: Ω_m = 0.25 = 1 − Ω_Λ for density parameters of matter and dark energy; Ω_b = 0.045 for the baryon density; σ₈ = 0.9 for the normalisation of the linear matter power spectrum at z = 0; a Hubble parameter of H₀ = h 100 km s⁻¹ Mpc⁻¹ with h = 0.73; and primordial fluctuations as in the Harrison-Zeldovich model. The fiducial cosmology in the Millennium Simulation is somewhat different compared to recent results by the Planck Collaboration XIII (2016). However, we expect that this recent update in cosmological parameters only mildly affects our model predictions, similar to the clustering properties of SAM galaxies in Henriques et al. (2017) where the authors rescale the Millennium Simulation results to the Planck cosmology for a comparison.

Regarding the synthetic lensing data, the simulated light-bundle distortions are an average over source distances with a probability distribution that is the same as that of CFHTLenS sources (see Fig. 5 in S13). Each of the 64 simulated light cones yields a 4 ° ×4° square with information on the theoretical lensing convergence and shear along the line-of-sight at 4096² pixel positions. The total area of the mock data is therefore 1024 deg². We mainly use the convergence grid for the estimator in Sect. 5.3. Only where we compare the convergence-stack method to the shear-stack method (Sect. 5.1), or where we quantify the model residuals in the maps, do we also generate synthetic shear catalogues. For these catalogues, we uniformly pick random source positions on the grid and assign each the shear value closest to that position as the observed source ellipticity. We therefore do not incorporate any shape noise.

For mock lens galaxies, we apply the H15 SAM-description, adjusted to the fiducial cosmology of the Millennium Simulation. We then follow the steps in S17 (specifically those listed in their Sect. 3.2), to obtain the mock lens samples for low-z and high-z, with a selection function that is consistent with that of galaxies in CFHTLenS, and including an emulation of statistical errors in both stellar masses and photometric redshifts. In particular, the mock samples have redshift distributions that are consistent with the CFHTLenS samples.

5. Estimators of excess mass

We present two estimators for both, the excess mass map (Eq. (16)) and the pair convergence map (Eq. (17)), which aim for different applications. In one version, we stack the shear field around pairs of galaxies and carry out a convergence reconstruction of the stack afterwards. This approach is suitable for observational weak-lensing data where only estimates of (reduced) shear are available for a set of discrete positions of source galaxies. In another version, we stack the convergence on contiguous grids directly, either for quick model predictions of the excess mass or to assess the accuracy of the shear-stacking method.

For practical implementations of both estimators, we note that frequently occurring phase factors e^2iφ_ij of separation vectors in our usual complex notation θ_ij = θ_ij e^iφ_ij are easily computed by ${\text{e}}^{2i{\phi }_{\text{i}j}}={\mathit{\theta }}_{ij}/{\mathit{\theta }}_{ij}^{\ast }$.

5.1. Shear stack

Let ${{\theta }}_i^{\text{d}}$ be positions of n_d lens galaxies on the sky, and (ϵ_i, w_i, ${{\theta }}_i^{\text{s}}$) the details of n_s source galaxies with ellipticities ϵ_i, statistical weights w_i, and positions ${{\theta }}_i^{\text{s}}$. From this, we compute the excess shear in Eq. (15) that is based on both the shear stack around galaxy pairs – the first term of the right-hand side of this equation – and the average shear around individual galaxies, which are the other terms on the right-hand side. Importantly, the excess shear is not identical to the average shear around galaxy pairs; it is usually only a small fraction of the latter.

The overall strategy for an estimator of the excess mass map is: (i) we estimate the excess shear by stacking source ellipticities around lens pairs in an appropriate reference frame and with weights 1 + ω(θ₁₂) (first term in Eq. (15)) and (ii) add the terms that involve ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$ to the stack (other terms); (iii) we apply the Kaiser–Squires inversion, Eq. (3), to obtain the excess mass map in Eq. (16); (iv) finally, we subtract a constant κ₀ from the map. The computation of the pair convergence is only slightly different as we explicitly set ω(θ₁₂)≡0 in this procedure. The following describes the details of the stacking and the convergence reconstruction. Therein we assume that estimates of ω(ϑ) and ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$ are already available; see the following section for estimators of those.

For stacking, we define a two-dimensional grid with N_p × N_p grid pixels (we choose N_p = 200); pixels shall have a square geometry. Each grid cell mn has the vector position p_mn = m + i n where (m, n) are its coordinates in a Cartesian stacking frame. Let (${{\theta }}_i^{\text{d}},{{\theta }}_j^{\text{d}}$) be the positions of a selected lens pair ij within the separation bin $\vartheta -\mathrm{\Delta }\vartheta /2\le |{\mathit{\theta }}_i^{\mathrm{d}}-{\mathit{\theta }}_j^{\mathrm{d}}|<\vartheta +\mathrm{\Delta }\vartheta /2$, and $({ϵ}_k,w_k,{\mathit{\theta }}_k^{\mathrm{s}})$ are the details of a source close to the lens pair. We map θ-coordinates to p-coordinates by a rotation α_ij and scaling |A_ij|, both encapsulated inside A_ij = |A_ij| e^iα_ij, and a translation B_ij,

$$\begin{array}{c}\hfill \mathit{p}=A_{\mathit{ij}}\mathit{\theta }+B_{\mathit{ij}}.\end{array}$$(25)

The complex-valued parameters A_ij and B_ij are determined by the mapping of the two lens positions to the fixed positions ${\mathit{p}}_1^{\mathrm{d}}$ and ${\mathit{p}}_2^{\mathrm{d}}$ in the stack,

$$\begin{array}{c}\hfill A_{\mathit{ij}}=\frac{{\mathit{p}}_2^{\mathrm{d}}-{\mathit{p}}_1^{\mathrm{d}}}{{\mathit{\theta }}_j^{\mathrm{d}}-{\mathit{\theta }}_i^{\mathrm{d}}};B_{\mathit{ij}}=\frac{{\mathit{\theta }}_j^{\mathrm{d}}{\mathit{p}}_1^{\mathrm{d}}-{\mathit{\theta }}_i^{\mathrm{d}}{\mathit{p}}_2^{\mathrm{d}}}{{\mathit{\theta }}_j^{\mathrm{d}}-{\mathit{\theta }}_i^{\mathrm{d}}},\end{array}$$(26)

(see Fig. 3 for an illustration). The positions of sources ${\mathit{\theta }}_k^{\mathrm{s}}$ are therefore ${\mathit{p}}_{\mathit{ijk}}^{\mathrm{s}}=A_{\mathit{ij}}{\mathit{\theta }}_k^{\mathrm{s}}+B_{\mathit{ij}}$ in the stacking frame. Additionally, the source ellipticities ϵ_k have to be rotated by ${ϵ}_k\mapsto {\text{e}}^{2\text{i}{\alpha }_{\mathit{ij}}}\times {ϵ}_k=A_{\mathit{ij}}/A_{\mathit{ij}}^{*}\times {ϵ}_k$ in the p-frame. Before mapping ϵ_k to the stacking frame we subtract off the average shear around each lens position in the θ-frame to obtain the excess shear. The complete stack for the excess shear at the grid pixel mn is then the weighted sum

$$\begin{array}{cc}& \mathrm{\Delta }{\gamma }_{\mathit{mn}}=\hfill \\ \hfill & \sum _{i,j,k=1}^{n_{\mathrm{d}},n_{\mathrm{s}}}\frac{{\mathrm{\Delta }}_{\mathit{ijk}}^{\mathit{mn}}{w}_k{\mathrm{e}}^{2\mathrm{i}{\alpha }_{\mathit{ij}}}}{W_{\mathit{mn}}}([1+\omega ({\theta }_{\mathit{ij}}^{\mathrm{dd}})]{ϵ}_k+{\mathrm{e}}^{2\mathrm{i}{\phi }_{\mathit{ik}}}{\overline{\gamma }}_{\mathrm{t}}({\theta }_{\mathit{ik}}^{\mathrm{ds}})+{\mathrm{e}}^{2\mathrm{i}{\phi }_{\mathit{jk}}}{\overline{\gamma }}_{\mathrm{t}}({\theta }_{\mathit{jk}}^{\mathrm{ds}})),\hfill \end{array}$$(27)

Fig. 3. Schematic of the stacking procedure. Positions θ on the sky are mapped to stacking frame positions p defined by the fixed positions of a selected lens pair (red points) in the stacking frame. The stacking-frame positions of the lenses are denoted by ${\mathit{p}}_1^{\mathrm{d}}$ and ${\mathit{p}}_2^{\mathrm{d}}$. The mapping involves a rotation so that source ellipticities in the stack have to be rotated.

where the total weight is

$$\begin{array}{c}\hfill W_{\mathit{mn}}=\sum _{i,j,k=1}^{n_{\mathrm{d}},n_{\mathrm{s}}}{\mathrm{\Delta }}_{\mathit{ijk}}^{\mathit{mn}}{w}_k,\end{array}$$(28)

and ${\mathrm{\Delta }}_{\mathit{ijk}}^{\mathit{mn}}=1$ flags if the source position ${\mathit{p}}_{\mathit{ijk}}^{\mathrm{s}}$ falls within the grid cell mn and ${\mathrm{\Delta }}_{\mathit{ijk}}^{\mathit{mn}}=0$ otherwise; by ${\theta }_{\mathit{ij}}^{\mathrm{dd}}$ we denote the separation between the lenses i and j, and by ${\mathit{\theta }}_{\mathit{ik}}^{\mathrm{ds}}$ the difference vector between the lens i and the source k; the angle φ_ik is the polar angle of ${\mathit{\theta }}_{\mathit{ik}}^{\mathrm{ds}}={\theta }_{\mathit{ik}}^{\mathrm{ds}}{\text{e}}^{\text{i}{\phi }_{\mathit{ik}}}$.

We then convert the stack Δγ_mn of excess shear into a map of the excess convergence. Owing to a sparse sampling of the shear stack by discrete source positions around lens pairs, an additional smoothing of this map is required. We apply this smoothing with a kernel K to Δγ_mn before the conversion to the excess mass map or the pair convergence map, namely by means of

$$\begin{array}{c}\hfill \mathrm{\Delta }{\gamma }_{\mathit{mn}}^{\mathrm{K}}=\frac{\sum _{n^{\prime },m^{\prime }=1}^{N_{\mathrm{p}},N_{\mathrm{p}}}{W}_{m^{\prime }{n}^{\prime }}K(m-m^{\prime },n-n^{\prime })\mathrm{\Delta }{\gamma }_{m^{\prime }{n}^{\prime }}}{\sum _{n^{\prime },m^{\prime }=1}^{N_{\mathrm{p}},N_{\mathrm{p}}}{W}_{m^{\prime }{n}^{\prime }}K(m-m^{\prime },n-n^{\prime })},\end{array}$$(29)

for the Gaussian kernel

$$\begin{array}{c}\hfill K(\delta m,\delta n)=exp(-\frac{1}{2}\frac{\delta m^2+\delta n^2}{{\sigma }_{\mathrm{rms}}^2}),\end{array}$$(30)

which has the kernel size σ_rms in units of our grid-pixel size. Using the weights W_mn for the smoothing ignores grid pixel with no shear information and gives more weight to pixels with a higher W_mn in the average of neighbouring pixels. We use a smoothing scale of σ_rms = 4 for our maps.

In the last step, we apply the algorithm by Kaiser & Squires (1993) to $\mathrm{\Delta }{\gamma }_{\mathit{mn}}^{\mathrm{K}}$ on the grid, employing Fast-Fourier Transformations, to obtain the a smoothed map ${\widehat{\mathrm{\Delta }\kappa }}_{\mathit{mn}}^{\mathrm{K}}=\mathrm{\Delta }{\kappa }_{\mathit{mn}}^{\mathrm{K}}+{\kappa }_0$ of the excess convergence with a constant offset κ₀. The real part in the excess convergence contains the E-mode of the signal, and the imaginary part is the B-mode. Applying the Kaiser–Squires technique on a finite field produces systematic errors which typically have the effect of increasing the signal towards the edges. We therefore remove 50 pixel from the outer edges of the grid in the final map. The inner cropped map has then the dimensions 100× 100 pixel². The constant offset κ₀ depends on the details of the implementation of the Kaiser–Squire algorithm and the number N_p of grid pixels. To have consistent maps in the following, we assert that the average excess-convergence over the cropped map has to vanish. We therefore subtract this average from the final map.

5.2. Galaxy-galaxy lensing and lens clustering

The second-order statistics ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$ and ω(ϑ) are estimated from the data by the following standard techniques. For the angular correlation function ω(ϑ) of the lens galaxies, we prepare a mock catalogue with n_r uniform random positions within the unmasked region of the survey. We then count the number DD(ϑ; Δϑ) of lens pairs within the separation bin [ϑ − Δϑ/2, ϑ + Δϑ/2), the number of random-galaxy pairs RD(ϑ; Δϑ), and the number of random-random pairs RR(ϑ; Δϑ). For the count rates, we consider all permutations of galaxy and mock positions, which means the total number of DD, RR, and DR for all separations equals n_d(n_d − 1), n_r(n_r − 1), and n_dn_r, respectively. According to Landy & Szalay (1993), we then estimate (for n_d, n_r ≫ 1)

$$\begin{array}{c}\hfill \omega (\vartheta ;\mathrm{\Delta }\vartheta )=\frac{n_{\mathrm{r}}^2}{n_{\mathrm{d}}^2}\frac{DD(\vartheta ;\mathrm{\Delta }\vartheta )}{RR(\vartheta ;\mathrm{\Delta }\vartheta )}-2\frac{n_{\mathrm{r}}}{n_{\mathrm{d}}}\frac{DR(\vartheta ;\mathrm{\Delta }\vartheta )}{RR(\vartheta ;\mathrm{\Delta }\vartheta )}+1,\end{array}$$(31)

for the angular clustering of lenses at separation ϑ.

To measure the mean tangential shear ${\overline{\gamma }}_{\mathrm{t}}$ and the cross shear ${\overline{\gamma }}_{\times }$ within the separation bin [ϑ − Δϑ/2, ϑ + Δϑ/2), we apply the estimator

$$\begin{array}{c}\hfill {\overline{\gamma }}_{\mathrm{t}}(\vartheta ;\mathrm{\Delta }\vartheta )+\mathrm{i}{\overline{\gamma }}_{\times }(\vartheta ;\mathrm{\Delta }\vartheta )=\frac{\sum _{d,s=1}^{n_{\mathrm{d}},n_{\mathrm{s}}}{\mathrm{\Delta }}_{\mathrm{ds}}(\vartheta ;\mathrm{\Delta }\vartheta )w_s(-{\mathrm{e}}^{-2\mathrm{i}{\phi }_{\mathrm{ds}}}{ϵ}_s)}{\sum _{d,s=1}^{n_{\mathrm{d}},n_{\mathrm{s}}}{\mathrm{\Delta }}_{\mathrm{ds}}(\vartheta ;\mathrm{\Delta }\vartheta )w_s},\end{array}$$(32)

where ${\text{e}}^{-2\text{i}{\phi }_{\mathrm{ds}}}={\mathit{\theta }}_{\mathrm{ds}}^{*}/{\mathit{\theta }}_{\mathrm{ds}}$ is the phase factor of ${\mathit{\theta }}_{\mathrm{ds}}={\mathit{\theta }}_s^{\mathrm{s}}-{\mathit{\theta }}_d^{\mathrm{d}}$, and Δ_ds(ϑ; Δϑ)=1 for ϑ − Δϑ/2 ≤ θ_ds <  ϑ + Δϑ/2 and Δ_ds(ϑ; Δϑ)=0 otherwise (e.g. Bartelmann & Schneider 2001).

Figure 4 is a summary of our measurements of ω(ϑ; Δϑ) and ${\overline{\gamma }}_{\mathrm{t}}(\vartheta ;\mathrm{\Delta }\vartheta )$ for both the CFHTLenS data and the H15 galaxy mocks. The measurements are sub-divided into the low-z and high-z redshift bins. We show the tight constraints from the mocks as 2σ regions, and the CFHTLenS measurements as large stars and square with 1σ error bars (obtained by jackknife re-sampling, see Sect. 5.5). The small data points in the bottom panel are the CFHTLenS cross-shear ${\overline{\gamma }}_{\times }(\vartheta ;\mathrm{\Delta }\vartheta )$ of the lens samples, which are plotted as absolute value in the logarithmic plot. The mean of all cross-shear data points is consistent with zero. Regarding a redshift dependence of the measurements, the low-z data points are somewhat higher than the high-z data points for both statistics, a trend which is also predicted by the H15 model. The model predictions for the tangential shear ${\overline{\gamma }}_{\mathrm{t}}$ are in very good agreement with the CFHTLenS data although the slope for the low-z profile is slightly shallower compared to H15 for ϑ ≲ 5 arcmin. A more detailed comparison of the galaxy-galaxy-lensing signal to H15 can be found in S17.

Fig. 4. Angular clustering (top panel) and mean tangential shear (bottom panel) of our different lens samples. The large data points show the CFHTLenS measurements for the low-z (stars) and high-z sample (squares) with 1σ error bars, and the coloured regions are 2σ predictions based on the H15 mocks. The grey regions are the predictions for the low-z samples, the magenta regions predict the amplitude of the high-z samples. The small data points with large errors bars at the bottom of the lower panel are the absolute values of the mean cross-shear for the CFHTLenS lenses.

The amplitude of the angular clustering ω(ϑ) in the model, on the other hand, is about 30% lower than observed. This might partly be explained by a distance distribution p_d( χ) of lenses that is actually narrower in CFHTLenS than the assumed distribution in the model. If so, galaxy-galaxy lensing ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$ would be little affected as long as the mean distance of lenses and sources is nevertheless sufficiently accurate. The clustering amplitude ω(ϑ), however, would be affected more strongly because it depends on $p_{\mathrm{d}}^2(\chi )$, although it is unlikely that a systematically broadened p_d( χ) in the model alone can fully explain the observed discrepancy. This would require the RMS variance of p_d( χ) in the model to be biased high by as much as 40%. We obtain this crude estimate by assuming a narrow top-hat shape for p_d( χ) with width Δr and centre r_c, as in Eq. (11) of Simon (2007). In this case, the amplitude of ω(ϑ) scales with (Δr)⁻¹.

More likely, the disagreement between the clustering of H15 galaxies and CFHTLenS galaxies reflects the current model uncertainties: in particular, Henriques et al. (2017) find a systematically too low amplitude of the galaxy clustering at z = 0.1 in the stellar mass range between $9\le {log}_{10}(M_{\mathrm{sm}}/h^{-2}{M}_{\odot })<10$ for the H15 model. They find an amplitude too low by ∼20 − 30% at a projected separation between the galaxies of a few 100 h⁻¹ kpc (see their Fig. 13). This comparison to galaxies from the Sloan Digital Sky Survey, albeit at z = 0.1 instead of 0.2 ≲ z ≲ 0.6, is consistent with our finding. The authors argue that the systematic error in the galaxy clustering is related to the treatment of supernova feedback and the gas reincorporation time in the model, affecting the clustering and prevalence of low-mass galaxies. We also refer to the recent work by Cohn (2017) for a thorough discussion on the impact of SAM and simulation parameters on several observational properties of galaxies. As to measurements of galaxy-galaxy-galaxy lensing in this study, we point out that the estimator for the excess mass map uses 1 + ω(ϑ) rather than ω(ϑ) at ϑ ≈ 1 arcmin for which the discrepancy is around 7%, and for the pair convergence map we do not use the angular clustering of lenses at all.

5.3. Convergence stack

In a variant of the previous estimator for maps of the excess mass, we stack the excess convergence $\mathrm{\Delta }{\overline{\kappa }}_{\mathrm{emm}}$ or $\mathrm{\Delta }\overline{\kappa }$ around lens pairs in simulated data directly. Let ${\mathit{\theta }}_i^{\mathrm{d}}$ be the positions of n_d lens galaxies and κ(θ) a simulated grid of convergence values. Similar to Sect. 5.1, we use a shear-free affine transformation to map θ positions around a given pair ij of lenses to the stacking frame with coordinates p. The estimation process proceeds in four steps: (i) we stack the convergence around a set of selected lens pairs from the separation bin $\vartheta -\mathrm{\Delta }\vartheta /2\le |{\mathit{\theta }}_i^{\mathrm{d}}-{\mathit{\theta }}_j^{\mathrm{d}}|<\vartheta +\mathrm{\Delta }\vartheta /2$ to obtain ${\overline{\kappa }}_{\mathrm{pair}}$; (ii) we stack the convergence around individual lenses to obtain ${\overline{\kappa }}_{\mathrm{ind},1}+{\overline{\kappa }}_{\mathrm{ind},2}$ in the stacking frame; (iii) we estimate 1 + ω averaged for the distribution of lens-lens separations in the sample of selected lens pairs; and (iv) we combine the steps (i) to (iii) to compute the excess mass map ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}=(1+\omega ){\overline{\kappa }}_{\mathrm{pair}}-{\overline{\kappa }}_{\mathrm{ind},1}-{\overline{\kappa }}_{\mathrm{ind},2}$ with an estimate of ω, which is averaged over the distribution of lens-lens separations in the stack. As before, setting ω ≡ 0 in (iii) yields the pair convergence map.

For the steps (i) and (ii), we use square grids with N_p × N_p pixels and coordinates p_mn = m + i n. The lens positions are defined to be at the fixed location ${\mathit{p}}_1^{\mathrm{d}}$ and ${\mathit{p}}_2^{\mathrm{d}}$. Using the definitions Eq. (26) for the parameters of the mapping between the p-frame and the θ-frame for a given lens pair ij, we obtain

$$\begin{array}{c}\hfill {\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}}=A_{\mathit{ij}}^{-1}({\mathit{p}}_{\mathit{mn}}-B_{\mathit{ij}})\end{array}$$(33)

for the position θ in the convergence grid that corresponds to the stack position p_mn. We then compute the stack ${\overline{\kappa }}_{\mathrm{pair}}$ for the grid pixel mn by the average

$$\begin{array}{c}\hfill {\overline{\kappa }}_{\mathrm{pair}}^{\mathit{mn}}=\frac{{\sum }_{i,j=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}\kappa ({\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}})}{{\sum }_{i,j=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}},\end{array}$$(34)

where ${\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}=1$ indicates if ${\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}}$ is inside the κ-grid or ${\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}=0$ otherwise. In case of ${\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}=1$, we choose the grid point in κ(θ) that is closest to ${\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}}$.

To obtain a map of the average convergence around individual lenses in step (ii), we have to factor in that the convergence around pairs in the stack is differently scaled for any new lens pair in the stack. Therefore, to obtain ${\overline{\kappa }}_{\mathrm{ind},1}+{\overline{\kappa }}_{\mathrm{ind},2}$ for a distribution of scale parameters |A_ij| we apply the following technique. For each lens galaxy at ${\mathit{\theta }}_i^{\mathrm{d}}$, we randomly pick a position ${\mathit{\theta }}_i^{\mathrm{rnd}}={\mathit{\theta }}_i^{\mathrm{d}}+\delta {\theta }_i{\text{e}}^{\text{i}{\varphi }_i}$ for an “imaginary” uncorrelated lens, where ϕ_i defines a uniformly random orientation ϕ_i ∈ [0, 2π) and δθ_i is a random separation from the ith lens. For δθ_i, we randomly pick a separation $|{\mathit{\theta }}_k^{\mathrm{d}}-{\mathit{\theta }}_l^{\mathrm{d}}|$ from a list of all separations of selected lens pairs in the data. This assures that (ii) consistently uses the same distribution of rescalings |A_ij| that are applied in step (i). Since the positions of imaginary lenses are uncorrelated with both κ(θ) and the lens positions, a stack around lens-imaginary pairs yields (on average) ${\overline{\kappa }}_{\mathrm{ind},1}$ if the lens is mapped to ${\mathit{p}}_1^{\mathrm{d}}$ and the imaginary lens mapped to ${\mathit{p}}_2^{\mathrm{d}}$; we obtain ${\overline{\kappa }}_{\mathrm{ind},2}$ if we swap the lens positions in the stack. Therefore, by adding the maps for both cases we obtain the stack ${\overline{\kappa }}_{\mathrm{rnd}}={\overline{\kappa }}_{\mathrm{ind},1}+{\overline{\kappa }}_{\mathrm{ind},2}$, which justifies the estimator

$$\begin{array}{c}\hfill {\overline{\kappa }}_{\mathrm{rnd}}^{\mathit{mn}}=\frac{{\sum }_{i=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}\kappa ({\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}})}{{\sum }_{i=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}}+\frac{{\sum }_{i=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ji}}^{\mathit{mn}}\kappa ({\mathit{\theta }}_{\mathit{ji}}^{\mathit{mn}})}{{\sum }_{i=1}^{n_{\mathrm{d}}}{\mathrm{\Delta }}_{\mathit{ji}}^{\mathit{mn}}},\end{array}$$(35)

where ${\mathrm{\Delta }}_{\mathit{ij}}^{\mathit{mn}}$ and ${\mathit{\theta }}_{\mathit{ij}}^{\mathit{mn}}$ are defined as before with the exception that we use ${\mathit{\theta }}_j^{\mathrm{d}}\equiv {\mathit{\theta }}_i^{\mathrm{rnd}}$ for the position of the second lens in the lens pair ij.

Finally, we combine the information from the previous steps to compute the excess mass around the lens pairs by

$$\begin{array}{c}\hfill {\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}^{\mathit{mn}}=[1+\omega (\vartheta ;\mathrm{\Delta }\vartheta )]{\overline{\kappa }}_{\mathrm{pair}}^{\mathit{mn}}-{\overline{\kappa }}_{\mathrm{rnd}}^{\mathit{mn}}.\end{array}$$(36)

For a consistent comparison with maps obtained by shear stacking, Sect. 5.1, we smooth the excess mass map with the same kernel as in Eq. (29),

$$\begin{array}{c}\hfill {\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},\mathrm{K}}^{\mathit{mn}}=\frac{\sum _{n^{\prime },m^{\prime }=1}^{N_{\mathrm{p}},N_{\mathrm{p}}}K(m-m^{\prime },n-n^{\prime }){\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}^{m^{\prime }{n}^{\prime }}}{\sum _{n^{\prime },m^{\prime }=1}^{N_{\mathrm{p}},N_{\mathrm{p}}}K(m-m^{\prime },n-n^{\prime })},\end{array}$$(37)

we apply the same cropping, and we subtract the average of ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},\mathrm{K}}^{\mathit{mn}}$ from the cropped map.

5.4. Combining measurements

Our data consist of i = 1…n_f separate fields: four fields W1 to W4 for the CFHTLenS data and 64 fields for the synthetic data. Separated fields means here that we ignore contributions from galaxy pairs or triples where not all galaxies are inside the same field. For a combined measurement, we apply the estimators described in the previous sections for each field individually and then average them as described in the following. Our strategy for performing measurements of galaxy-galaxy-galaxy lensing with CFHTLenS data is an improvement in comparison to Simon et al. (2013). In that work, measurements in n_f = 129 individual pointings were performed and combined afterwards for a final result. Here, using the continuous fields W1–4, each consisting of many adjacent pointings, allows us to include also galaxy tuples with galaxies from different pointings. We find that this new strategy can enhance the overall S/N in the CFHTLenS maps moderately by 10–30 per cent, depending on the lens samples and their redshift binning.

For a combined estimate of the mean tangential shear ${\overline{\gamma }}_{\mathrm{t}}(\vartheta )$, we imagine the application of Eq. (32) to a merged catalogue of all fields where positions between any pair of galaxies from separate fields are larger than the considered range of ϑ. For this merged catalogue, let $n_{\mathrm{d}}^i$ and $n_{\mathrm{s}}^i$ be the number of lenses and sources inside field i, and ${\overline{\gamma }}_{\mathrm{t}}^i+\text{i}{\overline{\gamma }}_{\times }^i$ the estimator Eq. (32) for galaxies in field i only. We then split the sums over lens-source pairs in Eq. (32) for the merged catalogue into additive contributions from each field to obtain

$$\begin{array}{c}\hfill {\overline{\gamma }}_{\mathrm{t}}(\vartheta ;\mathrm{\Delta }\vartheta )+\mathrm{i}{\overline{\gamma }}_{\times }(\vartheta ;\mathrm{\Delta }\vartheta )=\frac{{\sum }_{i=1}^{n_{\mathrm{f}}}{W}^i({\overline{\gamma }}_{\mathrm{t}}^i(\vartheta ;\mathrm{\Delta }\vartheta )+\mathrm{i}{\overline{\gamma }}_{\times }^i(\vartheta ;\mathrm{\Delta }\vartheta ))}{{\sum }_{i=1}^{n_{\mathrm{f}}}{W}^i},\end{array}$$(38)

with the field weights

$$\begin{array}{c}\hfill W^i={\displaystyle \sum _{s,d=1}^{n_{\mathrm{d}}^i,n_{\mathrm{s}}^i}{\mathrm{\Delta }}_{\mathrm{ds}}^i(\vartheta ;\mathrm{\Delta }\vartheta )w_s^i,}\end{array}$$(39)

for the combined estimate. Here $w_s^i$ are the statistical weights of sources s in field i, and the flag ${\mathrm{\Delta }}_{\mathrm{ds}}^i(\vartheta ;\mathrm{\Delta }\vartheta )$ applies to positions in field i only.

For the combined estimate of 1 + ω, we have to determine the count rates DD, DR, and RR in the merged catalogue. To this end, we cannot simply add the count rates of all individual fields. Instead we have to pay attention to the field variations of numbers $n_{\mathrm{d}}^i$ of observed galaxies and of unclustered random galaxies in the merged survey. For example, the number of random positions inside field i should depend on the effective area of the field. To quantify this, let p_i be the probability that a random position of unclustered galaxies is inside field i, and n_r is the total number of random points for all fields; the distribution of mocks in field i complies with the selection function in field i. Only if all fields have the same effective area (or selection function in general), we will find p_i = p_j = 1/n_f for i ≠ j. This, for example, is exactly the case for the 64 fields in our simulated data, and approximately for the 129 CFHTLenS pointings that make up the fields W1–4. For a measurement of the angular clustering of lenses in CFHTLenS, we combine the count rates inside the individual pointings, hence here n_f = 129; for all other correlation functions we perform measurements inside the large fields W1-4. Now, counting the total number of random-random pairs in the merged catalogue we find

$$\begin{array}{cc}& \frac{RR(\vartheta ;\mathrm{\Delta }\vartheta )}{n_{\mathrm{r}}^2}=\hfill \\ \hfill & \frac{1}{n_{\mathrm{r}}^2}{\displaystyle \sum _{i=1}^{n_{\mathrm{f}}}R{R}^i(\vartheta ;\mathrm{\Delta }\vartheta )=\sum _{i=1}^{n_{\mathrm{f}}}{p}_i^2\frac{R{R}^i(\vartheta ;\mathrm{\Delta }\vartheta )}{p_i^2{n}_{\mathrm{r}}^2}=:\sum _{i=1}^{n_{\mathrm{f}}}{p}_i^{2}r{r}^i(\vartheta ;\mathrm{\Delta }\vartheta ),}\hfill \end{array}$$(40)

where $r{r}^i:=R{R}^i/(p_i^2{n}_{\mathrm{r}}^2)$ is the count rate RRⁱ in field i normalised with the total number (p_in_r)² of random pairs in this field. Conveniently, the value of rrⁱ does not depend, on average, on the absolute number of mock positions that we actually in the individual measurement of field i. Therefore, for the combined result of normalised counts $rr=RR/n_{\mathrm{r}}^2$, we take the average of all individual rrⁱ weighted with $p_i^2$. Similarly, we obtain for the normalised count rate dr = DR/(n_dn_r) of DR pairs in the merged catalogue

$$\begin{array}{c}\hfill dr(\vartheta ;\mathrm{\Delta }\vartheta ):={\displaystyle \sum _{i=1}^{n_{\mathrm{f}}}\frac{p_i{n}_{\mathrm{d}}^i}{n_{\mathrm{d}}}\frac{D{R}^i(\vartheta ;\mathrm{\Delta }\vartheta )}{n_{\mathrm{d}}^i{p}_i{n}_{\mathrm{r}}}=\sum _{i=1}^{n_{\mathrm{f}}}{p}_i{f}_{i}d{r}^i(\vartheta ;\mathrm{\Delta }\vartheta ),}\end{array}$$(41)

where $f_i=n_{\mathrm{d}}^i/n_{\mathrm{d}}$ is the fraction of galaxies in field i, and drⁱ refers to the normalised rate in field i which, as before, does not depend on the number of mock positions used. Consequently, we compute the combined dr by taking a weighted sum of individual drⁱ. The normalised count $dd=DD/n_{\mathrm{d}}^2$ of DD pairs is

$$\begin{array}{c}\hfill dd(\vartheta ;\mathrm{\Delta }\vartheta ):={\displaystyle \sum _{i=1}^{n_{\mathrm{f}}}{f}_i^{2}d{d}^i(\vartheta ;\mathrm{\Delta }\vartheta ).}\end{array}$$(42)

In summary, we compute the count rates (ddⁱ, drⁱ, rrⁱ) for each separate field i and then perform the previous weighted sums for the combined rates (dd, dr, rr). By means of the estimator (31), we then get

$$\begin{array}{c}\hfill \omega (\vartheta ;\mathrm{\Delta }\vartheta )=\frac{dd(\vartheta ;\mathrm{\Delta }\vartheta )}{rr(\vartheta ;\mathrm{\Delta }\vartheta )}-2\frac{dr(\vartheta ;\mathrm{\Delta }\vartheta )}{rr(\vartheta ;\mathrm{\Delta }\vartheta )}+1,\end{array}$$(43)

for the combined clustering amplitude of all fields.

With regard to combining measurements of the shear stacks from a set of n_f separate fields we do the following. We apply the estimator Eq. (27) with identical grid parameters to each field i for $\mathrm{\Delta }{\gamma }_{\mathit{mn}}^i$ and the statistical weights $W_{\mathit{mn}}^i$. Therein, we constantly use the clustering amplitude 1 + ω as estimated once from the merged catalogue and the measurement of the tangential shear ${\overline{\gamma }}_{\mathrm{t}}^i$ in field i. This is consistent with the technique in Simon et al. (2008) but a slight variation compared to Simon et al. (2013) where 1 + ω is measured for each CFHTLenS pointing individually. We then combine the measurements of the excess shear inside the fields into

$$\begin{array}{c}\hfill \mathrm{\Delta }{\gamma }_{\mathit{mn}}=\frac{{\sum }_{i=1}^{n_{\mathrm{f}}}{W}_{\mathit{mn}}^i\mathrm{\Delta }{\gamma }_{\mathit{mn}}^i}{{\sum }_{i=1}^{n_{\mathrm{f}}}{W}_{\mathit{mn}}^i}·\end{array}$$(44)

Finally, we apply the smoothing Eq. (29) to the combined Δγ_mn and perform the remaining steps in Sect. 5.1 for a combined map of the excess mass.

For a combined measurement of the excess mass map with the convergence-stack technique in Sect. 5.3, we just take the equally weighted average of all grids ${\overline{\kappa }}_{\mathrm{emm}}^{\mathit{mn}}$ obtained from the individual n_f = 64 fields. An optimised weighting scheme is not necessary here because this simplistic approach already produces maps with negligible statistical noise for our simulated data.

5.5. Statistical errors

For an estimate of the statistical error in the CFHTLenS maps of the excess mass, we perform the jackknife technique (see e.g. Knight 1999). The basic idea is as follows. Let G(d) be an estimator for some quantity G based on the complete data set d. In our case, G(d) is our estimator for the smoothed excess mass map (or pair convergence map) ${\overline{\kappa }}_{\mathit{mn}}^{\mathrm{K}}$ at the grid pixel mn, and d comprises the merged catalogue of all separate fields. For the jackknife estimator, we split the complete data d = d₁ ∪ d₂ ∪ … ∪ d_njn into n_jn disjoint samples, and we compute a sample {G_−i : i = 1…n_jn} of estimates G_−i := G(d  ∖  d_i) based on d without the subset d_i. We define with

$$\begin{array}{c}\hfill {\sigma }^2(G)=\frac{n_{\mathrm{jn}}-1}{n_{\mathrm{jn}}}{\displaystyle \sum _{i=1}^{n_{\mathrm{jn}}}{(G_{°}-G_{-i})}^2,}\end{array}$$(45)

for

$$\begin{array}{c}\hfill G_{°}=\frac{1}{n_{\mathrm{jn}}}{\displaystyle \sum _{i=1}^{n_{\mathrm{jn}}}{G}_{-i},}\end{array}$$(46)

the estimator of the jackknife error-variance of G(d). For the following maps, we compute the jackknife error $\sigma ({\overline{\kappa }}_{\mathit{mn}}^{\mathrm{K}})$ for every pixel value ${\overline{\kappa }}_{\mathit{mn}}^{\mathrm{K}}$ and quote ${\overline{\kappa }}_{\mathit{mn}}^{\mathrm{K}}/\sigma ({\overline{\kappa }}_{\mathit{mn}}^{\mathrm{K}})$ for the S/N of a pixel value.

We expect that positional shot-noise and shape noise of the sources as well as cosmic variance are the dominating contributors of statistical noise in our measurement; see, for example, Kilbinger & Schneider (2005) for a discussion of statistical noise in related lensing correlation-functions. For the jackknife scheme, we remove individual pointings d_i from the merged catalogue; each CFHTLenS pointing has a square geometry and 1 × 1 deg² area. Since these pointings are significantly larger than our maps with typical angular scale of a few arcmin, we expect a sensible estimate for errors owing to cosmic variance at these scales (Shirasaki et al. 2017).

6. Results

In this section, we present maps of the excess mass ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm}}$ and pair convergence $\overline{\mathrm{\Delta }\kappa }$ for CFHTLenS lenses and for H15 galaxies in our synthetic data. We compute the maps for two photo-z bins and two angular separations of lenses on the sky: the separation bins 40″ ≤ θ₁₂ <  60″ (“close-θ”) and 60″ ≤ θ₁₂ <  80″ (“wide-θ”); and the two redshift bins 0.2 ≤ z_ph <  0.44 (“low-z”) and 0.44 ≤ z_ph <  0.6 (“high-z”). We only report E-mode maps in this section; the corresponding B-mode maps can be found in the Appendix. The B-modes are consistent with a vanishing signal. Moreover, we study the impact of a redshift slicing of lenses as means to reduce the fraction of chance pairs in a shear stack.

6.1. Code verification

For Fig. 5, we compare the reconstruction of the excess mass with synthetic data for two different methods: shear stacking (Sect. 5.1) and convergence stacking (Sect. 5.3). For convenience, in the high-z bin we show only the reconstructions for close-θ galaxies. The relative deviations of maps for other samples, including the maps of the pair convergence, are comparable to the single case shown here. The bottom panel employs, for ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},2}$, the shear-stacking that we apply to the CFHTLenS data. The top panel displays the map ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},1}$ using direct convergence stacking. Each of the 64 fields contains 2 × 10⁴ sources with no intrinsic shape noise, to reduce the noise in the maps for this comparison. Inside the panels, we indicate the lens positions ${\mathit{p}}_1^{\mathrm{d}}$ and ${\mathit{p}}_2^{\mathrm{d}}$ by black crosses. We overall find an excellent agreement for both approaches, with relative differences usually around five per cent or less, wherever the signal is larger than 5 × 10⁻⁴. However, the differences grow larger close to the lens positions inside the map, as indicated by the contours of $({\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},2}-{\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},1})/{\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},1}$, shown with the error levels in increments of 10%. Nevertheless, the errors are below 20% except within a couple of pixels separation from the lens positions. Presumably this error is owing to pixellation and binning of the correlation function ${\overline{\gamma }}_{\mathrm{t}}$. The 10% error contours at the edges of the map are likely just numerical noise, which becomes relevant for values of ${\overline{\mathrm{\Delta }\kappa }}_{\mathrm{emm},1}$ that are close to zero; the convergence stacking randomly picks a separation to an imaginary lens to subtract second-order correlations from the map. We reiterate that shear stacking is insensitive to constant offsets in κ so that this level of agreement here is only valid for a consistent definition of κ₀.