© ESO 2020

1. Introduction

Magnification bias is the apparent excess number of high redshift sources seen near low redshift structures, and it can be quantified using the cross-correlation between high and low redshift galaxies. Light rays coming from a distant source are deflected by the intervening gravitational field so that the apparent area of a given region of sky is stretched. The increased area decreases the apparent surface density of any background sources. The deflections also amplify the flux of the background sources, increasing their chances of being included in a flux-limited sample. The net effect is known as magnification bias, and it has been extensively described in the literature (see e.g. Schneider et al. 1992).

The cross-correlation between two source samples with non-overlapping redshift distributions is an unambiguous manifestation of magnification bias. The occurrence of such correlations has been tested and established in several contexts: an 8σ detection of cosmic magnification from the galaxy-quasar cross-correlation function (Scranton et al. 2005); a simultaneous detection of gravitational magnification and dust reddening effects due to galactic halos and large-scale structure (galaxy-quasar cross-correlation; Ménard et al. 2010); and a 7σ detection of a cross-correlation signal between z ∼ 3 − 5 Lyman-break galaxies and Herschel sources (Hildebrandt et al. 2013), among others.

The strength of the magnification bias depends on the size of gravitational deflection as light passes near low redshift galaxies, and since the deflection depends on the cosmological distances, and on the galaxy halo properties, we can use these measurements to set constraints on cosmological parameters. The aim of our current analysis is to investigate how our measurements of magnification bias can constrain parameters in the standard cosmological model.

The current standard cosmological model successfully accounts for anisotropies in the CMB (e.g. Hinshaw et al. 2013; Planck Collaboration XIII 2016; Planck Collaboration VI 2020), the accelerating expansion of the Universe observed from SNIa (e.g. Betoule et al. 2014), and Big Bang nucleosynthesis (e.g. Fields & Olive 2006), and also predicts some important features of large-scale structure (LSS) in the distribution of galaxies (e.g. Peacock et al. 2001) including baryon acoustic oscillations (BAOs; e.g. Ross et al. 2015). This means that independent and complementary constraints on the model parameters can be obtained by measuring such observables (e.g. Peacock & Dodds 1994). The consistency between different probes is a key reason that the model is considered to be so successful. However, with the increase in the quality and quantity of the measurements, some tensions and small-scale issues have arisen, such as the missing satellite problem (Moore et al. 1999; Klypin et al. 1999) and the cusp-core problem (Wechsler et al. 2002), that might indicate the necessity of modifications of the ΛCDM model.

One of the main tensions is in the value of the Hubble constant, H₀: local estimates suggest a relatively high value (74.03 ± 1.42 km s⁻¹ Mpc⁻¹, Riess et al. 2019), whereas CMB and BAOs analysis advise a relatively low value (67.4 ± 0.5 km s⁻¹ Mpc⁻¹Planck Collaboration VI 2020). Another important inconsistency lies in the usually degenerate relationship between the Ω_m and σ₈ parameters. Constraints from local observables, such as LSS in the galaxy distribution, tend to prefer lower values compared to those derived from high redshift observables such as the CMB (e.g. Heymans et al. 2013; Planck Collaboration XXIV 2016; Hildebrandt et al. 2017). This ∼2σ tension is discussed in Planck Collaboration VI (2020). With this in mind, any method that constrains cosmological parameters using a new set of observables is worth pursuing in the effort to resolve the tensions. At the moment our constraints from the magnification bias are relatively weak, but the use of new independent observables makes it a valuable new technique.

In this paper we use the cross-correlation function between low and high redshift galaxy samples to measure magnification bias. An optimal choice of foreground and background samples can enhance the cross-correlation signal, and in particular some features of the submillimetre galaxies (SMGs) make them close to the optimal background sample for lensing studies. Some of their most important characteristics, related to their exploitation in magnification bias studies, are summarised below.

Firstly, data collected before the advent of the European Herschel Space Observatory (Herschel; Pilbratt et al. 2010) and the South Pole Telescope (SPT, Carlstrom et al. 2011) suggested that the number density of SMGs drops off abruptly at relatively bright submillimetre flux densities (∼50 mJy at 500 μm), indicating a steep luminosity function and a strong cosmic evolution for this class of sources (e.g. Granato et al. 2004). Several authors have argued that the bright tail of the submillimetre number counts may contain a significant fraction of strongly lensed SMGs (Blain 1996; Negrello et al. 2007). The Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2010) and the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) are wide-field surveys (∼380 deg² and ∼550 deg², respectively) conducted by the Herschel satellite. Thanks to their sensitivity and frequency coverage both surveys have led to the discovery of several lensed SMGs (Negrello et al. 2017, and references therein). The selection of strongly lensed galaxies at these wavelengths is made possible by the (predicted) steep number counts of SMGs (Blain 1996; Negrello et al. 2007); in fact, almost all of those galaxies whose flux density has been boosted by an event of lensing can be observed above a certain threshold, namely ∼100 mJy at 500 μm (see Negrello et al. 2010; Wardlow et al. 2013).

Next, with Herschel data, Negrello et al. (2010) produced the first sample of five strongly lensed galaxies (SLGs) by means of a simple selection in flux density at 500 μm. Preliminary source catalogues derived from the full H-ATLAS were then used to identify the submillimetre brightest candidate lensed galaxies for follow-up observations with both ground-based and space telescopes to measure their redshifts (see Negrello et al. 2017, and references therein) with 80 SLG candidates and confirm their nature (Negrello et al. 2010; Bussmann et al. 2012, 2013; Fu et al. 2012; Calanog et al. 2014).

Using the same methodology, i.e. a cut in flux density at 500 μm, Wardlow et al. (2013) have identified 11 SLGs over 95 deg² of HerMES (Oliver et al. 2010), while more recently Nayyeri et al. (2016) have published a catalogue of 77 galaxy candidates at lenses with F₅₀₀ ≥ 100 mJy extracted from the HerMES Large Mode Survey (HeLMS; Oliver et al. 2012) and the Herschel Stripe 82 Survey (HerS; Viero et al. 2014) over an area of 372 deg². Altogether, the extragalactic surveys carried out with Herschel delivered a sample of more than a hundred submillimetre bright SLGs, which, as argued by González-Nuevo et al. (2012, 2019), might increase to over a thousand if the selection is based on the steepness of the luminosity function of SMGs (Lapi et al. 2011) rather than that of the number counts. Similarly, at millimetre wavelengths, the SPT survey has already discovered several tens of lensed galaxies (e.g. Vieira et al. 2013; Spilker et al. 2016) and other lensing events have been found in the Planck all-sky surveys (Cañameras et al. 2015; Harrington et al. 2016).

The properties of SMGs have been very well quantified by these previous analyses, and SMGs can be considered solid targets to be studied with lensing. In particular, the (sub)millimetre SLGs are very faint in the optical, whereas most foreground lenses are passive ellipticals (Auger et al. 2009; Negrello et al. 2010), essentially invisible at submillimetre wavelengths. This means that the foreground lens is “transparent” at (sub)millimetre wavelengths, i.e. it does not confuse the measurements of the background source and vice versa. In this way, the (sub)millimetre selection shares with spectroscopic searches the capability of detecting lensing events with small impact parameters, with the added advantage that in most cases there is no need to subtract the lens contribution to recover the source images within the effective radii of the lenses. A further advantage given by (sub)millimetre selection is that the typical source redshifts are above z >  1 − 1.5, and so the samples probe earlier phases of galaxy evolution, which typically have higher lensing optical depths. This makes the SMGs the perfect background sample for tracing the mass density profiles of galaxies and clusters, including baryonic and dark matter, and their time-evolution, using gravitational lensing.

A first attempt at measuring lensing-induced cross-correlations between Herschel/SPIRE galaxies and low-z galaxies was carried out by Wang et al. (2011), who found convincing evidence of the effect. With much better statistics, this possible bias was studied in detail in González-Nuevo et al. (2014) by measuring the angular cross-correlation function between selected H-ATLAS high-z sources, z >  1.5, and two optical samples with redshifts 0.2 <  z <  0.8, extracted from the Sloan Digital Sky Survey (SDSS; Ahn et al. 2012) and Galaxy and Mass Assembly (GAMA, Driver et al. 2011) surveys. The observed cross-correlation function was measured with high significance, > 10σ. Moreover, based on realistic simulations, it was concluded that the signal was entirely explained by a magnification bias produced by the weak lensing effect of low redshift cosmic structures (galaxy groups and clusters with halo masses in the range 10^13.2 − 10^14.5 M_⊙) signposted by the brightest galaxies in the optical samples.

González-Nuevo et al. (2017) used updated catalogues covering a wider area to make more accurate measurements of the cross-correlation function, with S/N >  5 above 10 arcmin and reaching S/N ∼ 20 below 30 arcsec. Thanks to the better statistics it was possible to split the sample into different redshift bins and to perform a tomographic analysis (with S/N >  3 above 10 arcmin and reaching S/N ∼ 15 below 30 arcsec). Moreover, a halo model was implemented to extract astrophysical information about the background galaxies and the deflectors that produce the lensing that links the foreground (lenses) and background (sources) samples. In the case of the sources, typical mass values were found to be in agreement with previous studies (10^12.1 − 10^12.9 M_⊙, Xia et al. 2012; Cooray et al. 2010; Mitchell-Wynne et al. 2012; Wilkinson et al. 2017). However, the lenses are massive galaxies or even galaxy groups and clusters, with minimum mass of M_lens ∼ 10^13.06 M_⊙, confirming the interpretation from González-Nuevo et al. (2014). The halo modelling also helped to identify a strong lensing contribution to the cross-correlation for angular scales below 30 arcsec. This interpretation is supported by the results of basic but effective simulations.

Following on from this, Bonavera et al. (2019) were able to measure the magnification bias produced by a sample of quasi-stellar objects QSOs acting as lenses, 0.2 <  z <  1.0, on the SMGs observed by Herschel at 1.2 <  z <  4.0, by applying the traditional cross-correlation function approach and the Davis-Peebles estimator through a stacking technique. Using the halo modelling of the cross-correlation function it was possible to estimate the halo mass where the QSOs acting as lenses are located: it was found to be greater than $M_{\min }>{10}^{13.6_{-0.4}^{+0.9}}{M}_{\odot }$, also confirmed by the mass density profile analysis (M_200c ∼ 10¹⁴ M_⊙). These mass values indicate that we are observing the lensing effect of a cluster size halo signposted by the QSOs.

In the rest of this paper, we describe our first attempts to estimate cosmological parameters using the cross-correlation function and data from González-Nuevo et al. (2017), and demonstrate the powerful potential of magnification bias as a cosmological tool. In particular, we focus on the estimation of the parameters Ω_m, σ₈, and H₀ in order to quantify what can actually be done with magnification bias, and what will be needed in the future to improve our cosmological analysis.

The work is organised as follows. In Sect. 2 the background and foreground samples are described and in Sect. 3 the methodology is presented, both for auto- and cross-correlation measurements. Our results and conclusions are discussed in Sects. 4 and 5, respectively. In Appendix A we show how the variations in the astrophysical and cosmological parameters affect the model behaviour and we compare the model with the data.

2. Data

2.1. Background sample

The background sample was selected from the H-ATLAS data set, the largest extra-galactic survey carried out by the Herschel space observatory (Pilbratt et al. 2010): its 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) cover about 610 deg² between 100 and 500 μm. The survey consists of five different fields: three located on the celestial equator (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) covering a total area of 161.6 deg² and two centred on the north and south galactic poles (NGP and SGP fields or DR2; Smith et al. 2017; Maddox et al. 2018) covering 180.1 deg² and 317.6 deg², respectively.

We selected our background sample from the sources detected in the three GAMA fields, covering a total area of ∼147 deg², and the part of the SGP that overlaps with the foreground sample (∼60 deg², see Sect. 2.2). This leads to ∼207 deg² of common area between the two samples.

It should be noted that an implicit 4σ detection limit at 250 μm (∼S250 >  29 mJy) was applied to both H-ATLAS DRs. At this frequency the estimated 1σ source detection total noise level (both instrumental and confusion) is 7.4 mJy for both DR1 and DR2 catalogues (Valiante et al. 2016; Maddox et al. 2018). In addition, a 3σ limit at 350 μm was introduced (González-Nuevo et al. 2017) to increase the robustness of the photometric redshift estimation.

In order to avoid overlap in the redshift distributions of lenses and background sources, only background sources with photometric redshift 1.2 < z < 4.0 were taken into account. Such redshifts were estimated with a minimum χ² fit of a template SED to the SPIRE data (using PACS data when possible). The adopted template is the SED of SMM J2135-0102 (“the Cosmic Eyelash” at z = 2.3; Ivison et al. 2010; Swinbank et al. 2010), proven to be the best overall template with Δz/(1 + z) = − 0.07 and dispersion of 0.153 (Ivison et al. 2016; González-Nuevo et al. 2012; Lapi et al. 2011). In the end, we are left with 57930 sources that correspond approximately to 24% of the initial sample.

Figure 1 shows the redshift distribution of the background sample in red, being $⟨z⟩=2.2_{-0.5}^{+0.4}$ the mean redshift of the sample (the uncertainty indicates the 1σ limits). This redshift distribution is the estimated p(z|W) of the galaxies selected by our window function (i.e. a top hat for 1.2 < z < 4.0). It takes into account the effect of random errors in photometric redshifts on the redshift distribution, as done in González-Nuevo et al. (2017).

Fig. 1. Redshift distributions of the background H-ATLAS sample (red) and the foreground GAMA sample (blue).

2.2. Foreground sample

As detailed in González-Nuevo et al. (2017), the foreground sources were drawn from the GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014; Liske et al. 2015) spectroscopic survey. Moreover, both H-ATLAS and GAMA II surveys were coordinated in order to maximise the common area coverage overlap. In particular, the three equatorial regions at 9, 12, and 14.5 h (referred to as G09, G12, and G15, respectively) were observed by both surveys, and the SGP region surveyed by H-ATLAS was partially covered by GAMA II. We are then left with a common area of about ∼207 deg² between the two surveys, which has been surveyed down to a limit of r ≃ 19.8 mag.

Our main sample consists of a selection of all GAMA II galaxies with reliable redshift measurements between 0.2 < z < 0.8. This main sample comprises a total of ∼150 000 galaxies, whose median redshift, z_spec, med = 0.28, is significantly lower than that of the background sample, as shown by the blue histogram in Fig. 1.

The lower limit was chosen because we do not expect a relevant gravitational lensing effect from mater density distribution at z < 0.2 for the background sample (see e.g. Lapi et al. 2012). For a similar reason, we can safely neglect any possible magnification bias affecting the lens sample: the lensing optical depth declines steeply for sources below z <  0.5.

3. Methodology

3.1. Halo occupation distribution (HOD)

Nowadays the power spectrum of the galaxy distribution is parametrised by the sum of a two-halo term, related to the correlation between different halos that dominates on large scales, and a one-halo term, more important on small scales, that depends on the distribution of galaxies within the same halo.

All halos above a minimum mass M_min host a galaxy at their centre, while any remaining galaxy is classified as satellite and such satellites are distributed proportionally to the halo mass profile (see e.g. Zheng et al. 2005). Halos host satellites when their mass exceeds the M₁ mass, and the number of satellites is a power-law (α) function of halo mass.

Therefore, the central galaxy occupation probability can be parametrised as a step function

$$\begin{array}{c}\hfill N_{\mathrm{cen}}(M_{\mathrm{h}})=\{\begin{array}{c}0\mathrm{if}{M}_{\mathrm{h}}<M_{\min }\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\end{array}$$(1)

and the satellite galaxies occupation as

$$\begin{array}{c}\hfill N_{\mathrm{sat}}(M_{\mathrm{h}})=N_{\mathrm{cen}}(M_{\mathrm{h}})·({\displaystyle \frac{M_{\mathrm{h}}}{M_1}}{)}^{{\alpha }_{\mathrm{sat}}},\end{array}$$(2)

where M_min, M₁, and α_sat are the free-parameters of the model.

3.2. Foreground angular auto-correlation function

In order to have an independent estimation of the HOD parameters that correspond to the foreground galaxies acting as lenses, we perform an auto-correlation analysis on the GAMA sample.

3.2.1. Estimation methodology

For this auto-correlation analysis, we focus on the three largest fields: G09, G12, and G15. They constitute ∼63% of the whole GAMA survey and they are the only fields with complete overlap with the H-ATLAS fields. We expect no relevant differences in their correlation properties compared to the rest of the survey. Moreover, we perform a further reduction in the foreground sample by selecting the normalising redshift quality (catalogue column labelled “NQ”) greater than 2 and the spectroscopic redshift (catalogue column labelled “Z”) in the range 0.05 <  z_spec <  0.6. The resulting median redshift of the sample of sources is z_med = 0.23.

We estimate the number density of sources by computing the volume of each field as V_field = A_field ⋅ Δd_C, where Δd_C = d_C(z_max)−d_C(z_min) is the difference in comoving distance between the maximum and minimum redshift selected. The area of the field is computed as A_field = d_hav(RA)⋅d_hav(Dec), where d_hav(RA) and d_hav(Dec) are the haversine distances

$$\begin{array}{c}\hfill d_{\mathrm{hav}}/r=arcsin(\sqrt{sin\left({\displaystyle \frac{{\varphi }_2-{\varphi }_1}{2}}\right)+cos{\varphi }_{1}cos{\varphi }_2\left({\displaystyle \frac{{\lambda }_2-{\lambda }_1}{2}}\right)})\end{array}$$(3)

along the RA and Dec directions, respectively.

In Table 1 we show the measures obtained for the three different selected fields with the cuts described above. The last row of the table shows the median value of n_g and maximum error value, σ_ng, among the three fields.

We compute the angular correlation function ω(θ) in each field by applying the modified Landy-Szalay estimator (Landy & Szalay 1993)

$$\begin{array}{c}\hfill w(\theta )={\displaystyle \frac{DD(\theta )-2DR(\theta )+RR(\theta )}{RR(\theta )}}\end{array}$$(4)

and estimate the error on this measure by dividing each field into 128 bootstrap samples (see Fig. 2).

Fig. 2. Measured angular auto-correlation function of the foreground sample (red circles). The best-fit halo model prediction is shown as a comparison (total, blue solid line; one-halo term, red dotted line; two-halo term, green dashed line).

3.2.2. Model

The traditional halo modelling (see e.g. Cooray & Sheth 2002) allows us to obtain estimates of the average surface density of sources and of the angular correlation function at a given redshift. The surface density is given by

$$\begin{array}{c}\hfill n_{\mathrm{g}}(z)\equiv {\displaystyle {\int }_{M_{\min }}^{M_{\max }}⟨N_{\mathrm{g}}⟩(M_{\mathrm{h}})n(M_{\mathrm{h}},z)\mathrm{d}{M}_{\mathrm{h}},}\end{array}$$(5)

where N_g(M_h) = N_cen(M_h)+N_sat(M_h) and where n(M_h, z) is the Sheth & Tormen mass function (Sheth & Tormen 1999). The angular correlation function is given by

$$\begin{array}{c}\hfill \omega (\theta ,z)={\displaystyle \int \mathrm{d}z{\displaystyle \frac{\mathrm{d}V(z)}{\mathrm{d}z}}{\mathcal{N}}^2(z)w[r_{\mathrm{p}}(\theta ,z),z],}\end{array}$$(6)

where w[r_p, z] is the projected two-point correlation function obtained using the Limber approximation, $\frac{\text{d}V(z)}{\text{d}z}$ is the comoving volume unit, and 𝒩(z) is the normalised redshift distribution of the target population. If we assume that w(r_p, z) is approximately constant in the redshift interval [z₁, z₂], we can then write

$$\begin{array}{c}\hfill \omega (\theta ,z)\approx [{\displaystyle {\int }_{z_1}^{z_2}\mathrm{d}z{\displaystyle \frac{\mathrm{d}V(z)}{\mathrm{d}z}}{\mathcal{N}}^2(z)]·w[r_{\mathrm{p}}(\theta ),\overline{z})]=A·w[r_{\mathrm{p}}(\theta ),\overline{z})]}\end{array}$$(7)

so that the normalisation factor, A, can be accounted for by adding another parameter to the model.

3.3. Angular cross-correlation function

For the cross-correlation measurements we use the data described in Sect. 2 and applied the same approach as in González-Nuevo et al. (2017) with slight changes, as described in the following section.

3.3.1. Estimation method

We use a modified version of the Landy & Szalay (1993) estimator (Herranz 2001),

$$\begin{array}{c}\hfill w_x(\theta )=\frac{D_1{D}_2-D_1{R}_2-D_2{R}_1+R_1{R}_2}{R_1{R}_2},\end{array}$$(8)

where D₁D₂, D₁R₂, D₂R₁, and R₁R₂ are the normalised data1-data2, data1-random2, data2-random1, and random1-random2 pair counts for a given separation θ.

In this work we do not apply exactly the same procedure used by González-Nuevo et al. (2017) and Bonavera et al. (2019). Instead, we follow the original estimation procedure of González-Nuevo et al. (2014), which divides the total common area between the two samples into several square mini-regions with minimal overlap. Therefore, the estimated cross-correlation function in each mini-region can be assumed to be statistically independent. In this way we can use the full available area (maximising the data statistics) without leaving portions of the mini-areas with no sources, as would be the case for the circular regions used by González-Nuevo et al. (2017). Taking into account the shape of the H-ATLAS areas (related with the satellite scanning strategy), we adopt square mini-regions of ∼4 deg². This division provided us with 62 usable almost independent mini-regions (15 in each GAMA field and 17 in the SGP field). We verified that using bigger square mini-regions (∼16 deg²) leads to results that are compatible within the uncertainties on large scales, but with higher uncertainties at intermediate scales due to the higher cosmic variance effect with only 16 independent mini-regions.

The estimated cross-correlation function is shown in Fig. 3 (red circles). Each measurement corresponds to the mean value of the cross-correlation functions estimated in each individual mini-region for a given angular separation bin. The uncertainties correspond to the standard error of the mean, i.e. ${\sigma }_{\mu }=\sigma /\sqrt{n}$ with σ the standard deviation of the population and n the number of mini-regions. For the cross-correlation represented by the red circles, the binning was chosen in the same way as in González-Nuevo et al. (2017), being the circle in Fig. 3 placed at the centre of the corresponding bin. This is the main data set, which can be directly compared with that used in González-Nuevo et al. (2017; shown in Fig. 3 as grey circles). The main difference between the two estimations are the two latest points at large angular scales that are slightly higher in the González-Nuevo et al. (2017) case. The much higher number of independent mini-regions makes us more confident in the current estimation with a more realistic uncertainties. In addition, we estimate again the cross-correlation function using the same procedure described above, but with a half-bin displacement (i.e. by using the same binning size as in González-Nuevo et al. 2017, but the bin centre is now shifted of half-bin size with respect to the previous measurement, shown in red). The cross-correlation measurements are shown as blue circles in Fig. 3. This additional set of results will provide us the opportunity to study potential dependence of our final results produced by the binning choice (important at large angular scales due to the much lower number of pairs in excess with respect to the random case and the cosmic variance).

Fig. 3. Cross-correlation data computed on Herschel SMGs (as background sample) and GAMA galaxies (as lenses sample). The grey points are exactly the results obtained in González-Nuevo et al. (2017). The red data set (obtained using the same binning as in González-Nuevo et al. 2017) and the blue data set (obtained using same the binning, but half-size bin shifted with respect to González-Nuevo et al. 2017) are the cross-correlation results computed for this work, exploiting a larger area. The solid lines represent the best fit with flat priors to the astrophysical and cosmological parameters using the default binning (red) and the shifted binning (blue). The green line is the best fit to the red case, but adopting Ωm = 0.3.

3.3.2. Model

Similar to the auto-correlation case, and as described in González-Nuevo et al. (2017) and summarised in Bonavera et al. (2019), in order to interpret a foreground-background source cross-correlation signal we adopt the halo model formalism from Cooray & Sheth (2002). We define halos as spherical regions whose mean overdensity with respect to the background at any redshift is given by its virial value, which is estimated following Weinberg & Kamionkowski (2003) assuming a ΛCDM model with Ω_m + Ω_Λ = 1. The density profile adopted is the one by Navarro et al. (1996; hereafter NFW) with the concentration parameter given in Bullock et al. (2001).

Weak lensing provides the connection between the low redshift galaxy-mass correlation and the cross-correlation between the foreground and background sources. The mass density field traced by the foreground galaxy sample causes weak lensing; the weak lensing field is imprinted onto the number counts of the background galaxy sample through magnification bias. Following mainly Cooray & Sheth (2002; see González-Nuevo et al. 2017, for details), we compute the correlation between the foreground and background sources adopting the standard Limber (Limber 1953) and flat-sky approximations (see e.g. Kilbinger et al. 2017, and references therein). It can be evaluated as

$$\begin{array}{cc}\hfill w_{\mathrm{fb}}=& 2(\beta -1){\displaystyle {\int }_0^{z_{\mathrm{s}}}\frac{\mathrm{d}z}{{\chi }^2(z)}\frac{\mathrm{d}{N}_{\mathrm{f}}}{\mathrm{d}z}{W}^{\mathrm{lens}}(z)}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }\frac{l\mathrm{d}l}{2\pi }{P}_{\mathrm{gal}-\mathrm{dm}}(l/{\chi }^2(z),z)J_0(l\theta ),}\hfill \end{array}$$(9)

where

$$\begin{array}{c}\hfill W^{\mathrm{lens}}(z)=\frac{3}{2}\frac{H_0^2}{c^2}{E}^2(z){\displaystyle {\int }_z^{z_{\mathrm{s}}}\mathrm{d}{z}^{\prime }\frac{\chi (z)\chi (z^{\prime }-z)}{\chi (z^{\prime })}\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{z}^{\prime }}}\end{array}$$(10)

with $E(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$, dN_b/dz, and dN_f/dz as the unit-normalised background and foreground redshift distribution, and z_s the source redshift. Because β is the logarithmic slope of the background sources number counts, N(S) = N₀S^−β, and χ(z) is the comoving distance to redshift z. It is common to assume β = 3 as that leads to the largest gravitational lensing effect in general, and the magnification bias in particular (Lapi et al. 2011, 2012; Cai et al. 2013; Bianchini et al. 2015, 2016; González-Nuevo et al. 2017; Bonavera et al. 2019). In the model β is a globally defined constant, but small variations in its value are almost completely compensated by changes in the M_min parameter. For example, an increase of 15% in β reduces log M_min in ∼1%.

The halo model formalism for the cross-correlation function implies that both the low and high redshift galaxy samples trace the same dark matter distribution. The foreground galaxies at the lenses redshift (z ∼ 0.4) directly trace the dark matter distribution, and the background galaxies (z >  1) indirectly trace it through the weak lensing effects. Following on from this consideration, the one-halo term relates to the subhalo correlation (by both samples) in the same halo and the two-halo term describes the correlation between one halo by the foreground galaxies and one halo by the background sources.

3.4. Estimation of parameters

To estimate the different sets of parameters, we performed a Markov chain Monte Carlo (MCMC) analysis using the open source emcee software package (Foreman-Mackey et al. 2013). It is an MIT licensed pure-Python implementation of Goodman & Weare (2010) Affine Invariant MCMC Ensemble sampler. In this analysis we generated at least 10 000 posterior samples to ensure a good statistical sampling after convergence.

The estimation of parameters in the foreground auto-correlation function analysis relies on the code¹ by Ronconi et al. (2020). We sampled the space of the model-defining parameters using MCMC to maximise a likelihood built as the sum of the χ² of the two measures we wanted to fit, namely the two-point angular correlation function at a given redshift, ω(θ, z), and the average number of sources at a given redshift, n_g(z):

$$\begin{array}{c}\hfill log\mathcal{L}\equiv -{\displaystyle \frac{1}{2}}({\chi }_{\omega (\theta ,z)}^2+{\chi }_{n_{\mathrm{g}}(z)}^2).\end{array}$$(11)

In the cross-correlation function analysis, we took into account both the astrophysical HOD parameters, and the cosmological parameters. The astrophysical parameters to be estimated are M_min, M₁, and α. The cosmological parameters we want to constrain are Ω_m, σ₈, and h. With the current samples, we do not have the statistical power to constrain Ω_B, Ω_Λ, and n_s in our analysis, and we keep them fixed to Planck cosmology, Ω_B = 0.0486, Ω_Λ = 0.6911, and n_s = 0.9667 (see Planck Collaboration VI 2020). A traditional Gaussian likelihood function was used in this case.

It should be noted that only the cross-correlation data in the weak lensing regime (θ ≥ 0.2 arcmin) are being taken into account for the fit since we are in the weak lensing approximation (see Bonavera et al. 2019, for a detailed discussion). Moreover, the foreground galaxy density n_g is not included in the cross-correlation analysis as an additional constraint.

3.5. Priors of parameters

The emcee package allows us to impose priors on the parameters to be estimated with MCMC, i.e. it allows MCMC to explore just the regions delimited by given boundaries. This is usually done in order to avoid spending computing time in searching physically unacceptable parameter space. In this work we want to study both a prior-free case (by imposing flat priors) and more constrained parameters regions (by applying Gaussian priors).

As for the auto-correlation method performed on the GAMA sample (see Sect. 3.2), we apply relatively extensive priors to the parameters in order to keep this analysis independent from the one with cross-correlation. In particular, the priors are set to log(M_min/M_⊙) and log(M₁/M_⊙) between 9 and 15, α between 0.1 and 1.5, and A between −5 and −1.

As for the cross-correlation method (see Sect. 3.3), we set the priors of the astrophysical parameters to min-max values (in the flat priors case) and mean-σ values (in the Gaussian priors case) in agreement with information available in the literature. In particular, after a review of the available literature on the foreground sample characteristics we extract the following flat (FPA) and Gaussian (GPA) priors for our astrophysical parameters.

As described above, the astrophysical parameters involved in our estimations are M_min, M₁, and α. In particular, we rely on Sifón et al. (2015) for M_min, on both Viola et al. (2015) and Sifón et al. (2015) for M₁, and on Viola et al. (2015) for α.

As for M_min, our choice of Sifón et al. (2015) is based on the fact that their sample is similar to ours: it consists of the 100 deg² overlapping region between KiDS and GAMA surveys, with a redshift range of 0 − 0.5 and an average of 0.25. In their work the average galaxy halo mass is determined for a sample of satellite galaxies in massive galaxy groups (M >  10¹³ h⁻¹ M_⊙). In particular, in their Fig. 2 (top right panel) the satellite distributions of stellar mass to the BCG are shown. Since M_⋆ are the satellite masses, they are a proxy of the minimum mass (galaxies with no satellites), but it has to be expressed in terms of halo mass in order to be considered a limit on M_min. We translate these stellar masses in the corresponding mass of the halo using the recipe by Pantoni et al. (2019) to switch from log M_⋆(h⁻¹ M_⊙) to log M_h. Sifón et al. (2015) gives values of log M_⋆(h⁻¹ M_⊙) between 9.5 and 11.5, which correspond to the range log M_h = 11.6 − 13.6. So, we use these values to set our flat priors on log M_min.

In order to apply a more constraining range on the allowed values, we also investigate the case with Gaussian priors by relying again on the work by Sifón et al. (2015). In particular, we focus on their Table 1; it gives values of log⟨M_⋆⟩(h⁻¹ M_⊙) between 10.36 and 10.81, the lowest possible value being in Bin 1 (taking into account the lower error) and the highest possible value in Bin 3 (taking into account the upper error). Such stellar masses translate (using Pantoni et al. 2019) into log M_h(M_⊙) between 12.3 and 12.5. So, we choose a mean value of μ = 12.4 and a dispersion of σ = 0.1 for the Gaussian priors. In this way the log M_min priors stay in the Sifón et al. (2015) more restrictive findings within the 68% confidence level (CL).

It should be noted that these flat and Gaussian priors are in agreement with Aversa et al. (2015). In particular, their Fig. 28 shows values of log M_h(M_⊙) between 11.5 and 13.5, and in their Fig. 5 they gave a value of log M_⋆(M_⊙) = 10.5 that translates to log M_h(M_⊙) = 12.3.

Concerning M₁, we rely on Viola et al. (2015; see their Fig. 16) to set the wide flat priors to log M₁(M_⊙) = 13 − 14.5. In order to limit this parameter to a more restricted range, we consider again Sifón et al. (2015) instead. The Gaussian priors are obtained using the information in their Table 2: we choose the range log M_host(h⁻¹ M_⊙) = 13.51 − 14.8, which spans from the lowest possible value (i.e. the one for Bin 1 including the lower error) and the highest possible value (i.e. the one for Bin 3 taking into account the upper error). By writing in the units we are using in this work, we get a typical halo mass with at least one satellite in the range log M₁(M_⊙) = 13.66 − 14.33. To set the Gaussian priors based on these values, we take the mean μ = 13.95 of this range and a dispersion of σ = 0.3, so that the Sifón et al. (2015) value corresponds to the 68% CL of our priors. We can argue that even smaller Gaussian priors can be set by only relying on their results for Bin 1 (in Sifón et al. 2015, Table 2). The choice of this bin is justified because we are looking for the minimum mass to have at least one satellite. However, we decided to use a less constraining prior in order to be more conservative.

With respect to α, we derive our priors by comparing our Eq. (2) to the expression for α in Viola et al. (2015) (Eq. (40)). From this comparison, we obtain a value of α = 0.92 ± 0.15. So, we set the Gaussian priors with mean μ = 0.92 and dispersion σ = 0.15. To define the flat priors, we consider instead a 3σ range around the mean value, obtaining α = 0.5 − 1.37. We also perform an additional analysis relying on the auto-correlation results, by setting the Gaussian priors to the μ and σ rounded values in Table 2. Finally, we set common flat priors for the Cosmological parameters: 0.1−0.8 for Ω_m, 0.6−1.2 for σ₈, and 0.5−1.0 for h.

4. Results

To gather as much information as possible about the derived cosmological constraints and their robustness, first we fix the cosmology to the Planck one and study the behaviour of the astrophysical parameters. We decide to use a non-informative, or uniform/flat, priors for our free astrophysical parameters and we perform the analysis for both data sets (the default and the shifted binning). Then we also allow the cosmological parameters to vary using flat priors for three different astrophysical parameter assumptions: (a) flat priors from the literature, (b) Gaussian priors from the literature, and (c) Gaussian priors from the auto-correlation analysis described in Sect. 3.2.

4.1. Astrophysical parameters

4.1.1. Auto-correlation data analysis

The posterior distributions for the model parameters (see Sect. 3.2.2) derived from the estimated auto-correlation function (see Sect. 3.2) through the MCMC analysis (see Sect. 3.4) are shown in Fig. 4². For this analysis we fix the cosmology to the Planck values: Ω_m = 0.3089, σ₈ = 0.8159, and h = 0.6774.

Fig. 4. Posterior distributions for the model parameters in the auto-correlation data analysis. The contours for these plots are set to 0.393 and 0.865 (see text for details).

The mean values at 68% CL, standard deviations, and peak values are summarised in Table 2. The auto-correlation function calculated using the best-fit parameter values is also shown in Fig. 2. The model parameters are well constrained and they provide a good fit to the observed auto-correlation function, obtaining a ${\chi }_{\omega }^2$ of 8.836 for 11 degrees of freedom (d.o.f.; corresponding to a reduced ${\chi }_{\omega }^2$ of 0.803) and a ${\chi }_{n_{\text{g}}}^2$ of 0.001 (0.0001). The value for the ${\chi }_{n_{\text{g}}}^2$ is so small because we estimated just one single value and this helps in its correct recovery.

Moreover, the values obtained for log(M_min/M_⊙) and log(M₁/M_⊙) ($12.{362}_{-0.007}^{+0.002}$ and $13.{64}_{-0.09}^{+0.03}$, respectively) are consistent with what can be found in the literature, as discussed in Sect. 4.1 (log(M_min/M_⊙) = 11.6 − 13.6 and log(M₁/M_⊙) = 13.0 − 14.5 for the flat priors cases). The derived α parameter ($1.{17}_{-0.05}^{+0.06}$) is also in agreement with the literature values (0.5 − 1.37 for the flat priors case), although slightly toward the upper limit.

Therefore, we can conclude that the auto-correlation function analysis results can be an additional alternative set of Gaussian parameters for the astrophysical parameters to be used in the full cosmological analysis.

4.1.2. Cross-correlation data analysis

We perform a similar analysis as in the auto-correlation case for the cross-correlation function. We also fix the cosmology to the Planck values, and perform the MCMC analysis on the astrophysical parameters only.

The mean values, CL, standard deviations, and the peak values are summarised in Table 3. The results of the MCMC runs for the two shifted data sets are listed in the last two columns. The posterior distributions shown in Fig. 5 give compatible results to the auto-correlation case (summarised in Table 2), in particular for M_min, the only one that is well constrained.

Fig. 5. Results obtained from the cross-correlation data with flat priors and fixed cosmology on the cross-correlation data in the two data sets described in Sect. 2 to obtain constraints on the astrophysical parameters. Priors and results are listed in Table 3. The contours for these plots are set to 0.393 and 0.865 (see text for details).

The value of M_min is well constrained with both data sets (default and shifted binning) with almost identical results ($log(M_{\min }/M_{\odot })=12.{42}_{-0.11}^{+0.20}$ and $12.{41}_{-0.07}^{+0.17}$, respectively) and in agreement with the results from González-Nuevo et al. (2017) and Bonavera et al. (2019), and to typical values in the literature, as described in Sect. 4.1.

It is worth mentioning here that matter distribution at z <  0.2 is not significantly biasing our lens sample (because of lensing, see Sect. 2.2). This is confirmed because the mass we estimate with the cross-correlation for our selected GAMA sample is in fairly good agreement with that found by the auto-correlation in the GAMA sample (Sect. 4.1.1), where only very few sources might be affected by lensing and the mass determination is dominated by the unlensed sources. Since the cross- and auto-correlation results are in agreement, it means that cross-correlation results are also not affected by magnification due to z <  0.2 matter distribution.

On the contrary, the M₁ and α parameters are both unconstrained with the two data sets. The results are mostly the mean values between the flat prior limits as expected in this situation. However, there is a hint in the shifted data set toward a similar log(M₁/M_⊙) value (with a peak value of 13.64), but lower α value (∼0.9) with respect to the auto-correlation case.

Figure A.1 is of great help in understanding this question. Once the M_min is fixed, for the flat ranges adopted for both M₁ and α, the variation is small enough to avoid a clear constraint, even more if we consider that both parameters are probably degenerate to a certain degree. Moreover, focusing on M₁, the analysis based on the shifted cross-correlation function is more constraining than the other one (see e.g. Fig. 6). What is causing the difference between these two cases might probably be the first cross-correlation point at θ ≳ 0.2 arcmin which appears to be lower in the shifted case than in the other. Looking again at Fig. A.1, it seems that lower values of M₁ predict lower amplitude of the cross-correlation function at those scales.

Fig. 6. Results obtained from the cross-correlation data with flat priors on astrophysical and cosmological parameters, in the two data set for the cross-correlation data. The limits on the priors for the astrophysical and cosmological parameters and the corresponding results are listed in Table 4. The contours for these plots are set to 0.393 and 0.865 (see text for details).

For these reasons, we decided to consider these results simply as a self-consistency check of the astrophysical parameters outputs. For the full cosmological analysis, we restricted ourselves to the priors discussed before and the auto-correlation results.

Finally, taking into account the results, we can also conclude that they are robust against small variations in the data set used.

4.2. Cosmological parameters

Having shown that using flat priors on astrophysical parameters provides values compatible with the literature, we extend the cross-correlation function analysis to allow for variations in the main cosmological parameters, Ω_m, σ₈, and H₀. We would like to note at this point that the main aim of this work is to exploit the cross-correlation measurements in order to constrain these parameters using flat priors on both astrophysical and cosmological parameters (FPA FPC). So we first present our results for this more general case, and then investigate the effects and possible improvements on estimating our cosmological parameters when setting the astrophysical parameters to the Gaussian priors derived from the literature and from auto-correlation results.

The best-fit posterior distributions of the free parameters is shown in Fig. 6 (red for the default binning and blue for the shifted). The best-fit values are listed in Table 4. The solid lines in Fig. 3 represent the halo model best fit, estimated using the MCMC approach. In particular, the red is for the main data set, which gives a χ² of 8.16 for 6 d.o.f. (corresponding to a reduced χ² of 1.36; the d.o.f. are the 12 data points minus the 6 parameters) and the blue is for the shifted data set, with a χ² of 9.92 (1.65, when taking into account the d.o.f.). The green line corresponds to the parameters obtained running on the main data set except for Ω_m that is set equal to 0.3. In this case, we get χ² of 13.76 (2.29).

First, the recovered astrophysical parameters are compatible with the results of the previous section (see Table 3). This also means that only M_min is well constrained. Therefore, the conclusion that the results are robust against small variations in the data sets can also be extended to the cosmological case.

With respect to the cosmological parameter results, they are not well constrained in this proof of concept, but there are interesting conclusions that can be drawn. The H₀ parameter is unconstrained, but for σ₈ there is a clear broad peak around 0.75 and its posterior distribution indicates an upper limit of < 1.0 at the 95% CL. The analysis of Ω_m is more complicated, but interesting at the same time. As clearly shown by Fig. A.1, Ω_m is most sensitive to the cross-correlation measurements on the largest scales, where the uncertainties are largest (the three last observations are almost upper limits). This translates into a wide posterior distribution toward higher Ω_m values driven mostly by the lower value of the last point. On the contrary, there is also a strong constraint on the lower Ω_m value. This analysis implies that Ω_m ≳ 0.24 at the 95% CL.

In order to investigate possible improvements on the cosmological constraints, we perform additional analyses with Gaussian priors on the astrophysical parameters and flat priors on the cosmological parameters.

The first choice for the Gaussian priors is the values from the literature described in Sect. 4.1. The priors for these runs are listed in the second column of Table 5. The remaining columns list the results of the MCMC for both data sets (same as Table 4). The best-fit posterior distributions of the free parameters are shown in Fig. 7 (red for the default binning and blue for the shifted).

Fig. 7. Results obtained from the cross-correlation data with Gaussian astrophysical and flat cosmological priors in the two shifted data sets for the cross-correlation data. The results are listed in Table 5. The Gaussian priors are set to the μ and σ values presented by Viola et al. (2015), Sifón et al. (2015), and Pantoni et al. (2019). The contours for these plots are set to 0.393 and 0.865 (see text for details).

We also analyse the case with Gaussian priors on the astrophysical parameters fixed with the auto-correlation analysis results described in Sect. 3.2. Our findings in this case are summarised in Table 6 and Fig. 8.

Fig. 8. Results obtained from the cross-correlation data with Gaussian astrophysical priors from auto-correlation analysis and flat cosmological priors in the two shifted data sets for the cross-correlation data. The results are listed in Table 6. The contours for these plots are set to 0.393 and 0.865 (see text for details).

As expected, the astrophysical parameters are prior dominated. In addition, by imposing Gaussian priors to the astrophysical parameters the uncertainties in the cosmological parameters are reduced. This is clearly shown by comparing the results using the more restrictive priors from the auto-correlation analysis with the Gaussian priors from the literature. For example, the σ₈ parameter improves from a ${\sigma }_8=0.{84}_{-0.10}^{+0.10}$ (Gaussian prior from the literature) to ${\sigma }_8=0.{74}_{-0.03}^{+0.03}$ (Gaussian prior from the auto-correlation analysis). In both cases there is also a slight correlation between σ₈ and M₁ parameters as both affect the cross-correlation function in an opposite way on intermediate angular scales (see Fig. A.1).

Both Ω_m and H₀ remain almost unconstrained, as in the flat prior case. There is slight tendency of H₀ toward values above 70 km s Mpc⁻¹, as is discussed in more detail in the next section. The strong lower limit constraint found with the flat priors for Ω_m is confirmed by additional analysis. Moreover, using the auto-correlation Gaussian priors, the posterior distribution of Ω_m shows a peak around 0.41 in the shifted case. This is an indication that reducing the uncertainties in the observed cross-correlation function and accurate astrophysical priors can provide strong independent constraints based on the magnification bias.

4.3. Comparison with other results

In this section we compare our findings with the cosmological constraints obtained using a different observable: the shear of the weak gravitational lensing. We consider both the CMB lensing from Planck (Planck Collaboration VIII 2020) and the cosmic shear of galaxies as measured by several surveys, focusing on those that are the most constraining to date. Specifically, we include results from the revised analysis of the cosmic shear tomography measurements of the Canada-France-Hawaii Telescope Lensing Survey presented in CFHTLenS (Joudaki et al. 2017), the first combined cosmological measurements of the Kilo Degree Survey and VIKING based on 450 deg² data (KV450, Hildebrandt et al. 2020), the first-year lensing data from the Dark Energy Survey (DES, Troxel et al. (2018)), together with recent constraints from the two-point correlation functions of the Subaru Hyper Suprime-Cam first-year data (HSC, Hamana et al. 2020). It is worth noting that these results span somewhat different redshift ranges and are affected by completely different systematic effects with respect to the magnification bias measurements analysed in this paper, which thus represent a complementary probe to weak lensing.

We focus on the comparison to the constraints in the Ω_m–σ₈ plane since weak lensing is most sensitive to these cosmological parameters. In particular, cosmic shear approximately constrains the combination ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.5}$, while ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$ is what is measured by CMB lensing. These combinations identify degeneracy directions that are clearly visible in Figs. 9 and 10, where we show the marginalised posterior contours (68% and 95% CL) for the above-mentioned data sets. We use MCMC results that have been publicly released by the different collaborations. It should be noted that in these comparative figures, in order to be consistent with the results shown in the literature for the experiments mentioned above, we set the contours to 0.68 and 0.95, which are different from the values of 0.393 and 0.865 used in the corner plots shown in the previous sections, and that actually corresponds to the relevant 1σ and 2σ levels for a 2D histogram of samples. The choice of 0.68 and 0.95 contours in this subsection is for a direct comparison with the literature; the 0.393 and 0.865 selection in the 2D histograms of the corner plots is for a direct comparison with the 1σ and 2σ levels in the 1D histograms in the upper part of the same corner plots.

Fig. 9. Ωm–σ8 plot comparing our results using flat priors in the astrophysical and cosmological parameters (both cases, grey and black) with results from Planck (blue and dark blue for the lensing only and all Planck cases, respectively), CFHTLenS (red), KV450 (orange), DES lensing (green), and HSC (cyan). The contours for these plots are set to 0.68 and 0.95.

Fig. 10. Ωm–σ8 plot comparing our results using Gaussian priors (left panel: from the literature; right panel: from the auto-correlation analysis) in the cosmological parameters (both cases, grey and black) with results from Planck (blue and dark blue for the lensing only and all Planck cases, respectively), CFHTLenS (red), KV450 (orange), DES lensing (green), and HSC (cyan). The contours for these plots are set to 0.68 and 0.95.

It should be noted that different studies adopt different cosmological parameters priors; however, we do not attempt to adjust them to our fiducial set-up as an in-depth comparison is beyond the scope of this paper. For completeness, in the figures we also show results from Planck CMB temperature and polarisation angular power spectra (dark blue). They have been known to be in tension with the CFHTLenS (red) and KV450 (orange) data, while a better agreement has instead been found with the displayed HSC (cyan) and DES (green) constraints. For the last, in particular, Troxel et al. (2018) made the conservative choice of excluding the non-linear scales where modelling uncertainties are higher. This also resulted in a reduction of the constraining power of the data set.

The cosmological constraints derived in this paper are shown for both of the data sets (grey filled contours and black dashed curves). The main results, based on the flat prior for both the astrophysical and cosmological parameters, are shown in Fig. 9. In Fig. 10, we show how the results can be improved using Gaussian priors on the astrophysical parameters. Even if our data sets might not seem constraining, two important conclusions can be drawn from the results. Firstly, the estimates provide a lower bound Ω_m >  0.24 at 95% CL that is robust against the choice of priors. Secondly, our results do not show the typical degeneracy that characterises the cosmic shear results. This means that future more constraining magnification bias samples can offer a complementary probe able to break the degeneracy (see e.g. Fig. 9).

The results on Ω_m can be better understood looking at Fig. A.1, which shows how the cross-correlation function depends on cosmological parameters. It is interesting to see how the largest angular separations are very sensitive to the matter density parameter, and how present data disfavour low values of Ω_m.

We note that the data in this paper alone cannot place useful constraints on the Hubble constant, as is evident from Fig. 11. This is usually the case even for weak lensing measurements. However, in most of our cases there is a mild pull towards higher H₀ values.

Fig. 11. Posterior distributions on H0 derived using the different data sets and assuming different priors for the astrophysical parameters. As a comparison the result by Riess et al. (2019; grey band) and Planck are also shown.

As a comparison, we also add the vertical grey band, corresponding to the Riess et al. (2019) results and the Planck constraints. However, due to our large error bars on the cross-correlation data, our results are not competitive with those in the recent literature and cannot give any hint on the Planck–Riess et al. (2019) tension problem yet.

4.4. Future developments

This work is meant to show that the magnification bias can become an additional independent cosmological probe to help to resolve the current tensions in modern Cosmology. Results are not yet competitive with current estimation of cosmological parameters, but some improvements are expected if smaller error-bars on the cross-correlation measurements can be achieved, especially on larger scales where the Ω_m and H₀ parameters mostly affect the cross-correlation function.

First of all, it should be noted that this work is based on just one (relatively large) redshift bin. As shown in González-Nuevo et al. (2017), a tomographic analysis can be done with such cross-correlation measurements. When we divide the sample into redshift bins, we reduce the sample size and we expect the error bars to increase. This clear disadvantage might be counterbalanced by easier constraining of the cosmological parameters, since only sources with very similar redshift and the known evolution of the cosmological parameters are being taken into consideration. This is a possibility that needs to be further explored.

Regarding the reduction of the error bars, the statistics can be increased enlarging both the lenses and the background source samples. As for the lenses, the sample might be increased with current data by considering sources with photometric redshifts. On the one hand, the number of sources with photometric redshifts is larger than that of sources with the spectroscopic redshifts. On the other hand, photometric redshifts are less accurate, and in this case the application of the cross-correlation that needs very clearly redshift-separated samples of foreground and background sources might not be straightforward. In the future more precise catalogues of lens sources might be available, such as from DES (Dark Energy Survey Collaboration 2016), the Javalambre Physics of the Accelerating Universe Astrophysical Survey (J-PAS, Benitez et al. 2014) and the Euclid (Laureijs et al. 2011) experiments/missions.

For the background sources, the number of SMGs might be increased by considering the catalogue of the whole area covered by Herschel (e.g. by investigating the possibility of using the Herschel Extragalactic Legacy Project, HELP; Shirley et al. 2019), and not just H-ATLAS data. Moreover, it will surely be increased by catalogues of new instruments in the near-infrared, such as the James Webb Space Telescope (JWST; Gardner et al. 2006).

5. Conclusions

The aim of this work was to test the capability of magnification bias with high-z SMGs as a cosmological probe. As demonstrated in González-Nuevo et al. (2017) and Bonavera et al. (2019), magnification bias is a powerful tool for constraining the free parameters of the HOD model. We now exploit cross-correlation data to constrain not only astrophysical parameters (M_min, M₁, and α), but also some of the cosmological ones. In particular, we realised that at this early stage the n_s, Ω_B, and Ω_Λ parameters cannot be constrained. We thus focus our analysis on Ω_m, σ₈, and H₀ for this proof of concept.

With our current results Ω_m and H₀ cannot be well constrained, although we showed that magnification bias can be a way to set constraints on these parameters, but we conclude that the most important issue regarding the Ω_m and H₀ parameters is the uncertainty on the cross-correlation data on large scales. So our future efforts must be focused on reducing the error bars on the data set. What we need is a larger area for the background sources to get a wider area and achieve a better statistics.

Moreover, we do not find a degeneracy between Ω_m and σ₈. If any, the correlation might be opposite to the common one (see Fig. 10): lower σ₈ implies in our case lower Ω_m in order to compensate the effect on large scales.

With all the cases studied in this work we can set a lower limit of 0.24 at the 95% CL on Ω_m, and we see a mild preference towards high H₀ values. Since Ω_m was allowed to vary in the range [0.1−0.8], this lower bound is a very robust conclusion and its importance is clearly shown in the comparison with the results from weak lensing data: our lower bound still imposes restrictions to the other weak lensing results.

For our constraints on σ₈, we obtained only a tentative peak around 0.75 but an interesting upper limit of σ₈ ≲ 1 at 95% CL. By imposing Gaussian priors to the astrophysical parameters, the constraints are improved: ${\sigma }_8=0.{74}_{-0.03}^{+0.03}$ using the Gaussian prior from the auto-correlation analysis.

Finally, our estimate of the H₀ parameter is still not accurate, being inclusive of the Planck and the Riess et al. (2019) results. Better results might be obtained if better measurements and smaller error bars for the cross-correlation points on large scales can be achieved.

Acknowledgments

The authors greatly thanks the anonymous referee for providing precious and appropriated comments. LB, JGN, MMC acknowledge the PGC 2018 project PGC2018-101948-B-I00 (MICINN/FEDER). LB and JGN acknowledge PAPI-19-EMERG-11 (Universidad de Oviedo). MM is supported by the program for young researchers “Rita Levi Montalcini” year 2015. LD and SJM acknowledge support from the European Research Council Advanced Investigator grant, COSMICISM and Consolidator grant, COSMIC DUST. AL acknowledges partial support from PRIN MIUR 2015 Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid, by the RADIOFOREGROUNDS grant (COMPET-05-2015, agreement number 687312) of the European Union Horizon 2020 research and innovation program, and the MIUR grant “Finanziamento annuale individuale attivita base di ricerca”. Finally, MN received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 707601. We deeply acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. In particular the projects “SIS19_lapi”, “SIS20_lapi” in the framework “Convenzione triennale SISSA-CINECA”. 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. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programs including GALEX MIS, VST KIDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA web-site is: http://www.gama-survey.org. In this work, we made extensive use of GetDist (Lewis 2019), a Python package for analysing and plotting MC samples. In addition, this research has made use of the python packages ipython (Pérez & Granger 2007), matplotlib (Hunter 2007) and Scipy (Jones et al. 2001).

References

Appendix A: Cross-correlation function sensitivity to the model’s parameters

Figure A.1 shows the dependence of the cross-correlation function with each of the model’s parameters: M_min, M₁, α, Ω_m, σ₈, and H₀ (from left to right and top to bottom).

Fig. A.1. Dependency of the cross-correlation function with each of the model’s parameters (maintaining the others fixed).

Of the astrophysical parameters, the most sensitive one is M_min. Considering the value ranges used in this work for M₁ and α, it is clear that the variation is located only at intermediate angular scales and it is small compared with the observational uncertainties. This explains that both parameters are almost unconstrained with the cross-correlation data alone.

With respect to the cosmological parameters, the most sensitive one is σ₈. Its effect on the cross-correlation function are most important on intermediate angular scales, but it also has a non-negligible effect on large angular scales. On the contrary, the effect of Ω_m is only important on large angular scales and negligible on lower angular scales. Moreover, higher values of σ₈ increase the cross-correlation signal, while the opposite occurs for Ω_m. This explains how the usual Ω_m–σ₈ degeneracy could be broken with the cross-correlation function once more accurate observation at large angular scales are available.

Finally, the effect produced by H₀ is relevant only on very large angular scales (even larger than in the Ω_m case). As currently these angular scales are those observed with the largest uncertainties, H₀ is unconstrained.

All Tables

All Figures

	Fig. 1. Redshift distributions of the background H-ATLAS sample (red) and the foreground GAMA sample (blue).
In the text

	Fig. 2. Measured angular auto-correlation function of the foreground sample (red circles). The best-fit halo model prediction is shown as a comparison (total, blue solid line; one-halo term, red dotted line; two-halo term, green dashed line).
In the text

Fig. 3. Cross-correlation data computed on Herschel SMGs (as background sample) and GAMA galaxies (as lenses sample). The grey points are exactly the results obtained in González-Nuevo et al. (2017). The red data set (obtained using the same binning as in González-Nuevo et al. 2017) and the blue data set (obtained using same the binning, but half-size bin shifted with respect to González-Nuevo et al. 2017) are the cross-correlation results computed for this work, exploiting a larger area. The solid lines represent the best fit with flat priors to the astrophysical and cosmological parameters using the default binning (red) and the shifted binning (blue). The green line is the best fit to the red case, but adopting Ωm = 0.3.

In the text

	Fig. 4. Posterior distributions for the model parameters in the auto-correlation data analysis. The contours for these plots are set to 0.393 and 0.865 (see text for details).
In the text

	Fig. 5. Results obtained from the cross-correlation data with flat priors and fixed cosmology on the cross-correlation data in the two data sets described in Sect. 2 to obtain constraints on the astrophysical parameters. Priors and results are listed in Table 3. The contours for these plots are set to 0.393 and 0.865 (see text for details).
In the text

	Fig. 6. Results obtained from the cross-correlation data with flat priors on astrophysical and cosmological parameters, in the two data set for the cross-correlation data. The limits on the priors for the astrophysical and cosmological parameters and the corresponding results are listed in Table 4. The contours for these plots are set to 0.393 and 0.865 (see text for details).
In the text

Fig. 7. Results obtained from the cross-correlation data with Gaussian astrophysical and flat cosmological priors in the two shifted data sets for the cross-correlation data. The results are listed in Table 5. The Gaussian priors are set to the μ and σ values presented by Viola et al. (2015), Sifón et al. (2015), and Pantoni et al. (2019). The contours for these plots are set to 0.393 and 0.865 (see text for details).

In the text

	Fig. 8. Results obtained from the cross-correlation data with Gaussian astrophysical priors from auto-correlation analysis and flat cosmological priors in the two shifted data sets for the cross-correlation data. The results are listed in Table 6. The contours for these plots are set to 0.393 and 0.865 (see text for details).
In the text

	Fig. 9. Ωm–σ8 plot comparing our results using flat priors in the astrophysical and cosmological parameters (both cases, grey and black) with results from Planck (blue and dark blue for the lensing only and all Planck cases, respectively), CFHTLenS (red), KV450 (orange), DES lensing (green), and HSC (cyan). The contours for these plots are set to 0.68 and 0.95.
In the text

Fig. 10. Ωm–σ8 plot comparing our results using Gaussian priors (left panel: from the literature; right panel: from the auto-correlation analysis) in the cosmological parameters (both cases, grey and black) with results from Planck (blue and dark blue for the lensing only and all Planck cases, respectively), CFHTLenS (red), KV450 (orange), DES lensing (green), and HSC (cyan). The contours for these plots are set to 0.68 and 0.95.

In the text

	Fig. 11. Posterior distributions on H0 derived using the different data sets and assuming different priors for the astrophysical parameters. As a comparison the result by Riess et al. (2019; grey band) and Planck are also shown.
In the text

	Fig. A.1. Dependency of the cross-correlation function with each of the model’s parameters (maintaining the others fixed).
In the text