© D. Crespo et al. 2022

1. Introduction

Magnification bias, one of the effects of gravitational lensing, is the magnification of the apparent flux densities of background sources that suffer from lensing due to the presence of a foreground lens. Such a boost in the flux might be enough to make the background source detectable in cases where otherwise it would be below the instrument detection limit (e.g. Bonavera et al. 2021a). In order to see excess flux density, the logarithmic slope of the integrated source number counts of the background sources needs to be steep. Estimation of the magnification bias can be performed through the cross-correlation function (CCF), which equates to the excess expected between two source samples with non-overlapping redshift distributions with respect to the case with no magnification (Scranton et al. 2005; Ménard et al. 2010; Hildebrandt et al. 2013; Bartelmann & Schneider 2001). This effect was directly measured by Dunne et al. (2020) with the Atacama Large Millimetre Array (ALMA).

The submillimetre galaxies (SMGs) are considered (González-Nuevo et al. 2014, 2017) to be an optimal sample of background sources for magnification bias studies thanks to their steep source number counts, n(> S) = n₀ S^−β with β ≳ 2, and their high redshift, z > 1. Moreover, the foreground lens (detected in the optical band) emission at submillimetre wavelengths can be considered negligible, while the SMGs are invisible in the optical band for a similar reason. Therefore, both samples are not cross-contaminated by each other. For these reasons, SMGs have already been successfully used in magnification bias analyses to investigate the projected mass density profile and concentration of foreground samples of quasi-stellar objects (QSOs; Bonavera et al. 2019), for cosmological studies (Bonavera et al. 2020, 2021b; González-Nuevo et al. 2021), and to observationally constrain the halo mass function (Cueli et al. 2021, 2022).

With respect to the type of lens discussed in this work, QSOs are extremely luminous active galactic nuclei (AGN) that can be detected over a very broad range of distances. They are very suitable as background objects in lensing studies in general and also when adopting CCF measurements for magnification bias studies (e.g. Bartelmann & Schneider 1994; Scranton et al. 2005; Ménard et al. 2010). For example, the Sloan Digital Sky Survey (SDSS) Quasar Lens Search (SQLS) identified 28 galaxy-scale multiply imaged quasars that undergo gravitational lensing (Oguri et al. 2006, 2008). On the other hand, the first strong gravitational lensing caused by a QSO (SDSS J0013+1523 at z = 0.120) was detected by Courbin et al. (2010). Additional cases of QSOs acting as foreground objects have since been identified (e.g. Courbin et al. 2012; Harvey & Courbin 2015; Danielson et al. 2017). QSOs are also used in shear measurements as in Mandelbaum et al. (2009) and Luo et al. (2022), where the authors measure the weak-lensing shear distortions taking AGNs as lenses, making these studies complementary to magnification bias measurements. QSOs acting as lenses on SMGs have also been used to extract information on the mass density profile (see Bonavera et al. 2019) with the cross-correlation approach.

Large galaxy surveys (as SDSS, York et al. 2000) are undoubtedly very useful in cosmology, (e.g. BOSS Dawson et al. 2013). A first attempt at measuring lensing-induced cross-correlations between SMGs and low-z galaxies was carried out by Wang et al. (2011), who found convincing evidence of the effect. With much better statistics, this bias was studied in detail with the CCF between SMGs and massive galaxies by González-Nuevo et al. (2014, 2017). More recently, the magnification bias on SMGs produced by such massive optical galaxies at z ≪ 1 was analysed to derive complementary and independent constraints on the main cosmological parameters (Bonavera et al. 2020, 2021b; González-Nuevo et al. 2021).

Galaxy clusters, instead, are massive bound systems that can be used to track the large-scale structure of the Universe because they are generally placed in the knots of filamentary structures. They are also used for cosmological studies (e.g. Allen et al. 2011) to investigate the evolution of galaxies (Dressler 1980; Butcher & Oemler 1978, 1984; Goto et al. 2003) and the lensed high-redshift galaxies (e.g. Blain et al. 1999). Moreover, the correlation between galaxy clusters and selected background objects has been exploited to investigate possible lensing effects (Myers et al. 2005; Lopez et al. 2008) and, in particular, the CCF measurements using the SMGs as background sample to estimate the masses and concentration of the galaxy clusters (Fernandez et al. 2022).

Although the signal produced by weak lensing events is fainter than that produced by strong lensing events, the former are more frequent, allowing them to be studied by applying stacking techniques. Such techniques consist in co-adding the signal of weak or undetectable objects in order to enhance the signal and lower the background emission so that the overall statistical properties of this signal can be obtained. Stacking has been used in the recovery of the integrated signal of the Sachs–Wolfe effect with Planck data (Planck Collaboration XIX 2014; Planck Collaboration XXI 2016), and in studies of the faint polarised signal of radio and infrared sources (see Stil et al. 2014; Bonavera et al. 2017a,b). Moreover, the stacking technique has been used to obtain the mean spectral energy distribution (SED) of optically selected quasars (Bianchini et al. 2019) and to recover the weak gravitational lensing of the cosmic microwave background in the Planck lensing convergence map (Bianchini & Reichardt 2018). Stacking has also been used to probe star formation in dense environments of z ∼ 1 lensing halos (Welikala et al. 2016) and more recently in Fernandez et al. (2022) to recover the CCF signal between galaxy clusters for different richness ranges and SMGs.

In this work, the stacking technique is exploited to estimate and compare the mass density profiles of three different types of lenses (QSOs, galaxies, and galaxy clusters) and is organised as follows. Sections 2 and 3 present the data used and the methodology we are applying (both the stacking and the CCF estimation), respectively. Section 4 describes the theoretical framework of the CCF, weak gravitational lensing, and halo mass density profiles. The results are discussed in Sect. 5 and the conclusions are summarised in Sect. 6. The cosmological model used throughout the paper is the flat Λ cold dark matter (ΛCDM) with the cosmological parameters estimated by Planck Collaboration VI (2021, Ω_m = 0.31, σ₈ = 0.81 and h = H₀/100 km s⁻¹ Mpc⁻¹ = 0.67).

2. Data

2.1. Foreground samples

In this work, we estimate the CCF between two different lens samples (QSOs and galaxies) and the SMGs. Our initial sample of QSOs was obtained from that used in Bonavera et al. (2019). This QSO sample is selected from the Sloan Digital Sky Survey (SDSS), in particular, SDSS-II and SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) catalogues of spectroscopically confirmed QSOs detected over 9376 deg². We use the 7th¹ (DR7, Schneider et al. 2010) and 12th² (DR12, Pâris et al. 2017) SDSS data releases (see Ross et al. 2012, for a detailed discussion of the QSO target-selection process). The selection is based mainly on the DR7 catalogue, which mostly includes ‘low-z’ sources at z < 2.5. The DR12 sample specifically targeted QSOs at z > 2.5. Even so, there is a secondary maximum around z ∼ 0.8 as a consequence of colour degeneracy in the target selection of photometric data that leads to the observation of a fraction of low-redshift QSOs. Approximately 4% of the DR7 objects at z > 2.15 were re-observed for DR12, and so a combined sample was created joining the QSO information contained in the DR12 sample, namely any QSO present in both DR7 and DR12 is included only once.

In order to minimise the potential cross-contamination due to the redshift overlap between the foreground and background samples, we selected only QSOs with redshift between z = 0.2 and z = 1.0. This restriction leaves a total of 1546 QSOs in the common area with the background sample. The redshift distribution of the selected QSOs is shown in Fig. 1 through the blue histogram with a mean redshift $⟨z⟩=0.7_{-0.2}^{+0.1}$ (the uncertainty indicates the 1σ limits).

Fig. 1. Redshift distributions of the background sources from the H-ATLAS sample (red histogram), specifically the DR1 and NGP zones, and the QSO sample acting as lenses from the DR7 (Schneider et al. 2010) and DR12 (Pâris et al. 2017) catalogues (blue histogram).

On the other hand, the foreground galaxies sample was extracted from the GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014) survey, which was coordinated with the H-ATLAS one (our background sample; see following subsection; Pilbratt et al. 2010; Eales et al. 2010) to maximise the common area. Indeed, they both covered the three equatorial regions at 9, 12, and 14.5 h and part of the south Galactic pole, amounting to a common area of ∼207 deg².

The foreground galaxy sample consists of GAMA II sources within the interval 0.2 < z < 0.8, which results in approximately 102 672 galaxies (significantly higher than the case of QSO) with a mean spectroscopic redshift of $⟨z⟩=0.3_{-0.1}^{+0.1}$. Their redshift distribution is shown in Fig. 2 in blue. The red histogram corresponds to the background sample and is not exactly the same as in Fig. 1 due to the different zones covered by both foreground samples.

Fig. 2. Redshift distributions of the background sources from the H-ATLAS sample (red histogram), only the DR1 zone, and the galaxies sample acting as lenses from the GAMA catalogue (Liske et al. 2015, blue histogram).

2.2. Background sample

The background sample used in this work is the same for all the different foreground samples, and consist of the officially detected galaxies from Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS) data. These observations were collected by the Herschel Space Observatory (Pilbratt et al. 2010), cover about 610 deg², and are composed of five different fields. Three of these fields are taken from the first data release on the celestial equator, known as Galaxy and Mass Assembly (GAMA) fields or Data Delivery 1 (DR1, Valiante et al. 2016; Bourne et al. 2016; Rigby et al. 2011; Pascale et al. 2011; Ibar et al. 2010), and correspond to equatorial regions at 9, 12, and 14.5 h.

The other fields correspond to the North and South Galactic Poles, and are referred to as Data Delivery 2 (DR2), or NGP and SGP, respectively (Smith et al. 2017; Maddox et al. 2018). These latter two have total areas of 180.1 deg² and 317.6 deg², respectively. We use the three H-ATLAS GAMA fields (G09, G12 and G15) for both foreground samples. In the case of the QSO sample, only the NGP field is covered and only the SGP field is covered in the case of the galaxy sample. However, we decided to discard the SGP measurements due to the large uncertainty, which is much larger than that estimated in the GAMA fields. The main reason is that the galaxy distribution is very inhomogeneous within the field. Although taking into account this effect in the random simulations following González-Nuevo et al. (2021) improved the results, we considered the results from this field to be model dependent and not sufficiently robust to be considered for further analysis.

The H-ATLAS is composed of two instruments, the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010), which work in five photometric bands: 100, 160, 250, 350, and 500 μm. In both H-ATLAS DRs, there is an implicit 4σ detection limit at 250 μm (∼S₂₅₀ ≳ 29 mJy). The 1σ noise for source detection, including both confusion and instrumental noise, is 7.4 mJy at 250 μm (Valiante et al. 2016; Maddox et al. 2018). In addition, as in González-Nuevo et al. (2017), a 3σ limit at 350 μm was applied to increase the robustness of the photometric redshift estimation.

To avoid any overlap in the redshift distribution between the lens and background samples, only those SMGs with photometric redshift between 1.2 and 4.0 are taken into account. Such redshifts were estimated by taking the mean of a minimum χ² fit of a template SED to SPIRE data and also using PACS data when possible. The SED of SMM J2135–0102 (‘The Cosmic Eyelash’ at z = 2.3; Ivison et al. 2010; Swinbank et al. 2010) was shown to be a good template; it was found to be the best overall template with Δz/(1 + z) = − 0.07 and a dispersion of 0.153 (Ivison et al. 2016; González-Nuevo et al. 2012; Lapi et al. 2011). Following this selection process, we have 49 293 and 27 517 SMGs for the QSOs and galaxy sample analysis, respectively. They constitute approximately 21 and 23% of the corresponding initial background samples (the sources in the common areas between the SMGs sample and the QSOs and galaxy samples, respectively).

The redshift distribution of the background samples (in red) is shown in Fig. 1 for the QSO sample and Fig. 2 for the galaxy sample. Although the associated redshift distribution is not exactly the same for both cases, they have the same mean redshift $⟨z⟩=2.2_{-0.5}^{+0.4}$ (the uncertainty indicates the 1σ limits).

3. Measurements

3.1. Stacking

The technique known as ‘stacking’ is a method that can be used to determine the mean flux density of a large set of sources that are too weak to be analysed individually (Dole et al. 2006; Marsden et al. 2009; Béthermin et al. 2012). This tool is of value in our particular science case because it provides us with valid statistical information in cases where the noise is comparable to the signal, preventing detection. It consists in adding up patches of interest of the sky in order to enhance the signal, which would be otherwise undetectable in single events, not only for statistical flux density measurements.

A slight variation of this method was previously used by Bonavera et al. (2019) and Fernandez et al. (2022) to study the CCF signal due to magnification bias between SMGs (as the background sample) and QSOs and galaxy clusters (as the foreground sample), respectively. In these cases, the authors stacked the position of the sources around the lenses, not their flux densities, because the signal of interest is the number of background sources near the lens positions. Bonavera et al. (2019) derived the stacked magnification bias of lensed SMGs in lens positions signposted by QSOs and Fernandez et al. (2022) study the stacked magnification bias produced by galaxy clusters acting as lenses on the background SMGs.

The same procedure is applied in this work and follows a similar approach to the traditional CCF estimator, but with the advantage that it accounts for positional errors and can be used to identify the foreground–background pairs in the stacked map. Here we use this procedure to search for background sources in a circular area centred on the lens position. The search radius is set to be 200 arcsec, which means we must build a square map of 0.5 arcsec pixel size and 400 × 400 pixels centred at the lens position (QSO or galaxy, depending on the case being analysed) containing the paired background sources to the lens. The maps obtained for all the lenses are added up and normalised to the total number of lenses (102 672 targets for the galaxy sample and 1546 for the QSO sample) to produce the final stacked map. To account for the positional accuracy of the catalogues, a Gaussian filter of σ = 2.4 arcsec is applied. This value corresponds to the positional accuracy of the H-ATLAS catalogue (SMGs; Bourne et al. 2016; Maddox et al. 2018) and is more than an order of magnitude greater than that for the QSO or galaxy samples.

The resulting maps for QSO and galaxy samples are shown in the left and right panels of Fig. 3, respectively. The colour scale indicates the relative excess probability with respect to the random mean value (stacked pairs/random mean-1) of finding a background–foreground pair. The mean and standard deviation per pixel of the random stacked image are 2.25 × 10⁻⁶ and 1.21 × 10⁻⁶ for QSOs and 2.32 × 10⁻⁶ and 1.15 × 10⁻⁶ for the galaxies. Figure 4 shows the expected signal in the absence of lensing where QSOs and galaxy positions are randomly simulated (left and right panel, relatively) following the same procedure as for real data. The random stacked maps were simulated using approximately 6000 and 200 000 targets for the QSOs and galaxy samples for each GAMA region and 14 000 targets for the NGP region for QSOs (ten times the round-up value of the available possible foreground lenses).

Fig. 3. Relative excess probability (stacked image/random mean-1) for QSO (left) and galaxy (right) pairs considering the background sources within an angular radius of 200 arcsec from the position of the target. In both cases, the pixel size is 0.5 arcsec and we apply a 2.4σ Gaussian filter to take into account the positional uncertainties (see text for more details). The mean and standard deviation per pixel of the random stacked image for QSOs are 2.25 × 10−6 and 1.21 × 10−6, respectively, whereas the mean and standard deviation per pixel of the random stacked image for the galaxy sample are 2.32 × 10−6 and 1.15 × 10−6, respectively.

Fig. 4. Relative excess probability (stacked image/mean-1) of random pairs for QSOs (left) and galaxies (right) using the same radius and pixel size as for the data case and smoothed with the same Gaussian filter.

3.2. Cross-correlation function

We applied the CCF estimator as in Bonavera et al. (2019) to the stacked images in order to extract statistical information about the lens sample and the lensing system in general. The adopted estimator is that provided by Davis & Peebles (1983):

$$\begin{array}{c}\hfill {\stackrel{\sim }{w}}_x(\theta )=\frac{\mathrm{DD}}{\mathrm{RR}}-1,\end{array}$$(1)

where the measurements are computed by drawing a set of circles the centres of which are at the central position of the map and the radii of which logarithmically increase by 0.05 starting from 1 arcsec. The values of the pixels contained in the initial circle and in each one of the rings are added up as DD and the same procedure is repeated for the random map to produce RR.

To compute the CCF errors, each one of the first four rings is divided into sections of equal area with no more than eight pixels per section. The following rings instead have enough pixels to be split in 15 sections of equal area. Following Bonavera et al. (2019), in each ring, a Jackknife method is applied to estimate the DD and RR uncertainties, which are then propagated according to Eq. (1) to get the errors on ${\stackrel{\sim }{\mathit{w}}}_x$. Moreover, we analyse the possible correlation between the data points by computing the covariance matrix when performing the Jackkinfe analysis; this results in a negligible correlation.

Figures 5 and 6 show the estimated CCFs using stacking (black dots) and the CCF estimated following González-Nuevo et al. (2017, red dots) for the QSOs and galaxy samples, respectively. The angular scales range from 1 to 1000 arcsec, allowing us to study the mass density profile over a wide spectrum of physical scales (∼5 kpc to ∼4 Mpc for z = 0.3 or ∼8 kpc to ∼7 Mpc for z = 0.7): from the central region (mainly related with strong lensing events) to the outer part of the lens object (where weak lensing is most likely to happen). This is shown in Fig. 3 and measured in Figs. 5 and 6, where the CCF signal is stronger at the smaller angular scales and decreases towards the larger ones.

Fig. 5. Stacking data (black points) and those of the cross-correlation function (red points) computed using the QSOs acting as lenses as in Bonavera et al. (2019). On the left, the black line represents the NFW+SIS fits, where the corresponding SIS (green dashed line) and NFW (blue dashed line) fits are also separately shown. In the same way, on the right, the NFW+NFW profile is represented by the black line. In this case, the green dashed line corresponds to the fit to the points at the small scales only and the blue dashed line is for large scales. The grey dashed vertical line at ∼10 arcsec represents the separation between small and large angular scale regimes. Grey points are considered outliers and are not taken into account for the analysis.

Fig. 6. Stacking data (black points) and those of the cross-correlation function (red points) computed using the galaxies acting as lenses as in González-Nuevo et al. (2017). On the left, the black line represents the NFW+SIS fits, where the corresponding SIS (green dashed line) and NFW (blue dashed line) fits are also shown separately. In the same way, on the right, the NFW+NFW profile is represented with the black line. In this case, the green and blue dashed lines correspond to the fits to the points at the small scales only and large scale only, respectively. The grey dashed vertical line at ∼20 arcsec represents the separation between small and large scales. Grey points are considered outliers and are not taken into account for the analysis.

4. Theoretical framework

The cross-correlation measurements described above are strongly related to the physical properties of the mass distribution of the lenses. Indeed, the magnification field due to the density profile of a lens can be computed and linked directly to the cross-correlation signal by interpreting the physical meaning of the latter.

4.1. Gravitational lensing and the cross-correlation function

One of the consequences of the phenomenon of gravitational lensing is the modification of the integrated number counts of a sample of background sources due to a mass distribution between this latter sample and the observer. Indeed, due to the competing effects of a magnification that enhances fainter sources and a dilution that enlarges the solid angle in the sky, the number of background sources per solid angle and redshift with observed flux density larger than S is modified at every two-dimensional angular position θ in the sky as (Bartelmann & Schneider 2001)

$$\begin{array}{c}\hfill n(>S,z;\mathit{\theta })=\frac{1}{\mu (\mathit{\theta })}{n}_0(>\frac{S}{\mu (\mathit{\theta })},z),\end{array}$$

where μ(θ) is the magnification field and n₀ denotes the integrated number counts in the absence of lensing. If the (unlensed) background source number counts are assumed to follow a redshift-independent power-law behaviour, that is, n₀(> S, z) = A S^−β, then it is clear that

$$\begin{array}{c}\hfill \frac{n(>S,z;\mathit{\theta })}{n_0(>S,z)}={\mu }^{\beta -1}(\mathit{\theta }).\end{array}$$

The connection between the magnification field of a sample of lenses and the cross-correlation observable that we aim to measure can be elucidated when we understand the physical meaning of the above ratio. Indeed, the quantity n(> S, z_b; θ)/n₀(> S, z_b) represents the excess (or lack) of background sources (in direction θ as viewed by a lens at z_l) at redshift z_b > z_l with respect to the case without lensing. The angular CCF between a sample of foreground lenses at redshift z_l and a sample of background objects at redshift z_b is defined as

$$\begin{array}{c}\hfill w_x(\mathit{\theta };z_{\mathrm{l}},z_{\mathrm{b}})\equiv ⟨\delta n_{\mathrm{f}}(\mathit{\varphi })\delta n_{\mathrm{b}}(\mathit{\varphi }+\mathit{\theta })⟩,\end{array}$$

where δn_b and δn_f denote the background and foreground object density contrast, respectively. It follows from the above argument that

$$\begin{array}{c}\hfill w_x(\mathit{\theta };z_{\mathrm{l}},z_{\mathrm{b}})={\mu }^{\beta -1}(\mathit{\theta })-1,\end{array}$$

as we are stacking the lenses at a fixed position, which we take as the origin. We thus obtained a relation between our observable and a quantity that depends on the physical properties of the mean density profile of the lenses.

4.2. Mass density profiles

If a lens located at an angular diameter distance D_d (from the observer) deflects the light coming from a source at an angular diameter distance D_s, the convergence field at an angular position θ = xi/D_d on the image plane is defined as

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })=\frac{\mathrm{\Sigma }(D_{\mathrm{d}}\mathit{\theta })}{{\mathrm{\Sigma }}_{\mathrm{cr}}},\end{array}$$

where Σ(ξ) denotes the mass density projected onto a plane perpendicular to the light ray and

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{cr}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{d}}{D}_{\mathrm{ds}}}\end{array}$$

is the critical mass density, where D_ds is the angular diameter distance between the lens and the background source.

In the case of axially symmetric lenses, choosing the origin as the symmetry centre yields κ(θ) = κ(θ) and the magnification field μ(θ) is related to the convergence via (Bartelmann & Schneider 2001)

$$\begin{array}{c}\hfill \mu (\theta )=\frac{1}{(1-\overline{\kappa }(\theta )(1+\overline{\kappa }(\theta )-2\kappa (\theta ))},\end{array}$$

where $\overline{\kappa }(\theta )$ is the mean surface mass density within the angular radius θ. We now present the expressions of the magnification field for the two models we consider for the mean density profile of the lenses.

4.2.1. Navarro-Frenk-White profile

The Navarro-Frenk-White (NFW, Navarro et al. 1996) profile is a two-parameter model given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{NFW}}(r;r_{\mathrm{s}},{\rho }_{\mathrm{s}})=\frac{{\rho }_{\mathrm{s}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2},\end{array}$$

where r_s and ρ_s are the so-called scale radius and density, respectively. If halos are identified as spherical overdense regions with a mean density value of ρ_h, then r_s and ρ_s are related via

$$\begin{array}{c}\hfill {\rho }_{\mathrm{s}}=\frac{{\rho }_{\mathrm{h}}}{3}\frac{C^3}{ln(1+C)-C/(1+C)},\end{array}$$

where C ≡ r_h/r_s is referred to as the concentration parameter and r_h is the (truncation) radius of the identified halo.

It can be shown that the NFW profile satisfies (Schneider et al. 2006)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{NFW}}(\theta )=\frac{2r_{\mathrm{s}}{\rho }_{\mathrm{s}}}{{\mathrm{\Sigma }}_{\mathrm{cr}}}f(\theta /{\theta }_{\mathrm{s}}){\overline{\kappa }}_{\mathrm{NFW}}(\theta )=\frac{2r_{\mathrm{s}}{\rho }_{\mathrm{s}}}{{\mathrm{\Sigma }}_{\mathrm{cr}}}h(\theta /{\theta }_{\mathrm{s}}),\end{array}$$

where θ_s ≡ r_s/D_d is the angular scale radius,

$$\begin{array}{c}\hfill f(x)\equiv \{\begin{array}{cc}\frac{1}{x^2-1}-\frac{arccos(1/x)}{{(x^2-1)}^{3/2}}\hfill & \mathrm{if}x>1\hfill \\ \frac{1}{3}\hfill & \mathrm{if}x=1\hfill \\ \frac{1}{x^2-1}+\frac{\mathrm{arccosh}(1/x)}{{(1-x^2)}^{3/2}}\hfill & \mathrm{if}x<1\hfill \end{array}\end{array}$$

and

$$\begin{array}{c}\hfill h(x)\equiv \{\begin{array}{cc}\frac{2}{x^2}(\frac{arccos(1/x)}{{(x^2-1)}^{1/2}}+log\frac{x}{2})\hfill & \mathrm{if}x>1\hfill \\ 2(1-log2)\hfill & \mathrm{if}x=1\hfill \\ \frac{2}{x^2}(\frac{\mathrm{arccosh}(1/x)}{{(1-x^2)}^{1/2}}+log\frac{x}{2})\hfill & \mathrm{if}x<1\hfill \end{array}.\end{array}$$

4.2.2. Singular isothermal sphere profile

The singular isothermal sphere (SIS) model is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{SIS}}(r)=\frac{{\sigma }_{\mathrm{v}}^2}{2\pi G{r}^2},\end{array}$$

corresponding to a system of particles with a Maxwell velocity distribution at every radius with one-dimensional velocity dispersion σ_v, given by

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{v}}=\sqrt{\frac{\mathit{GM}}{2r_{\mathrm{h}}}},\end{array}$$

where M is the mass of the halo described by the profile. The convergence and mean surface density inside θ are easily shown to be (Schneider et al. 2006)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{SIS}}=\frac{{\theta }_{\mathrm{E}}}{2|\theta |}{\overline{\kappa }}_{\mathrm{SIS}}(\theta )=\frac{{\theta }_{\mathrm{E}}}{|\theta |},\end{array}$$(2)

where

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E}}=4\pi {\left(\frac{{\sigma }_{\mathrm{v}}}{c}\right)}^2\frac{D_{\mathrm{ds}}}{D_{\mathrm{s}}}\end{array}$$(3)

is the so-called Einstein radius of the model.

5. Results

As in Fernandez et al. (2022), with a view to analyse the measured CCFs and to extract physical conclusions, the data have been fitted by a combination of the two mass density profiles introduced above. The different fits to the data produced in this work clearly show that a single mass density profile is unable to fit the data at all scales. On the whole, using a combination of two common mass density profiles gives us a better description of the results.

In particular, in the following sections we present our results with CCF measurements using QSOs (Sect. 5.1) and galaxies (Sect. 5.2) as lenses. The background objects are from the same SMG sample for both cases.

5.1. Quasi-stellar objects

Although Bonavera et al. (2019) have already performed an analysis of CCF results using stacking with QSOs as lenses, we are driven to repeat this analysis in light of the fact that Bonavera et al. (2019) and Fernandez et al. (2022) use slightly different theoretical frameworks. The Bonavera et al. (2019) analysis with QSOs relies on a weak-lensing approximation which is no longer valid at small angular separations (the strong lensing regime). Fernandez et al. (2022), on the other hand, attempted to include the strong lensing regime in their study of galaxies clusters and the weak-lensing approximation is removed from their theoretical framework. In the present work, we are applying this latter approach – which does not require any kind of approximation and is therefore valid for all the lensing regimes – in order to study and compare the cases of galaxy clusters, galaxies, and QSOs acting as lenses, including those in the small angular scales regime. Therefore a re-analysis of QSOs is needed in order to apply a consistent methodology for all three cases.

The results obtained with stacking are shown in Fig. 5 with black dots. The red points are the CCF measurements derived with the traditional estimator described in González-Nuevo et al. (2017) (which does not make use of the stacking technique). Overall, there is good agreement between the two CCF estimations, with larger uncertainties at larger scales for the stacking approach. The fact that the CCF with no stacking is more robust in the low-signal regime was pointed out by Bonavera et al. (2019). The dotted vertical line indicates the two angular separation regions that divide the analysis: the inner part that extends until 10 arcsec and the outer part for angular scales greater than 10 arcsec. This particular angular separation was chosen as representative of the transition between the weak and strong lensing regimes. The stacking estimations have larger error bars in the outer region than in the inner one. In particular, this means that the stacking measurements (black points) are less informative above 100 arcsec than the CCF measurements (red points). For this reason, at angular scales above 100 arcsec we rely only on the red points to perform our analysis. Moreover, we exclude the CCF points at the largest scales from the analysis (grey points in Figs. 5 and 6, θ ≳ 200 arcsec or ∼1.4 Mpc at z ∼ 0.7 for the QSOs and θ ≳ 300 arcsec or ∼1.3 Mpc at z ∼ 0.3 for the galaxies), because they correspond to the two-halo term. In the case of QSOs, given the shortage of points above 100 arcsec we opt to keep the point at θ ∼ 250 arcsec in order to perform the fit, even if it is placed at a scale where the two-halo term might already come into place. Therefore, the estimated concentration value can be considered a lower value of the real concentration. Larger samples, which imply smaller uncertainties, and a proper two-halo term contribution estimation will be needed to improve the accuracy on the concentration at larger scales.

The left panel of Fig. 5 shows the combined SIS and NFW profiles with the black solid line. The corresponding individual NFW and SIS profiles are also shown with dashed red and gray lines. In the overall fit is good, with the SIS part fitting the inner region and the NFW fitting the outer region. This behaviour was not imposed by the methodology, but comes naturally. The major disagreement is in the region between 6 and 10 arcsec where the fit is below the data, as in Fernandez et al. (2022).

Table 1 summarises our results: from left to right, the columns are the profile, the characteristic halo mass and its error, and the concentration C and its error. In the inner region we obtain a mass value of log₁₀(M/M_⊙) = 12.934 ± 0.001 for the SIS profile. For the outer part (using a NFW profile) the results are log₁₀(M/M_⊙) = 13.44 ± 0.17 and C = 0.36 ± 0.1804, with a negligible correlation between the two parameters (3%).

It the right panel of Fig. 5, the fit is obtained with two independent NFW profiles (black solid line). In this case, the fit perfectly explains the measurements, although there is a hint of a small overfit of the data in the region between 10 and 30 arcsec even if still compatible with the large error bars. The NFW fit in the inner region gives a mass value of log₁₀(M/M_⊙) = 13.51 ± 0.04 and a concentration of C = 6.85 ± 0.34, with a negligible correlation between the two parameters (1%). The outer fit using a NFW profile is the same as in the previous case.

As a comparison, Bonavera et al. (2019) estimated a total mass of $M=1.7_{-0.5}^{+2.1}\times {10}^{14}{M}_{\odot }$ using the NFW profile that is three to four times our current estimates. This difference confirms the importance of improving the theoretical background in order for the analysis to function without any approximation, and demonstrates the effectiveness of the methodology, namely using two different independent mass density profiles for the inner and outer regions.

5.2. Galaxies

We repeated the same analysis for the sample of galaxies acting as lenses. However, this time the separation between the inner and outer angular separation regions is set at 20 arcsec (indicated as a dotted vertical line), because the mean redshift of this sample is different from that of QSOs (⟨z⟩∼0.7 for QSOs and ∼0.3 for the galaxies) implying different angular separations for roughly the same physical transition scale, ∼80 kpc.

As in the case of QSOs, with the stacking technique (black points in Fig. 6) we obtain large error bars for the outer region that limit our measurements above 50 arcsec. We therefore use the CCF measurements (red points) to better constrain the outer fit. The grey points are those excluded from the analysis. There is good agreement between the measurements obtained with the stacking and the traditional CCF approaches, especially at the smallest scales (which cannot be fully appreciated in our analysis of QSOs) even if no smoothing is applied to the CCF. This is confirmation of the necessity to introduce the smoothing step in the stacking procedure and of the positional uncertainty value used.

We again perform a first fit to the cross-correlation data, combining a SIS and an NFW profile. The results are shown in the left panel of Fig. 6 with the black solid line. The dashed red and grey lines show the parts of the fit corresponding to the NFW and the SIS, respectively.

As can be seen in the left panel of Fig. 6, the SIS+NFW does not produce a good fit to the measurements in both regions. The best fit is clearly below the measurements in the inner region, approximately between 5 and 10 arcsec, as in Fernandez et al. (2022), and is above them in the outer region. The parameter values for the profiles are listed in Table 2. The mass obtained for the inner SIS fit is log₁₀(M/M_⊙) = 12.784 ± 0.003 and the mass and concentration for the outer NFW fit are log₁₀(M/M_⊙) = 12.78 ± 0.21 and C = 1.21 ± 1.01, respectively. The correlation between the two parameters in this case is about 20%.

The right panel of Fig. 6 shows the fit resulting from the combination of two independent NFW profiles. Overall, this fit is much better than the fit with the SIS+NFW profile. There is still disagreement with the data in the region of transition between the inner and outer part at about 20 arcsec where the best fit is clearly above the measurements (which is similar to the results derived from the galaxy cluster sample analysis by Fernandez et al. 2022). As listed in Table 2, the value of the mass and concentration for the inner part are log₁₀(M/M_⊙) = 13.32 ± 0.08 and C = 8.23 ± 0.77, with a negligible correlation between the two parameters (6%). The results for the outer part are the same as the SIS+NFW fit.

5.3. Discussion

As described in Sects. 5.1 and 5.2, the measured CCFs cannot be explained using a single mass density profile. This mirrors the conclusions of Fernandez et al. (2022) and those of previous works (Bauer et al. 2014; Johnston et al. 2007; Okabe et al. 2016), but now with two additional independent samples with different kinds of lenses.

In this work, we follow the approach of Fernandez et al. (2022), fitting the measured CCFs for both samples with a combination of two different mass density profiles (SIS+NFW and NFW+NFW). Although both combinations provide a reasonable fit, the SIS+NFW profile predicts a lower cross-correlation at angular separations between ∼5 and 10 arcsec, which is more obvious for the galaxy sample (> 5σ difference). This issue is produced by the SIS mass density profile that is naturally chosen to fit the inner region. The same problem was also identified by Fernandez et al. (2022) using the galaxy clusters sample. Therefore, we discard the NFW+SIS profile and focus on the NFW+NFW profile for the rest of the discussion.

In the inner region, the estimated lens average mass is almost the same for both samples, with log₁₀(M/M_⊙) = 13.51 ± 0.04 for the QSOs and log₁₀(M/M_⊙) = 13.32 ± 0.08 for the galaxies. These masses are in agreement with those expected for a typical large red galaxy (LRG): an effective halo mass of M_eff = 2.8 − 4.4 × 10¹³ M_⊙ is derived from the modelling of the LRG angular correlation function (Blake et al. 2008) and consistent masses are also estimated from the analysis of their large-scale redshift-space distortions $M=3.5_{-1.4}^{+1.8}\times {10}^{13}{M}_{\odot }$, (Cabré & Gaztañaga 2009; Bauer et al. 2014). The estimated masses are also in agreement with the central mass of M_NFW = 3 − 4 × 10¹³ M_⊙ derived for the galaxy clusters and this agreement is more or less independent of richness (Fernandez et al. 2022). Therefore, we can conclude that our lenses are primarily LRGs, and that in the case of the galaxy clusters, these correspond to the brightest central galaxy (BCG), which has similar physical characteristics to the LRGs.

In the outer region, there is a clear difference between the two lens samples in their estimated lens average mass. For the galaxies, the derived outer mass, log₁₀(M/M_⊙) = 12.78 ± 0.21, is about five times lower than the inner mass. We believe this is because the galaxies acting as lenses in this sample tend to be isolated, or at least they are the central galaxy of a small group with dwarf galaxies as satellites. However, for the QSO sample, the estimated outer mass, log₁₀(M/M_⊙) = 13.44 ± 0.17, is similar to the inner mass. Therefore, we can expect the QSOs to be the central galaxies of a small group of smaller galaxies. Bonavera et al. (2019) came to the same conclusion after analysing the same QSO lens sample.

We can also compare the outer mass derived for the QSO sample with the masses obtained by Fernandez et al. (2022) for the different richness bins of the galaxy cluster sample (see Table 3). The estimated outer mass for the QSOs is five to six times smaller than that obtained for the lowest richness bin (12–17 members). Considering that the outer mass increases with the richness, this result implies that each of the QSOs is the central galaxy of a group of galaxies with less than about ten members, which is in agreement with our previous conclusion.

With respect to concentration, in the outer region, all the lens samples (QSOs, galaxies, and also galaxy clusters) show similar very low concentration values, namely C ≲ 1.2. As already noted in Fernandez et al. (2022), these values are generally lower than the ones retrieved from the most common mass–concentration relationships (e.g. Mandelbaum et al. 2008; Dutton & Macciò 2014; Child et al. 2018). However, Bonavera et al. (2019) estimate a concentration value of $C=3.5_{-0.3}^{+0.5}$ for the same QSO sample, which is in closer agreement with the mass–concentration relationships. The fact that the Bonavera et al. (2019) results are obtained by analysing the whole range of angular separations, which is the most common procedure, implies that the overall concentrations could be overestimated due to the influence of the highly concentrated central regions of the halos. Moreover, for the inner regions, we find much higher concentrations for both lens samples: C = 6.85 ± 0.34 for the QSOs and C = 8.23 ± 0.77 for the galaxies. These values are higher than the one estimated by Fernandez et al. (2022) for the lowest richness bin of the galaxy cluster sample, namely C = 3.63.

It is very interesting that the findings of Luo et al. (2022) for the satellite fraction can be extended to our case to explain the lower concentration values we obtain. In their work on dark-matter halos of luminous AGNs from galaxy–galaxy lensing, these latter authors modelled the halo occupation distribution (HOD) following Mandelbaum et al. (2005) and set the satellite fraction as a free parameter. They find that a larger satellite fraction leads to an upturn at larger scales and flattens the overall profile which translates in a lower concentration value. Furthermore, at relatively larger scales, the excess surface density drops down towards the profile of less massive galaxies, which is explained by the fact that the satellite fraction increases causing a flattening of the profile. Analogously, for the samples adopted in this work, we treat the relatively large scales (which still belong to the one-halo term) separately from the small scales, which enhances the flattening effect, because in this region we are isolating the satellites and obtaining low concentration values. This flatting is what drove us to split the fit into two parts to better describe such different behaviours between the two angular scale ranges.

Figure 7 shows the comparison between the cross-correlation measurements with stacking for the two different lens samples studied in this work. The data corresponding to the galaxy sample are represented with the red solid line, with the red coloured areas being the uncertainties, and those of the QSOs lens sample are shown with the blue solid line with their uncertainty represented by the blue coloured area. In addition, the cross-correlation obtained by Fernandez et al. (2022) for the three bins with the lowest richness is shown by the green solid lines (uncertainties are given by the green area). This figure illustrates our previous conclusions: the galaxies can mostly be considered isolated halos, the QSO measurements are similar to those for galaxy clusters in the outer region (i.e. they are the main members of a smaller group of galaxies), and the BCG halos of the lowest richness are less concentrated with respect the QSOs and galaxies. This last point suggests an evolutionary and dynamical variation of the inner region concentration. The stacking results describe the average characteristics, and not the individual ones, and therefore the following discussion does not apply to a single galaxy but should be considered as a plausible average scenario.

Fig. 7. Comparison of the different stacking results by means of continuous lines with their corresponding error bands. The stacking results obtained in this work are shown in blue for the case of QSOs and in red for the galaxy sample. Such results are compared to the cross-correlation data obtained in Fernandez et al. (2022), where the galaxy clusters are acting as lenses: in particular, the BIN 1 case is shown in dark green, BIN 2 in green, and BIN 3 in light green, corresponding to the 12–17, 18–25, and 26–40 richness intervals, respectively.

We can use the galaxy sample results as a benchmark, considering that the inner mass is similar to the other cases and assuming that they are mainly isolated halos. On the other hand, the QSOs show a lower concentration and, by definition, have an AGN in their nuclei. Moreover, they show a lower signal with respect to the galaxy clusters between 10 and 30 arcsec (∼70−210 kpc at z = 0.7), which corresponds to the outskirts of the host galaxy. Therefore, a plausible scenario to explain all these characteristics is that the host galaxy, as the BCG of a small group of galaxies, recently underwent an interaction with a surrounding satellite galaxy, which initiated the AGN episode.

On the other hand, the BCGs of the lowest richness bin of the galaxy clusters sample have similar mass but much lower concentration. In addition, the BCGs are mostly elliptical galaxies that are supposed to be the final product after the QSO phase. It is well established that the dynamical effects produced by an AGN affect the entire host galaxy (jets, outflows, star formation rate, etc.) and even its surroundings. Therefore, our results indicate that the BCGs of the smaller groups of galaxies could have less concentrated galaxy halos as a remnant of the QSO phase and AGN effect.

Moreover, Fig. 7 compares the cross-correlation function derived for the galaxy clusters in the different richness bins. As already discussed in Fernandez et al. (2022), the concentration increases with the mass or richness in our case. For the more massive galaxy clusters, which are older than the less massive ones, the BCG halos have had enough time to revert the AGN dynamical effects to a more concentrated mass density distribution due to gravity itself.

From this discussion, it is clear that a comparative analysis of mass density profiles estimated through the magnification bias for different kinds of lenses can be performed. Here, we demonstrated that this interesting and novel methodology can be used to study the average physical properties of the lenses and their surroundings and can potentially link these characteristics to the different galaxy evolution stages.

6. Conclusions

In this work, we exploit the magnification bias – a gravitational lensing effect produced on the SMGs observed by Herschel – in order to analyse the average mass density profiles of different types of lenses (QSOs, galaxies, and galaxy clusters). The magnification bias can be studied from the estimation of the CCF function, which we measure by stacking the background–foreground pairs with the Davis-Peebles estimator. One of the advantages of this approach is that positional uncertainties are easily taken into account, which are especially important at small angular scales. Moreover, with stacking data and including traditional CCF data at the largest angular scales, we can estimate the mass density profile for a large range of angular scales, ∼1 − 1000 arcsec, considering both weak and strong lensing effects and including the inner part < 100 kpc with high radial resolution. One the other hand, our results are average values obtained for samples and do not apply to single objects.

As already concluded by Fernandez et al. (2022), the CCF measurements cannot be explained with a single mass density profile. Therefore, we combine two different mass density profiles: we use an NFW for the outer part and we fit the inner part both to a SIS and an NFW profile. The SIS profile always underestimates the data in the central part, especially in the galaxy sample case, and the two NFW profiles produce a better fit. For this reason, we focus our discussion on the NFW+NFW case.

In the inner region, the average mass of the lenses is similar for both QSOs and galaxy samples log₁₀(M/M_⊙) = 13.32 and 13.51, respectively, and agrees with the values found in literature for typical LRGs. Moreover, this value also agrees with the average central mass estimated by Fernandez et al. (2022) for galaxy clusters, which is almost richness independent.

In the outer region instead, there is a clear difference between the masses estimated for the two samples. In the galaxy case, the outer estimated mass is 100 times lower than the inner mass, which is probably due to the fact that these galaxies are isolated or are the central galaxies of galaxy groups with a few dwarf galaxies. On the other hand, the QSO sample provides outer estimated masses that are similar to those of the inner part. A comparison with the results of Fernandez et al. (2022) as a function of the richness of the galaxy clusters suggests that the QSOs are the central galaxies of galaxy groups (as concluded in Bonavera et al. 2019) with less than ten members.

Our results show that the concentrations of the different sources are very similar for the three samples (QSOs, galaxies, and galaxy clusters) in the outer region but are lower than the majority of values reported in literature, including those obtained by Bonavera et al. (2019) for the same QSO sample. This can be explained by the fact that these latter studies consider only one mass density profile for all angular scales, including the more concentrated central halo regions. This explanation is also supported by our higher concentration estimations in the small-scale region. This is also analogous to the findings of Luo et al. (2022) for the satellite fraction: the increase in the satellite fraction at relatively large scales flattens the profile, which leads to a lower concentration. In our case, where the fit is done by isolating the data points at such scales, the effect on concentration is magnified.

Finally, the results from both angular scale regimes considered in our work suggest a scenario where the QSO host galaxies recently interacted with a close satellite galaxy, probably starting the AGN activity. Furthermore, the lower concentration values for the BCG halos in the galaxy clusters (mainly for the lowest richness bin) could be interpreted as the effects of the previous evolutionary activity produced during the earlier QSO and AGN phases acting on the surrounding environment. However, the galaxies in the galaxy sample and the BCGs in the more massive galaxy clusters, which are older than the less massive ones, had had enough time to weaken the AGN dynamical effects and drift back their halo into a more concentrated mass density distribution due to gravity itself.

Acknowledgments

We would like to warmly acknowledge the anonymous referee for the very useful comments that largely improved the work. D.C., J.G.N., L.B., M.M.C., J.M.C. acknowledge the PID2021-125630NB-I00 project funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. M.M.C. acknowledges PAPI-20-PF-23 (Universidad de Oviedo). The Herschel-ATLAS is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The H-ATLAS web-site is http://www.h-atlas.org. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. This research has made use of the python packages ipython (Pérez & Granger 2007), matplotlib (Hunter 2007) and Scipy (Jones et al. 2001).

References

All Tables

All Figures

	Fig. 1. Redshift distributions of the background sources from the H-ATLAS sample (red histogram), specifically the DR1 and NGP zones, and the QSO sample acting as lenses from the DR7 (Schneider et al. 2010) and DR12 (Pâris et al. 2017) catalogues (blue histogram).
In the text

	Fig. 2. Redshift distributions of the background sources from the H-ATLAS sample (red histogram), only the DR1 zone, and the galaxies sample acting as lenses from the GAMA catalogue (Liske et al. 2015, blue histogram).
In the text

Fig. 3. Relative excess probability (stacked image/random mean-1) for QSO (left) and galaxy (right) pairs considering the background sources within an angular radius of 200 arcsec from the position of the target. In both cases, the pixel size is 0.5 arcsec and we apply a 2.4σ Gaussian filter to take into account the positional uncertainties (see text for more details). The mean and standard deviation per pixel of the random stacked image for QSOs are 2.25 × 10−6 and 1.21 × 10−6, respectively, whereas the mean and standard deviation per pixel of the random stacked image for the galaxy sample are 2.32 × 10−6 and 1.15 × 10−6, respectively.

In the text

	Fig. 4. Relative excess probability (stacked image/mean-1) of random pairs for QSOs (left) and galaxies (right) using the same radius and pixel size as for the data case and smoothed with the same Gaussian filter.
In the text

Fig. 5. Stacking data (black points) and those of the cross-correlation function (red points) computed using the QSOs acting as lenses as in Bonavera et al. (2019). On the left, the black line represents the NFW+SIS fits, where the corresponding SIS (green dashed line) and NFW (blue dashed line) fits are also separately shown. In the same way, on the right, the NFW+NFW profile is represented by the black line. In this case, the green dashed line corresponds to the fit to the points at the small scales only and the blue dashed line is for large scales. The grey dashed vertical line at ∼10 arcsec represents the separation between small and large angular scale regimes. Grey points are considered outliers and are not taken into account for the analysis.

In the text

Fig. 6. Stacking data (black points) and those of the cross-correlation function (red points) computed using the galaxies acting as lenses as in González-Nuevo et al. (2017). On the left, the black line represents the NFW+SIS fits, where the corresponding SIS (green dashed line) and NFW (blue dashed line) fits are also shown separately. In the same way, on the right, the NFW+NFW profile is represented with the black line. In this case, the green and blue dashed lines correspond to the fits to the points at the small scales only and large scale only, respectively. The grey dashed vertical line at ∼20 arcsec represents the separation between small and large scales. Grey points are considered outliers and are not taken into account for the analysis.

In the text

Fig. 7. Comparison of the different stacking results by means of continuous lines with their corresponding error bands. The stacking results obtained in this work are shown in blue for the case of QSOs and in red for the galaxy sample. Such results are compared to the cross-correlation data obtained in Fernandez et al. (2022), where the galaxy clusters are acting as lenses: in particular, the BIN 1 case is shown in dark green, BIN 2 in green, and BIN 3 in light green, corresponding to the 12–17, 18–25, and 26–40 richness intervals, respectively.

In the text