© The Authors 2024

1. Introduction

The excess of high-redshift sources near low-redshift mass structures, known as magnification bias, has the potential to serve as a valuable tool for investigating the large-scale structure of the universe and the distribution of dark matter when applied to background submillimeter galaxies (Bonavera et al. 2022, and references therein). This can be done by providing independent constraints on cosmological parameters, whereby measurements of the magnification bias can complement other cosmological probes. The phenomenon caused by gravitational lensing, namely, the distribution of matter in the Universe (including galaxies and clusters of galaxies) bending light as it travels through space, which results in the apparent brightness and size of distant objects being magnified or distorted (e.g. Schneider et al. 1992). The magnification bias effect manifests as a significant cross-correlation function between two source samples with different redshift ranges, indicating that the large-scale structure traced by the foreground sources is amplifying the background ones. This phenomenon has been observed in various situations, such as the correlation between galaxies and quasars (e.g. Scranton et al. 2005; Ménard et al. 2010), Herschel sources and Lyman-break galaxies (Hildebrandt et al. 2013), and between the cosmic microwave background (CMB) and other sources (Bianchini et al. 2015, 2016). In particular, submillimeter galaxies (SMGs) are very suitable as background sources due to their exceptional properties, which include a steep luminosity function, high redshifts, and faint optical band emission (see e.g. González-Nuevo et al. 2012, 2017; Blain 1996; Negrello et al. 2007, 2010, 2017; Bussmann et al. 2012, 2013; Fu et al. 2012; Wardlow et al. 2013; Calanog et al. 2014; Nayyeri et al. 2016; Bakx et al. 2020). Previous studies have demonstrated the magnification bias effect on SMGs and have measured it to derive both astrophysical and cosmological information (González-Nuevo et al. 2017, 2021; Bonavera et al. 2020; Cueli et al. 2021). Additionally, the ability to divide the foreground sample into distinct redshift bins allows for a more detailed tomographic analysis (Bonavera et al. 2021; Cueli et al. 2022).

Indeed, similarly to shear studies (Bacon et al. 2000; Van Waerbeke et al. 2000; Wittman et al. 2000; Rhodes et al. 2001), cosmological parameters could benefit from a tomographic analysis, especially by allowing for better estimation of the halo occupation distribution (HOD) parameters. Slight variations in these parameters are expected when explored in distinct redshift bins. González-Nuevo et al. (2017) demonstrated this by performing tomographic studies on the HOD parameters, separating the foreground sample into four redshift bins: 0.1 − 0.2, 0.2 − 0.3, 0.3 − 0.5, and 0.5 − 0.8. It was discovered that the values of M_min, the mean minimum mass for a halo to host a (central) galaxy, exhibited a clear evolution, increasing with redshift as predicted by theoretical calculations. Furthermore, in González-Nuevo et al. (2021), a simplified tomographic analysis is conducted using a background sample of Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS) galaxies with photometric redshifts of z > 1.2, and two independent foreground samples: Galaxy And Mass Assembly (GAMA) survey galaxies with spectroscopic redshifts and Sloan Digital Sky Survey (SDSS) galaxies with photometric redshifts, with 0.2 < z < 0.8. The study obtained constraints on the cosmological parameters Ω_M ($0.{50}_{-0.20}^{+0.14}$) and σ₈ ($0.{75}_{-0.10}^{+0.07}$) as a 68% confidence interval (C.I.).

Bonavera et al. (2021) investigated the possibility of improving constraints on cosmological parameters Ω_M, σ₈, and H₀ by using unbiased measurements of the H-ATLAS/GAMA cross-correlation function, as concluded by the methodological analysis of González-Nuevo et al. (2021) and a tomographic analysis. These authors found a trend towards higher values of M_min at higher redshifts. They also explored different values of the dark energy equation of state parameter, ω. Their results showed maximum posterior values, at a 68% C.I., of $0.{26}_{-0.09}^{+0.15}$ for Ω_M, and $0.{87}_{-0.12}^{+0.13}$ for σ₈ in the ΛCDM model. The study also examined a more general ω₀CDM model, allowing for values of ω₀ ≠ −1, and found that the constraints on Ω_M and σ₈ were similar. However, the maximum posterior value for ω₀ was found to be at $-1.{00}_{-0.56}^{+0.53}$, with a 68% C.I. In the ω₀ω_aCDM model, where a possible z dependence of the ω parameter is allowed, the results showed a maximum posterior value of $-1.{09}_{-0.63}^{+0.43}$ with a 68% C.I. for ω₀ and a maximum posterior value of $-0.{19}_{-1.69}^{+1.29}$ with a 68% C.I. for ω_a. The tomographic analysis improved the constraints on the σ₈ − Ω_M plane and demonstrated that the magnification bias does not to show the degeneracy found with cosmic shear measurements. Finally, the study revealed a trend of higher ω₀ values for lower H₀ values concerning dark energy.

This work is the companion paper to Cueli et al. (2023), also referred as Paper I, that focuses on deriving cosmological and astrophysical constraints using a single broad foreground redshift bin. Additionally, the theoretical modelling of the signal is reexamined with respect to the work of Bonavera et al. (2021) and González-Nuevo et al. (2021) to assess the significance of the logarithmic slope of background galaxy number counts and a specific numerical correction.

The goal of this second work is to assess the performance of this observable in a tomographic setup after the methodological refinement described in Paper I. Given the numerous parameters to be jointly analysed in tomography and the numerical complexity of the theoretical model, the computational time required for this task is large. The aim is to investigate the possibility to optimise CPU time and strategies, such as selecting different redshift bins, to improve the stability and accuracy of the results. This work suggests various strategies to streamline CPU time without compromising the precision of the measurements. It also delves into the selection of different redshift bins and their effect on the ultimate outcome.

In this manuscript, we have structured the content as follows. Section 2 summarises the methodology of this work (data, theoretical framework, and parameter estimation). The results obtained for different combinations of redshift bins and ranges are summarised in Sect. 3. Finally, in Sect. 4, we draw conclusions based on our findings.

2. Methodology

2.1. Data and redshift subsample

The details of the data used for computing the cross-correlation signal are discussed in detail in Paper I (Cueli et al. 2023). The background sample primarily consists of sources with photometric redshifts, while the foreground sample comprises sources with spectroscopic redshifts. The analysis involves dividing the foreground sample into redshift bins, necessitating the use of spectroscopic measurements for this particular sample in order to avoid overlapping between bins. Photometric errors for the background sample are considered to smooth the background redshift distribution required for the theoretical model but do not present any concerns in this context.

Specifically, the background sample is selected from the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) sources in the Galaxy And Mass Assembly (GAMA) fields (G09, G12 and G15) which covers approximately ∼147 deg² of common area and part of the south Galactic pole (SGP), totalling around 60 deg². A photometric redshift selection of 1.2 < z < 4.0 was applied to avoid redshift overlap with the lenses. As a result, we obtained a sample of approximately 66 000 sources with an average photometric redshift of 2.20. The redshift distribution of the background sources used in the study is estimated as p(z|W) for galaxies selected by a window function with a redshift range of 1.2 < z < 4.0. The distribution takes into account the effect of random errors in photometric redshifts, following the methodology described in González-Nuevo et al. (2017).

The foreground or lens sample used in this work is referred to as the “z_spec sample” in González-Nuevo et al. (2021), consisting of approximately 150 000 galaxies with spectroscopic redshifts in the range of 0.2 < z_spec < 0.8. These redshift values were obtained from the GAMA II spectroscopic survey. The total common area used in this study was approximately 207 deg². To facilitate our analysis, we divided the foreground sample into different sub-samples based on their redshift. Following the previous results from Bonavera et al. (2021), we adopted the “base case” which involved dividing the entire redshift range into four bins: 0.1 < z < 0.2 (bin 1), 0.2 < z < 0.3 (bin 2), 0.3 < z < 0.5 (bin 3), and 0.5 < z < 0.8 (bin 4). Since the spectroscopic redshifts for the foreground sources have negligible errors, there is no overlap in the redshift distribution among bins, and no source is misplaced in the wrong bin, although this detail is not critical for cross-correlation analyses with a non-overlapping background sample.

To investigate the possible dependence on the number of redshift bins but maintaining the same redshfit range, we conducted a joint analysis using six redshift bins, referred to as the “six bins” scenario. In this scenario, the entire range was divided into 0.1 < z < 0.2 (bin 1), 0.2 < z < 0.25 (bin 2a), 0.25 < z < 0.3 (bin 2b), 0.3 < z < 0.4 (bin 3a), 0.4 < z < 0.5 (bin 3b), and 0.5 < z < 0.8 (bin 4). Finally, we also conducted a joint analysis using six bins of redshift, but with a wider range, denoted as the “six bins-WR” case. In this case, the range extended from 0.01 to 0.9, with the following bin divisions: 0.01 < z < 0.1 (bin 0), 0.1 < z < 0.2 (bin 1), 0.2 < z < 0.3 (bin 2), 0.3 < z < 0.5 (bin 3), 0.5 < z < 0.7 (bin 4), and 0.7 < z < 0.9 (bin 5).

The redshift distributions of the foreground and background sources are shown in Fig. 1. The bins of redshift adopted in the foreground sample are highlighted with varying shades of the same colour. The three different sets of bins, consisting of four bins, six bins and six bins in a wider range, are displayed from top to bottom in blue, indigo, and purple, respectively.

Fig. 1. Redshift distributions of the background SMGs (in red) and foreground lenses with the adopted redshift bins indicated by different shades of the same colour. From top to bottom, the “base case” is shown in blue, the “six bins” case in the 0.1 < z < 0.8 range in indigo and the “six bins-WR” case in the wider 0.01 < z < 0.9 range in purple, respectively.

2.2. Measurements

In previous works, the available area is divided in small regular regions, mini-tiles, due to its computational efficiency and accuracy in measuring the angular clustering of galaxies (e.g. González-Nuevo et al. 2017, 2021; Bonavera et al. 2019, 2020, 2021, 2022; Cueli et al. 2021, 2022). However, it requires an integral constraint (IC) correction to obtain an unbiased estimate of the true function, particularly crucial for constraining cosmological parameters on the largest angular scales. Moreover, using mini-tiles poses a practical challenge as the vertices of each of them must be manually defined, which is poorly scalable when expanding the study to other wide field surveys observed by Herschel that may include hundreds of mini-tiles.

To overcome these challenges, we have adopted a statistically rigorous approach that leverages the full field area and combines the number of different foreground-background pairs from each field into a single estimation. This new approach, described in Paper I (see also González-Nuevo et al. 2023, for a more detailed discussion), reduces statistical uncertainty and effectively utilises all available information in the data, as opposed to separately measuring cross-correlation functions for individual fields. The cross-correlation function is then measured using the modified version of the Landy & Szalay (1993) estimator by Herranz (2001),

$$\begin{array}{c}\hfill {\widehat{w}}_{\mathrm{fb}}(\theta )=\frac{D_{\mathrm{f}}{D}_{\mathrm{b}}(\theta )-D_{\mathrm{f}}{R}_{\mathrm{b}}(\theta )-D_{\mathrm{b}}{R}_{\mathrm{f}}(\theta )+R_{\mathrm{f}}{R}_{\mathrm{b}}(\theta )}{R_{\mathrm{f}}{R}_{\mathrm{b}}(\theta )}·\end{array}$$(1)

For the angular separation, θ, the values of D_fD_b, D_fR_b, D_bR_f, and R_fR_b represent the normalised foreground-background, foreground-random, background-random, and random-random pair counts coming from the whole available area. The background random samples are generated taking into account the surface density variation due to the Herschel scanning strategy.

Furthermore, to estimate the covariance matrix, we employed a Bootstrap method by dividing each field into at least five patches, ensuring a sufficient number of measurements. Instead of using the already mentioned mini-tiles, we have implemented a more adaptable procedure for defining subregions that retains similar characteristics to the previous approach but can be easily scaled for future studies. As explained in Paper I, these patches were created using a k-mean clustering algorithm, ensuring that the spatial dependence structure is retained during resampling. For defining the subregions, we took inspiration from TreeCorr, a widely used software package for gauging galaxy clustering (Jarvis 2015). TreeCorr employs a k-means clustering algorithm to divide the data into subregions referred to as “patches”, which are similar in size and shape to our mini-tiles. Specifically, we utilise the k-means algorithm available in the SciPy library. The objective of the algorithm is to minimise the sum of squared distances between data points and their assigned centroid. To determine the optimal number of clusters (i.e. patches), we enforced a minimum area requirement for each cluster. Additionally, we introduced an extra step by repeating the procedure ten times with different random initial centroids and selecting the run that demonstrated the most consistent number of data points across the various clusters.

This new methodology has demonstrated enhanced robustness and significantly reduced uncertainties in the obtained measurements. In essence, the overall area was divided into 22 subregions and N_b = 10 000 Bootstrap samples were created with an oversampling factor of 3. Therefore, the covariance matrix is given by

$$\begin{array}{c}\hfill \mathrm{Cov}({\theta }_i,{\theta }_j)=\frac{1}{N_{\mathrm{b}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{b}}}[{\widehat{w}}_k({\theta }_i)-\overline{\widehat{w}}({\theta }_i)][{\widehat{w}}_k({\theta }_j)-\overline{\widehat{w}}({\theta }_j)],}\end{array}$$(2)

where ${\widehat{\mathit{w}}}_k$ denotes the measured cross-correlation function from the kth Bootstrap sample and $\overline{\widehat{\mathit{w}}}$ is the corresponding average value over all Bootstrap samples. It should be noted that in this work, no attempt of correction of the super-sample covariance (SSC, e.g. Lacasa & Kunz 2017) is made since our main source of uncertainty at the moment is cosmic variance. It is due to the non-negligible variations that we have among the spatially separated fields (G09, G12, G15, and SGP), which is studied and commented on in Paper I.

In this work, we conduct a tomographic analysis by calculating the cross-correlation functions for each redshift bin defined in the previous section. The cross-correlation measurements, along with the outcomes of each case, are presented in Sect. 3. Specifically, in Fig. 2, we compare the cross-correlation measurements obtained by Bonavera et al. (2021; light grey circles) with those derived using our new methodology (black circles), thereby confirming the improved results also in the context of tomography.

Fig. 2. Cross-correlation data (black points) and the results obtained with the tomographic case: the posterior sampling of the cross-correlation function is shown in red and the best fit with the black dashed line. Bins 1 to 4 are shown from left to right and top to bottom. For comparison, the cross-correlation data by Bonavera et al. (2021) are shown in grey.

2.3. Theoretical framework

The low redshift galaxy-mass correlation and the cross-correlation between foreground and background sources are connected through weak lensing. The mass density field traced by the foreground galaxy sample causes weak lensing and this affects the number counts of the background galaxy sample through magnification bias, as described in detail in Paper I. To compute the correlation between the foreground and background sources, we adopt the halo model formalism described in Cooray & Sheth (2002) and the Limber and flat-sky approximations. The correlation can be evaluated using the following equation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill w_{\mathrm{fb}}(\theta )=& 2(\beta -1){\displaystyle {\int }_0^{\infty }\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{g}-\mathrm{m}}(l/\chi (z),z)J_0(l\theta ),}\hfill \end{array}\end{array}$$(3)

where W^lens(z) is defined as:

$$\begin{array}{c}\hfill W^{\mathrm{lens}}(z)=\frac{3}{2}\frac{H_0^2}{c^2}{E}^2(z){\displaystyle {\int }_z^{\infty }\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}$$(4)

Here, $E(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$, dN_b/dz, and dN_f/dz are the unit-normalised background and foreground redshift distributions, χ(z) is the comoving distance to redshift, z, and P_g − m is the galaxy-matter cross-power spectrum. The logarithmic slope of the background sources number counts is denoted by β, and the details of our handling of such parameter in the analysis are described in Paper I.

Due to the considerable number of integrals that need to be computed for each evaluation of the model, a mean-redshift approximation has been utilised. Such an approximation speeds up the computational time dramatically, by a factor of 10, enabling us to perform a variety of different tests in a reasonable time that would have otherwise been impossible due to our computer resources. Moreover, given the small redshift ranges used in tomography, we anticipate that this approximation will remain useful in the future, where improved parameter constraints are expected due to advancements in data, such as those from new surveys.

This approximation refers to the evaluation of the outermost integrals in Eq. (3) not being carried out directly for the entire redshift distribution, but instead done through computation at the sample’s mean redshift:

$$\begin{array}{c}\hfill w_{\mathrm{fb}}(\theta )\approx 2(\beta -1)\frac{W^{\mathrm{lens}}(\overline{z})}{{\chi }^2(\overline{z})}{\displaystyle {\int }_0^{\infty }\frac{l\mathrm{d}l}{2\pi }{P}_{\mathrm{g}-\mathrm{m}}(l/\chi (\overline{z}),\overline{z})J_0(l\theta ).}\end{array}$$(5)

As detailed in Paper I, the cross-power spectrum of galaxies and matter is expressed as the sum of two terms: the 1-halo term, which describes galaxy-matter correlations within the same halo, and the two-halo term, which refers to the cross-correlation between different halos.

The cross-correlation function depends on the value of beta and the cosmological parameters of the assumed underlying model, but it also needs information about the way galaxies populate dark matter halos. This is described by the halo occupation distribution (HOD), for which we assume the simple three-parameter model by Zehavi et al. (2005). According to it, when the mass of a halo exceeds a certain threshold, M_min, a galaxy is located at its centre. Any other galaxies are considered satellites and their distribution is directly proportional to the halo mass profile, as detailed in works such as Zheng et al. (2005). If the mass of a halo exceeds a different threshold, M₁, it will host satellites, and the number of satellites present is described by a power-law function with a coefficient of α.

Consequently, it is possible to represent the probability of a central galaxy being present as a step function:

$$\begin{array}{c}\hfill N_{\mathrm{c}}(M)=\{\begin{array}{c}0\mathrm{if}M<M_{\min },\hfill \\ 1\mathrm{otherwise}.\hfill \end{array}\end{array}$$(6)

The satellite galaxies occupation can be described as:

$$\begin{array}{c}\hfill N_{\mathrm{s}}(M)=N_{\mathrm{c}}(M)·{\left(\frac{M}{M_1}\right)}^{\alpha },\end{array}$$(7)

with M_min, M₁, and α the free-parameters of the model.

It is important to emphasise that the analysis only considers cross-correlation data in the weak lensing regime, namely, θ ≥ 0.5 arcmin since we are working within such approximation in our theoretical model (for a more comprehensive discussion, refer to González-Nuevo et al. 2017; Bonavera et al. 2019). Moreover, as discussed in Paper I, the validity of the standard halo model at small scales should initially be taken with caution. Physical effects such as baryonic feedback could affect the 1-halo cross-correlation function via the modification of the galaxy-matter power spectrum (van Daalen et al. 2014, 2020; Renneby et al. 2020; Amon et al. 2023). An approach that has been used in the literature to address this issue implies adopting a different normalisation for dark matter halo profiles and for the satellite galaxy distribution (Cacciato et al. 2013; Viola et al. 2015; van Uitert et al. 2016; Dvornik et al. 2018, 2023). We checked that introducing this modification affects only the sub-arcmin data points. However, these variations can be regarded as negligible when compared to the current uncertainties in the observational data. Therefore, in this work we do not discuss the effect of the baryonic feedback, which should be taken into account in future works, when more constraining results will be achieved.

2.4. Parameter estimation

In order to estimate the parameters, we use the open source software package called emcee, which is licensed under MIT, and utilises a Markov chain Monte Carlo (MCMC) algorithm to estimate the parameters. It is based on the Affine Invariant MCMC Ensemble sampler by Goodman & Weare (2010) and is implemented purely in Python. The MCMC runs in this work are performed with a number of walkers set to three times the number of parameters to be estimated, N, and a fixed number of iterations set to 5000. This generates 15 000 N posterior samples per run, which are then reduced by discarding the initial iterations and introducing a burn-in when flattening the chain. We examined each parameter’s walkers separately to identify the burn-in value. We selected the iterations where all walkers had moved to explore a specific area rather than wandering randomly. The number of iterations we kept depended on the specific case. Finally, we confirmed the chains’ convergence using the common tests provided by GETDIST (Lewis 2019).

In our analysis, we aim to estimate both astrophysical and cosmological parameters. The astrophysical parameters of interest are M_min, M₁, and α, with different values for each redshift bin. With respect to Bonavera et al. (2021), the parameter β is also added to the analysis, since, as studied in Paper I, fixing it to the previously used value of 3 could affect the results. On the other hand, the cosmological parameters include Ω_M, σ₈, and h, which is defined as H₀ = 100 × h km s⁻¹ Mpc⁻¹. We assume a flat universe where Ω_DE = 1 − Ω_M, and we hold the values of Ω_b and n_s fixed at their best-fit values from Planck Collaboration VI (2021), which are Ω_b = 0.0486 and n_s = 0.9667, respectively. In the weak lensing approximation, only the cross-correlation function data with angular scales greater than or equal to 0.5 arcmin are considered, as discussed in detail in Bonavera et al. (2019).

When using the Bayesian method for parameter estimation, it is necessary to define both a prior and a likelihood distribution. For the likelihood distribution, a conventional Gaussian function is employed. In terms of selecting prior distributions for the astrophysical parameters, we used uniform priors. In particular, in the scenarios using the four bins, that is “base case”, “w/o G15”, “bins 1+4”, and “bins 2+3”, we adopted the log M_min and log M₁ priors listed in Table 1. In the “six-bin” case we use the priors specified in Table 2. For the “six bins-WR” case, we give those in Table 3. For all bins and all cases, α ∼ 𝒰[0.5, 1.5] and β priors are Gaussian with β ∼ 𝒢[2.9, 0.04], following the discussion in Paper I.

The prior distributions for the cosmological parameters that we aim to estimate are the same as those employed in Bonavera et al. (2020, 2021). More specifically, Ω_M ∼ 𝒰[0.1 − 0.8], σ₈ ∼ 𝒰[0.6 − 1.2], and h ∼ 𝒰[0.5 − 1.0].

3. Results

In this section we describe and discuss the results obtained by jointly analysing different combination of bins of redshift for the foreground sample, the lenses. In particular, Sect. 3.1 describes the “base case”, while Sect. 3.2 presents a study of the potential effect in the results due to the anomalous measurements in the G15 field, pointed out in Paper I. Finally, in Sect. 3.3 the effect on the constraints of using different number and ranges of bins for the tomographic analysis is discussed. On one hand, we address the possibility of reducing the joint analysis to only a pair of bins (“bins 1+4” and “bins 2+3”), which would diminish the computational time of the parameters estimation. On the other hand, we increase the number of bins to six to explore both the effect of a thinner slicing (“six bin” case) of the foreground sample in the same redshift range of the “base case” and the effect of introducing two more bins, one at each end of the four bin cases (“six bins-WR” case).

3.1. Base case: Four redshift bins

The first step of our study consists in the comparison between the results obtained with the “full model” and the mean-redshift approximation, as described in Paper I. This approximation reduces the computational time by about a factor of 10 and it is of crucial importance when performing tomographic analyses, given the large number of parameters to be estimated in such a configuration. The same improved data and methodology described in Paper I are used in both cases. The comparative plots are shown in Appendix A. The HOD parameters prior distributions used for both cases are summarised in Table 1.

Figures A.1 and 2 present the posterior sampling of the cross-correlation function (in red), the best fit (dashed black line), and the data (black dots) in the four redshift bins for the “full model” and the mean-redshift approximation, respectively. The panels represent the results for the different bins of redshift: bins 1 to 4 from left to right, top to bottom, respectively. The redshift range is indicated in every panel. The data have comparable errorbars in the four bins – except in bin 1, which present larger errorbars at intermediate and large angular scales. Overall, the best fit and sampled area area in good agreement with the data in both cases.

Moreover, Fig. A.2 illustrates the marginalised posterior distribution and probability contours¹ for all the estimated parameters. Given the large number of involved parameters, in order to clearly show the results, we divide the corner plot according to the number of bins. In this way we depict the posterior distribution of the HOD parameters of each bin together with the cosmological parameters (that are the same for every bin). With the same aim of making the presentation of the results more intelligible, we decide not to include in the plots the β parameter, since its posterior distribution barely deviates from its (Gaussian) prior. As expected, given the small bin ranges used for tomography, the mean-redshift approximation (“base case”) results are in good agreement with the “full model” case.

The top panel of Fig. 3 displays the marginalised posterior distributions of the HOD parameters for the “base case”. The solid black, dashed red, dotted blue, and dash-dotted green lines correspond to bin 1, bin 2, bin 3, and bin 4, respectively. A notable trend can be observed, with the value of log M_min increasing as a function of redshift from bin 1 to bin 4. The distribution peaks are located at 11.35, 11.45, 11.87, and 12.71, respectively². This behaviour is consistent with previous tomographic analyses (González-Nuevo et al. 2017; Bonavera et al. 2021; Cueli et al. 2022).

Fig. 3. HOD parameters for the tomographic case (top). HOD parameters for the tomographic case excluding the data from the G15 field (bottom). The solid black, dashed red, dotted blue and dash-dotted green lines are the bin 1, bin 2, bin 3, and bin 4, respectively. The mass values are expressed in M⊙ h−1 units.

Regarding the log M₁ parameter, the marginalised posteriors appear relatively homogeneous across the four bins, with peaks ranging from 13.22 to 13.58; with the exception of bin 1, which exhibits slightly higher values peaking at 14.08 along with a broader distribution. On the other hand, the α parameter remains poorly constrained in all bins.

Comparing our results to those of Bonavera et al. (2021), we generally find narrower distributions and consistent findings. The main discrepancy lies in the log M_min distributions for bin 4, where the values obtained in the “base case” are lower. Additionally, some subtle distinctions can be observed in the opposite trends of the α distributions for bin 2, although it should be noted that α is not constrained in either case.

Figure 4 presents the estimated cosmological parameters for the “base case” represented by the blue curves. The marginalised distributions exhibit peaks at 0.16 for Ω_M and 0.70 for σ₈. The results from the non-tomographic analysis in Paper I and the tomographic analysis by Bonavera et al. (2021) are also displayed in red and green, respectively. The impact of the new methodology, particularly the reduction in measurement uncertainties, is evident for Ω_M, as indicated by the significantly narrower posterior distributions compared to the results obtained by Bonavera et al. (2021). However, it is noteworthy that there is a preference for lower values of Ω_M, with the posterior distributions deviating considerably from the standard model. This issue will be further explored in the subsequent subsection.

Fig. 4. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the cosmological parameters obtained in the “base case” in blue, with single bin from Paper I in red and by Bonavera et al. (2021) in green.

Regarding σ₈, there is also a preference for lower values, with the newer analysis exhibiting a slightly narrower posterior distribution. The higher value observed in the non-tomographic case is consistent with the comparison between the non-tomographic results from González-Nuevo et al. (2021) and the tomographic analysis by Bonavera et al. (2021). Notably, the parameter h remains unconstrained in all three cases.

3.2. Without G15: Four redshift bins without the G15 field

The inclusion or exclusion of the G15 field has a significant impact on the measured cross-correlation function at larger angular scales, as previously discussed in Paper I. Since the cosmological constraints heavily rely on these angular scales, the effect of the G15 field on the results becomes crucial, as demonstrated in the same work for a broad single redshift bin. In light of this, we explore the potential effects of excluding the G15 field in the tomographic analysis in this study.

The posterior sampling of the cross-correlation function (shown in red in Fig. 5) confirms the effect of removing the G15 field on the largest angular scales, resulting in a steeper decrement of the signal. Due to the reduced statistics compared to the case where all four fields are included, the uncertainties of the data are larger. Comparing both sets of measurements, it is evident that the removal of the G15 field leads to a worsening of the best fit and posterior sampling for bin 1 and, especially, for the highest redshift bin. This can be attributed to the statistical weight of the central bins in the overall best fit due to the relatively better uncertainties in the measurements.

Fig. 5. Cross-correlation data (black points) and the results obtained with the tomographic case excluding G15: the posterior sampling of the cross-correlation function is shown in red and the best fit with the black dashed line. Bins 1 to 4 are shown from left to right, top to bottom, and the redshift ranges are indicated in each panel.

Figure B.1 displays the posterior distribution of the estimated parameters, organised per bin, similar to Fig. A.2. The marginalised posterior distributions of the estimated HOD parameters, presented in the bottom panel of Fig. 3, reveal that (similarly to the “base case”) the log M_min parameters exhibit increasing peak values with the bin range. Specifically, the distributions peak at 10.70, 11.94, 12.05, and 12.74 for bins 1 to 4, respectively. Comparing with the “base case” in the upper panel, the main difference lies in the log M_min distribution of bin 2, which shifts to higher values and of bin 1, which shifts to lower values, when the G15 field is excluded. This shift may be related to the substantial difference in the signal levels observed in bin 1 and bin 2 between these two cases. The parameter log M₁ displays a worse but similar behaviour to that of the “base case”, while the α distributions are mostly unconstrained and with slightly opposite trend than in the “base case”.

Figure 6 presents the marginalised posterior distributions and probability contours of the cosmological parameters estimated for the “w/o G15” case (in green). The marginalised posterior distribution for Ω_M peaks at a higher value of 0.32 compared to the “base case” (in blue), while it peaks at a lower value of 0.60 for σ₈. The parameter h remains unconstrained. In the same plot, we also compare the results with the corresponding non-tomographic analysis from Paper I (in red). The main difference lies in the higher value of σ₈, which peaks at approximately 0.73 in the non-tomographic case, but the qualitative behaviour is the same: the removal of the G15 region produces lower values of σ₈ and higher values of Ω_M both in the non-tomographic and tomographic cases.

Fig. 6. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) of the cosmological parameters obtained in the “w/o G15” case (in green) and in the corresponding case with single bin from Paper I (in red). The “base case” results are also shown in blue.

In summary, we confirm in this study that the exclusion of the G15 field leads to a decrease in the cross-correlation signal at the largest angular separation distances, as observed in Paper I. While this effect is noticeable in all redshift bins, it is particularly evident across all angular scales in bins 1 and 2. Consequently, the main difference with respect to the “base case” is observed in the log M_min distribution for those bins, which shifts towards significantly lower and higher values, respectively. These findings suggest that the anomaly associated with the G15 field may be related to a sample variance problem within the large-scale structure in the line of sight, but also within a specific range of redshifts.

3.3. Tomographic analysis with different bins of redshift

In this section, we assess the impact of considering different numbers of redshift bins and exploring different redshift ranges. One interesting aspect to investigate from a computational perspective is the possibility of using only two bins of redshift, which would effectively halve the number of HOD parameters to be estimated compared to the “base case” (see Sect. 3.3.1). Additionally, we aim to discern the major contributor to cosmological parameter estimation: whether it is more influenced by the extreme redshift bins (“bins 1+4” case), where the redshift evolution of the measurements plays a prominent role, or the central bins (“bins 2+3” case), which offer better statistical robustness.

In Sect. 3.3.2, despite not reducing computational time, we further explore the possibilities with a finer slicing of the default redshift range. Specifically, we investigate the potential benefits of jointly analysing six bins of redshift (“six bin” case). Lastly, we extend the total redshift range by incorporating two additional bins at the endpoints of the original “base case” in the foreground sample. Through these analyses, we aim to gain insights into the implications of different redshift binning strategies and explore the trade-offs between computational efficiency and precision in cosmological parameter estimation.

3.3.1. Two-bin analysis

The panels in Fig. 7 depict the data (black dots), the posterior sampling of the cross-correlation function (in red), and the best fit (dashed black line) obtained in the joint analysis of bin 1 and bin 4 (“bins 1+4”, left panel) as well as in the joint analysis of bin 2 and bin 3 (“bins 2+3”, right panel). In both cases, the upper panel refers to the lowest bin in the pair and the bottom to the highest. The redshift range is indicated in every panel. In most cases the best fit and sampled area are in agreement with the data. The marginalised and probability distribution are shown in Figs. B.2 (for the joint bin 1 and bin 4 case) and B.3 (for the joint bin 2 and bin 3 case). As with the previous plots, the results are divided according to the number of jointly analysed bins – two in this case – for an improved visualisation.

Fig. 7. Cross-correlation data (black points), best fit (black dashed line) and the posterior sampling (in red) in the “bins 1+4” scenario (left panel) and in the “bins 2+3” case (right panel).

The HOD marginalised posterior distributions for the two-bins cases are shown in Fig. 8: the red dotted line indicates the results of “bins 1+4” and the blue dashed line those of “bins 2+3”. These cases are compared with the “base case” (black solid line) and the common bins in the “six bins-WR” case. This figure shows the three analysed HOD parameters log M_min, log M₁, and α in the columns from left to right, respectively. The four rows represents the four bins (bins 1 to 4) from top to bottom.

Fig. 8. Comparison of the marginalised posterior distributions of the HOD parameters (log Mmin, log M1 and α from left to right, respectively) derived in the different analysed cases, for bin 1, bin 2, bin 3, and bin 4, respectively (from top to bottom). The “base case” is shown with the black solid line, “bins 1+4” with the red dotted line, “bins 2+3” with the blue dashed line and the “six bins-WR” with the green dash-dotted line. The mass values are expressed in M⊙ h−1 units.

For the “bins 1+4” case, the M_min parameter is slightly shifted towards lower values with that of the “base case”, peaking at 11.16 for bin 1 and 12.46 for bin 4. The log M₁ parameter aligns with the “base case” peaking at 13.38 in bin 4 but it is slightly higher, 14.91, in bin 1. The α parameter remains unconstrained in bin 1 and provides a lower limit in bin 4. In the “bins 2+3” case, log M_min agrees with the “base case”, with peaks at 11.44 in bin 2 and 11.86 in bin 3. The log M₁ results are also consistent, peaking at 13.10 in bin 2 and in bin 3. The parameter α exhibits a lower limit in bin 2 and shows a tentative peak at 1.23 in bin 3.

The posterior distributions of cosmological parameters are shown in Fig. 9, with the “bins 1+4” case displayed in red and the “bins 2+3” case in green. For comparison, the “base case” is plotted in blue, with uncertainties only slightly improved compared to the other two cases. The marginalised posterior distribution of Ω_M shows good agreement among the three cases, peaking at 0.16. The σ₈ posteriors reveal an intriguing displacement: the distribution for “bins 2+3” is shifted to lower values (peaking at 0.61) compared to “bins 1+4” (peaking at 0.82), while the “base case” distribution falls in between (peaking at 0.70). The parameter h remains unconstrained as expected. Tomographic analysis achieves substantial constraining power even with only two redshift bins. However, the discrepancy in σ₈ across the different cases suggests a sample variance issue that can be mitigated by increasing the number of redshift bins. Thus, the sample variance problem identified in Paper I, which affected the spatial distribution, is also manifested as sample variance in redshift, as hinted at in the analysis of the G15 effect in Sect. 3.2.

Fig. 9. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) of the cosmological parameters obtained in the “base case” (in blue) and in the two bins cases, “bins 1+4” (in red) and “bins 2+3” (in green).

3.3.2. Six-bin analysis

Considering the conclusions of the previous section, we investigate the effect of increasing the number of redshift bins in two specific cases: maintaining the same redshift range (“six bins”) or extending it by adding additional redshift bins (“six bins-WR”). The definition of the redshift bins for both cases was already discussed in Sect. 2.1 and the prior distributions used are summarised in Tables 2 and 3, respectively.

The panels in Fig. 10 display the posterior sampling of the cross-correlation function (in red), the best fit (dashed black line), and the data (black dots) for the joint analysis of six bins in the original redshift range (“six bins” case, left plot) and the wider redshift range from 0.01 to 0.9 (“six bins-WR” case, right plot). Each panel represents the results for a specific redshift bin, and the range of each bin is indicated. In the “six bins” case, the data are visibly affected by poorer statistics compared to the “base case” due to the same total number of sources being spread over more bins. In the “six bins-WR” case, the central bins show similar behaviour to the “base case”, while the newly added first and last bins have large error bars due to lower lensing probability (see Lapi et al. 2012) in the first bin and significantly poorer statistics in the last bin. In most cases, the best fit and sampled area align well with the data, except for the last bin in the “six bins-WR” case, which exhibits poor agreement.

Fig. 10. Cross-correlation data (black points), best fit (black dashed line) and sampling posterior (in red) obtained in the “six bin” case (left) and the “six bins-WR” case (right).

The marginalised and probability distribution of the six-bins analysis are illustrated in Figs. B.4 (“six bins” case) and B.5 (“six bins-WR” case). As in the previous cases, the corner plot is divided according to the number of jointly analysed bins (six in this case) for a better visualisation.

Figure 11 displays the marginalised posterior distributions of the HOD parameters for the “six bins” cases: the “six bins” case in the upper panel and the “six bins-WR” case in the lower panel. The parameters log M_min, log M₁, and α are shown from left to right, and different redshift bins are indicated with various line styles. In the “six bins-WR” case, the most similar bins to the original ones are represented with the same line style. The additional bins added at the extremes of the original range are depicted with the yellow dash-dotted line (bin 0) and the purple dash-dot-dotted line (bin 5) in the bottom panel.

Fig. 11. HOD parameters marginalised posterior distributions: log Mmin, log M1 and α in the first, second, and third column, respectively. The top panel shows the results for the “six bin” case and the bottom panel those for the “six bins-WR” case. The mass values are expressed in M⊙ h−1 units.

The M_min parameter generally increases with redshift, except for bin 2a in the “six bins” case, which shows a distribution similar to bin 1, as well as bin 5 in the “six bins-WR” case, positioned between bins 3 and 4. This might simply be due to the much larger errobars and less significant signal in such bins, especially bin 5 in the “six bins-WR” case, which is the one with the worst fit. All the distributions are well constrained showing peaks (from lower to higher bin) at 11.23, 11.07, 11.66, 11.77, 12.31 and 12.66 for the “six bins” case and at 10.32, 11.53, 11.75, 11.95, 13.36 and 12.09 for the “six bins-WR” case. Notably, in the “six bins-WR” case, the values tend to converge around 12.0 for all redshift bins except bin 0 and bin 4. Figure 8 compares, for the redshift bins in common, the HOD marginalised posterior distributions of the “six bins-WR” case (green dash-dotted line) with the other three cases studied in this work. The figure shows again how the log M_min values tend to converge toward a value near 12.0.

The marginalised posterior distributions of log M₁ for the “six bins” show very similar peaks that range from 13.85 to 14.55, except for the last bin. In the “six bins-WR” case, the distributions do not converge and they do not exactly follow the same behaviour with redshift as in the log M_min ones. Their peaks range from 12.37 for bin 0 to 14.96 for bin 5. The parameter α does not show a clear trend, as the posterior distributions peak at different values within the available range for both cases.

The marginalised posterior distributions of the cosmological parameters are presented in Fig. 12 in red for the “six bins” case and in green for the “six bins-WR” case. The “base case” results are plotted in blue for comparison. The marginalised posterior distribution of Ω_M clearly shifts towards higher values, with the “six bins” case peaking at 0.20 and the “six bins-WR” case peaking at 0.32, compared to the “base case” peak at 0.16. The posterior distributions in both cases are wider than in the “base case”.

Fig. 12. Comparison of the marginalised posterior distribution and probability contours for the cosmological parameters obtained in the “base case” (in blue) and in the six bins cases: “six bins” (in red) and “six bins-WR” (in green).

Increasing the number of redshift bins has a similar effect to removing the G15 field, especially when additional bins are added instead of dividing the existing ones. This suggests that the evolution with redshift, i.e. the cosmological model, starts to dominate over the sample variance effect. However, the overall fit is not perfect for each individual bin, resulting in poorer fits for the most biased samples. Extending the redshift range seems to work better than simply increasing the number of redshift bins, as it improves statistics per bin and produces tighter constraints compared to the “w/o G15” case.

The σ₈ posterior distribution for the “six bins-WR” case is compatible with the “base case” results but shows a wider distribution towards lower values, with a peak (mean) value at 0.66 (0.69). In contrast, the “six bins” case exhibits a clear peak at 0.99. This difference can be attributed to the splitting of the two central bins, which reduces their statistics compared to bin 1 and bin 4, resulting in a lower weight in the parameter estimation and producing estimates similar to the “bins 1+4” case. The secondary peak observed between 0.7 and 0.8 in the “six bins-WR” case is of particular interest. We have confirmed the stability of the full distribution, ruling out any convergence issues. Instead, we interpret the secondary peak as a local valid minimum within the constraints of our current data.

Finally, the parameter h remains unconstrained, but for the first time in this work, the posterior distributions show tight upper limits of 0.69 at 95% C.I., in the “six bins-WR”. These results align with the findings in Paper I, where a joint analysis of the cross-correlation function and galaxy clustering measurements was performed. Although these results should be interpreted with caution, they indicate that with improved data and better-constrained parameters, there is hope for obtaining constraints on h.

4. Conclusions

This study and its companion paper examine various aspects associated with the study of the submillimeter galaxy magnification bias, with a specific focus on the tomographic scenario. It is akin to the work carried out by Bonavera et al. (2021), but using the updated data from Paper I. The objective is to refine the methodology employed in constraining the free parameters of the HOD model, that is, M_min, M₁, and α, and the cosmological parameters Ω_M, σ₈, and h of a flat ΛCDM model on a tomographic scenario.

Due to the large number of parameters that must be jointly analysed in tomography, a significant amount of computational time is required to accomplish this kind of task. Therefore, one of the aims of this work is to explore ways to optimise CPU time and strategies for analysing different redshift bins without compromising the precision of the measurements. Various data sets have been explored to this end, which involve the combination of different redshift bins and the inclusion or exclusion of the G15 field suspected to have an anomalous behaviour. Each data set consists in the measurement of the cross-correlation function between background H-ATLAS submillimeter galaxies and foreground GAMA galaxies within a tomographic setup.

As in Bonavera et al. (2021) and Cueli et al. (2022), we divided the entire redshift range into four bins: bin 1 (0.1 < z < 0.2), bin 2 (0.2 < z < 0.3), bin 3 (0.3 < z < 0.5), and bin 4 (0.5 < z < 0.8). In the “base case” and “w/o G15” scenarios, we examine all four bins together, including SGP and all the GAMA fields and excluding the G15 field, respectively.

We first compare the “base case” which makes use of the mean-redshift approximation with the “full model” case. This approximation reduces the computational time by about a factor of 10 and it is of crucial importance when performing tomographic analyses, given the large number of parameters to be estimated in such a configuration. The comparative plots confirm that given the thinness of the redshift bins used for tomography, the “base case” results are in good agreement with the “full model” case.

The marginalised posterior distributions exhibit a systematic increase in log M_min from bin 1 to bin 4 as redshift increases, while the distributions of log M₁ show relatively consistent values across the four bins. In addition to the HOD parameters, we compare the inferred cosmological parameters with the results obtained from the non-tomographic runs in Paper I and the tomographic analysis conducted by Bonavera et al. (2021). The impact of our new methodology, particularly the reduction in measurement uncertainties, becomes evident when examining Ω_M, as indicated by the considerably narrower posterior distributions compared to the findings of Bonavera et al. (2021). However, it is worth noting that our results exhibit a preference for lower values of Ω_M, with the posterior distributions deviating noticeably from the predictions of the standard model.

In relation to the “w/o G15” scenario, our study confirms the findings of Paper I, which demonstrate that excluding the G15 field results in a reduction of the cross-correlation signal at larger angular separations. While this effect is observed across all redshift bins, it is particularly pronounced in bins 1 and 2, affecting all angular scales. As a result, the most significant deviation from the “base case” is observed in the log M_min distribution for these bins, which shifts towards lower and higher values, respectively. These observations suggest that the peculiar behavior associated with the G15 field might be attributed to sample variance within the large-scale structure along the line of sight, particularly within a specific range of redshifts.

For the sake of computational efficiency, we narrow our focus to the joint analysis of two bin pairs: “bins 1+4” (the first and fourth bin) and “bins 2+3” (the second and third bin). Remarkably, the HOD parameters obtained from these analyses are consistent with those of the “base case”. The marginalised posterior distribution of Ω_M exhibits strong agreement across all three cases, peaking at 0.16. However, the σ₈ posteriors reveal an intriguing disparity: the distribution for “bins 2+3” is shifted towards lower values (peaking at 0.61), whereas “bins 1+4” yields higher values (peaking at 0.82), and the “base case” falls in between (peaking at 0.70). This discrepancy suggests the presence of sample variance in redshift, which can be mitigated by increasing the number of redshift bins. Thus, the sample variance issue identified in Paper I, which affected the spatial distribution, also manifests as sample variance in the redshift domain, as implied by the analysis of the G15 effect.

We investigate the impact of increasing the number of redshift bins in two cases: “six bins” (maintaining the same redshift range) and “six bins-WR” (extending the range with additional bins). In the “six bins” case, the data are affected by poorer statistics due to the same number of sources being spread over more bins. In the “six bins-WR” case, the central bins resemble the “base case”, while the newly added first and last bins have larger error bars and significantly poorer statistics.

The log M_min distributions are well-constrained, with peaks increasing with redshift, as expected. In the “six bins-WR” case, the values show a convergence around the value 12.0 for most bins. The marginalised posterior distributions of log M₁ tend to peak around the value 14 for all the bins in both cases, although there is more dispersion in the “six bins-WR” case. The Ω_M distribution shifts towards higher values, peaking at 0.20 in the “six bins” case and 0.32 in the “six bins-WR” case, compared to the “base case” peak at 0.16. Increasing the number of redshift bins has a similar effect as removing the G15 field, suggesting the dominance of the cosmological model over sample variance. However, individual bin fits are not perfect, resulting in poorer fits for the most biased samples. Extending the redshift range improves statistics and produces tighter constraints compared to the “w/o G15” case.

The σ₈ posterior distribution for the “six bins-WR” case is compatible with the “base case” results but shows a wider distribution towards lower values, with a peak (mean) value of 0.66 (0.69). In contrast, the “six bins” case exhibits a clear peak at 0.99. This difference is due to reduced statistics in the two central bins compared to bins 1 and 4, producing estimates more similar to the “bins 1+4” case. The parameter h remains unconstrained, but the posterior distributions show tight upper limit for the “six bins-WR” case, consistent with previous findings in Paper I.

As a concluding remark, this study helps to refine the methodology for tomographic analyses and highlights the importance of considering sample variance and redshift binning in obtaining precise measurements of HOD and cosmological parameters. The spatial sample variance as well as (most likely) the redshift one can be minimised by increasing the number of independent fields, for instance, by incorporating additional Herschel wide area surveys such us the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) or the Great Observatories Origins Deep Survey-Herschel (GOODS-H; Elbaz et al. 2011). However, our results also indicate that enhancing both the number of redshift bins and the redshift range yields more accurate and robust outcomes, mitigating the potential impact of both types of sample variance. To implement this second approach effectively, it is important to adopt updated foreground catalogues with improved statistics for our purposes, such as the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005) or the future Euclid mission (Euclid Collaboration 2022).

Acknowledgments

L.B., J.G.N., J.M.C., and D.C. acknowledge the PID2021-125630NB-I00 project funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. L.B. also acknowledges the CNS2022-135748 project funded by MCIN/AEI/10.13039/501100011033 and by the EU “NextGenerationEU/PRTR”. J.M.C. also acknowledges financial support from the SV-PA-21-AYUD/2021/51301 project. We deeply acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. In particular the projects “SIS22_lapi”, “SIS23_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 website 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 website is: http://www.gama-survey.org. 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: Comparison "full model" and "base case"

This appendix presents the results obtained in the "base case" and the original model, "full model" which is used in previous works (e.g. González-Nuevo et al. 2017; Bonavera et al. 2020, 2021, 2022; Cueli et al. 2021, 2022). The jointly analysed bins of redshift and the priors are the same in both runs. Figure A.1 shows the cross-correlation data (black points), the best fit (black dashed line), and the sampled area (in red) obtained with the MCMC. The agreement with the data is good as well as with the findings in Fig. 2.

Fig. A.1. Cross-correlation data (black points) and the results obtained with the full model with no approximations: the black dashed line shows the best fit and the red lines are the sampling posteriors.

Moreover, we show in red the marginalised posterior distributions and probability contours (set to 0.393 and 0.865) in Fig. A.2. Even if the analysis is jointly performed among the four bins, due to the large amount of parameters to be shown in the figure, we split such results into four panels. Each panel shows the HOD estimated parameters for one bin together with the cosmological parameters, that are the same for the four panels. In particular, bins 1 to 4 are shown from left to right, top to bottom, respectively. In this figure we also show in blue the results of the "base case" where we adopt the mean-redshift approximation (detailed in Paper I), which consists of computing the outer integrals in Eq. 3 at the mean redshift of the sample (see Eq. 5). In the tomographic analysis performed in this work, we can safely adopt such approximation, given the small range of redshift in each in bin.

Fig. A.2. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "base case" (in blue) and in the "full model" case (in red).

The results we find based on these two cases are in very good agreement and the subsequent analyses in this work can be safely performed with it. It provides great advantages from the computational point of view, especially for tomographic analyses, which involve a large number of parameters.

Appendix B: Posterior distributions of the individual cases

This section presents the marginalised posterior distributions and probability contours (set to 0.393 and 0.865) of all the cases addressed in this work and compared to the "base case" presented in Sect. A. In particular, Fig. B.1 show the case corresponding to Sect. 3.2 in which the joint analysis is performed on the four bins bin 1, bin 2, bin 3, and bin 4 without taking into account the G15 field.

Fig. B.1. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "w/o G15" case.

The "two bins" case described in Sect. 3.3.1, where only two bins are jointly analysed, are shown in Figs. B.2 and B.3. In particular, Fig. B.2 corresponds to the "bins 1+4" case, where just bin 1 and bin 4 are jointly analysed and Fig. B.3 shows the results obtained in the "bins 2+3" case, that is, considering only bin 2 and bin 3.

Fig. B.2. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "bins 1+4" case.

Fig. B.3. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "bins 2+3" case.

The "six bins" cases described in Sect. 3.3.2 are depicted in Figs. B.4 and B.5, while Fig. B.4 describes the findings of the "six bins" case, which are attained using bin 1, bin 2a and bin 2b (by splitting bin 2), bin 3a and bin 3b (by splitting bin 3), and bin 4. Figure B.5 shows the "six bins-WR" case, obtained adding two extra bins to the "base case", that is: bin 0 (to the lower end of the original case, 0.01 > z > 0.1), bin 1, bin 2, bin 3, bin 4 (but only up to z = 0.7 to increase the statistics of bin 5), and bin 5 (0.7 > z > 0.9).

Fig. B.4. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "six bins" (i.e. 0.1 < z < 0.8) case.

Fig. B.5. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "six bins-WR" (i.e. 0.01 < z < 0.9) case.

All Tables

All Figures

	Fig. 1. Redshift distributions of the background SMGs (in red) and foreground lenses with the adopted redshift bins indicated by different shades of the same colour. From top to bottom, the “base case” is shown in blue, the “six bins” case in the 0.1 < z < 0.8 range in indigo and the “six bins-WR” case in the wider 0.01 < z < 0.9 range in purple, respectively.
In the text

	Fig. 2. Cross-correlation data (black points) and the results obtained with the tomographic case: the posterior sampling of the cross-correlation function is shown in red and the best fit with the black dashed line. Bins 1 to 4 are shown from left to right and top to bottom. For comparison, the cross-correlation data by Bonavera et al. (2021) are shown in grey.
In the text

	Fig. 3. HOD parameters for the tomographic case (top). HOD parameters for the tomographic case excluding the data from the G15 field (bottom). The solid black, dashed red, dotted blue and dash-dotted green lines are the bin 1, bin 2, bin 3, and bin 4, respectively. The mass values are expressed in M⊙ h−1 units.
In the text

	Fig. 4. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the cosmological parameters obtained in the “base case” in blue, with single bin from Paper I in red and by Bonavera et al. (2021) in green.
In the text

	Fig. 5. Cross-correlation data (black points) and the results obtained with the tomographic case excluding G15: the posterior sampling of the cross-correlation function is shown in red and the best fit with the black dashed line. Bins 1 to 4 are shown from left to right, top to bottom, and the redshift ranges are indicated in each panel.
In the text

	Fig. 6. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) of the cosmological parameters obtained in the “w/o G15” case (in green) and in the corresponding case with single bin from Paper I (in red). The “base case” results are also shown in blue.
In the text

	Fig. 7. Cross-correlation data (black points), best fit (black dashed line) and the posterior sampling (in red) in the “bins 1+4” scenario (left panel) and in the “bins 2+3” case (right panel).
In the text

Fig. 8. Comparison of the marginalised posterior distributions of the HOD parameters (log Mmin, log M1 and α from left to right, respectively) derived in the different analysed cases, for bin 1, bin 2, bin 3, and bin 4, respectively (from top to bottom). The “base case” is shown with the black solid line, “bins 1+4” with the red dotted line, “bins 2+3” with the blue dashed line and the “six bins-WR” with the green dash-dotted line. The mass values are expressed in M⊙ h−1 units.

In the text

	Fig. 9. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) of the cosmological parameters obtained in the “base case” (in blue) and in the two bins cases, “bins 1+4” (in red) and “bins 2+3” (in green).
In the text

	Fig. 10. Cross-correlation data (black points), best fit (black dashed line) and sampling posterior (in red) obtained in the “six bin” case (left) and the “six bins-WR” case (right).
In the text

	Fig. 11. HOD parameters marginalised posterior distributions: log Mmin, log M1 and α in the first, second, and third column, respectively. The top panel shows the results for the “six bin” case and the bottom panel those for the “six bins-WR” case. The mass values are expressed in M⊙ h−1 units.
In the text

	Fig. 12. Comparison of the marginalised posterior distribution and probability contours for the cosmological parameters obtained in the “base case” (in blue) and in the six bins cases: “six bins” (in red) and “six bins-WR” (in green).
In the text

	Fig. A.1. Cross-correlation data (black points) and the results obtained with the full model with no approximations: the black dashed line shows the best fit and the red lines are the sampling posteriors.
In the text

	Fig. A.2. Comparison of the marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "base case" (in blue) and in the "full model" case (in red).
In the text

	Fig. B.1. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "w/o G15" case.
In the text

	Fig. B.2. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "bins 1+4" case.
In the text

	Fig. B.3. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "bins 2+3" case.
In the text

	Fig. B.4. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "six bins" (i.e. 0.1 < z < 0.8) case.
In the text

	Fig. B.5. Marginalised posterior distribution and probability contours (set to 0.393 and 0.865) for the HOD and cosmological parameters obtained in the "six bins-WR" (i.e. 0.01 < z < 0.9) case.
In the text