© ESO 2022

1. Introduction

The concordance cosmological framework explains the process of galaxy formation via the hierarchical growth of structure from primordial Gaussian density perturbations. Indeed, virialized dark matter halos break away from the expansion of the background and provide the potential wells for gas to cool down and for galaxies to form. An accurate understanding of the abundance and evolution of dark matter halos is therefore of paramount importance to accurately explain the early Universe and to compute the nonlinear effects of structure formation.

The seminal paper of Press & Schechter (1974) laid the conceptual foundations of dark matter halo formation by relating their collapse to the moment that the spherically smoothed linear density field goes over a certain threshold value that depends on cosmology, thus deriving the first analytical model of the halo mass function (HMF). The subsequent work by Bond et al. (1991) established a new analytical approach to dark matter collapse known as the excursion set formalism, which reproduced the results from Press & Schechter but also provided a much wider framework.

However, given the highly nonlinear character of dark matter halos, N-body simulations have traditionally been the main arena in which to study their formation and abundance. Indeed, a long series of different fits for the HMF has arisen from simulations since the beginning of the century (Sheth et al. 2001; Jenkins 2001; Reed et al. 2003, 2007; Warren et al. 2006; Tinker et al. 2008; Crocce et al. 2010; Courtin et al. 2011; Watson et al. 2013; Despali et al. 2016; Diemer 2020), each of them for different mass and redshift ranges. Moreover, baryonic processes such as radiative cooling or feedback from active galactic nuclei are thought to cause non-negligible effects on the abundance of dark matter halos (Sawala et al. 2013; Bocquet et al. 2015; Castro et al. 2020). Consequently, work toward an observational determination of the abundance of dark matter halos may shed some light on the underlying physics and validate the results from simulations.

With regard to direct measurements of the HMF, it is well known that they present several challenges, such as the Eddington bias of mass-observable relations (Eddington 1913). Weak lensing-based stacked background shear distributions have thus far been a popular method for calibrating these quantities by means of linking the mean halo mass to observables such as luminosity or stellar mass (Simet et al. 2017; Leauthaud et al. 2017; Murata et al. 2018). Mass-observable relations can then be exploited to constrain the halo mass function from observations (Castro et al. 2016; Li et al. 2019). However, also within the context of gravitational lensing, we propose the angular foreground-background cross-correlation function as an observable of the weak lensing magnification bias effect and a method for shedding light on the direct measurement of the halo mass function.

Although traditionally considered suboptimal under usual circumstances by comparison with tangential shear measurements, we argue that magnification bias is a promising candidate for cosmological analysis when combined with a sample of background submillimeter galaxies (González-Nuevo et al. 2017; Bonavera et al. 2019, 2020). When a nonzero signal of this cross-correlation is found between a foreground and a background sample of massive objects with nonoverlapping redshift distributions, it signals the physical effect of magnification bias, by which an excess of background sources is measured around foreground sources with respect to what would be found in the absence of the latter (Bartelmann & Schneider 2001). Given its direct dependence on the halo mass of the foreground lenses, we propose this observable as a useful observational probe of the HMF under a halo model that only needs additional statistical information about the halo occupation distribution (HOD).

This paper constitutes a first step forward from the recent work of Cueli et al. (2021), where the cross-correlation between a foreground sample of GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014; Liske et al. 2015) sources and a foreground sample of H-ATLAS (Pilbratt et al. 2010; Eales et al. 2010) submillimeter galaxies was measured with the aim of providing constraints on the halo mass function. By performing a theoretical halo modeling and parametrizing the halo mass function according to the well-known Sheth and Tormen (ST) and Tinker fits, Bayesian posterior probability distributions for their parameters were obtained, and it was concluded that the former model can potentially be well constrained when given enough parameter freedom. This provided a method to obtain tabulated credible intervals of the halo mass function at the redshift of the lenses. However, the present work aims to obtain stronger constraints by performing a tomographic analysis (analogously to Bonavera et al. 2021) that divides foreground sources into four redshift bins in order to incorporate the likely time evolution of the galaxy content of dark matter halos. This approach can also provide information on the universality of the halo mass function by deriving posterior distributions for its parameters at different redshifts. Furthermore, and in light of the work by González-Nuevo et al. (2021), the cross-correlation function will be measured using two different tiling schemes, the first to keep coherence with Cueli et al. (2021) and the second because it provides more meaningful statistical conclusions.

The paper has been structured as follows. Section 2 sets the theoretical foundations of this work. We describe the usual framework to characterize the HMF and the parametrization we use as well as the halo model prescription for the angular cross-correlation function. Section 3 details the methodology we follow, namely a description of the data, the measurement process, the different tiling schemes, and the Markov chain Monte Carlo (MCMC) algorithm to sample the posterior distribution. Section 4 presents the results we obtained for both the non-tomographic and tomographic cases and for an additional analysis on universality. Our conclusions are summarized in Sect. 5.

2. Theoretical framework

2.1. The halo mass function

The (differential) halo mass function n(M, z) is defined as the comoving mean number density of halos at redshift z per unit mass. In other words, ${\int }_{M_1}^{M_2}n(M,z)\text{d}M$ is the comoving mean number density of halos of mass between M₁ and M₂ at redshift z.

A useful way of parametrizing this quantity is

$$\begin{array}{c}\hfill n(M,z)=\frac{{\rho }_0}{M^2}f(\nu ,z)|\frac{\mathrm{d}ln\nu (M,z)}{\mathrm{d}lnM}|,\end{array}$$(1)

where

$$\begin{array}{c}\hfill \nu (M,z)\equiv [\frac{{\widehat{\delta }}_{\mathrm{c}}(z)}{\sigma (M,z)}{]}^2,\end{array}$$

${\widehat{\delta }}_{\mathrm{c}}(z)$ being the linear critical overdensity at redshift z for a region to collapse into a halo at that same redshift and ρ₀ the current mean background density. It should be noted that the redshift dependence of ${\widehat{\delta }}_{\mathrm{c}}(z)$ is known to be very weak, but we decided to take it into account through the fit by Kitayama & Suto (1996) given its potential influence on the non-universality of the HMF (Courtin et al. 2011). Moreover, σ²(M, z)≡D²(z) σ²(M, 0), where D(z) is the linear growth factor for a ΛCDM universe (normalized at z = 0) and σ²(M, 0) is the mass variance of the filtered linear overdensity field at present. The function f(ν, z) depends on the HMF model we consider and, for the purposes of this work, we fix it to the ST fit, that is,

$$\begin{array}{c}\hfill f(\nu ,z)=A\sqrt{\frac{a\nu }{2\pi }}[1+(\frac{1}{a\nu }{)}^p]e^{-a\nu /2}.\end{array}$$(2)

Several comments should be made about this. Firstly, the choice of the ST fit over other models follows from the results of Cueli et al. (2021): for our study, this model seems to compare favorably with other fits, where an important degeneration among parameters is bound to arise together, thus compromising the statistical conclusions, especially when they depend on the prior range. Secondly, the original work by Sheth & Tormen (1999) imposed a normalization condition accounting for the fact that all mass should be bound up in halos, which allowed them to fix A in terms of p (and resulted in the commonly used best-fit triple A = 0.33, a = 0.707 and p = 0.3). Following the results from Cueli et al. (2021), we opted to leave the A parameter free in most cases, provided however that the condition ${\int }_0^{\infty }f(\nu )/\nu \text{d}\nu \le 1$ is met so as to have a meaningful physical interpretation of the HMF. Lastly, it should be noted that the original ST function f(ν, z) does not depend on redshift, but given the tomographic character of this work, we aim to study the possible time evolution of its parameters to quantify deviations from universality.

2.2. The foreground-background cross-correlation function

Aside from galaxy tangential shear measurements, the galaxy-mass correlation can be probed via the foreground-background cross-correlation function in the context of weak gravitational lensing. In the case of nonoverlapping redshift distributions, this quantity can be expressed as (Cooray & Sheth 2002)

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

where

$$\begin{array}{c}\hfill W^{\mathrm{lens}}(z)=\frac{3}{2}\frac{H_0^2}{c^2}[\frac{E(z)}{1+z}{]}^2{\displaystyle {\int }_z^{\infty }\mathrm{d}{z}^{\prime }\frac{\chi (z)\chi (z^{\prime }-z)}{\chi (z^{\prime })}{n}_{\mathrm{b}}(z^{\prime }),}\end{array}$$

n_b(z) (n_f(z)) being the unit-normalized background (foreground) source distribution, χ(z) the comoving distance at redshift z, J₀ the zeroth-order Bessel function of the first kind, and β the logarithmic slope of the background source number counts, commonly fixed to 3 for submillimeter galaxies (Lapi et al. 2011; Cai et al. 2013; Bianchini et al. 2015, 2016; González-Nuevo et al. 2017; Bonavera et al. 2019, 2020).

The nonlinear galaxy-dark matter cross-power spectrum is parametrized according to the halo model formalism as described in detail in Appendix A. Besides cosmology, it depends on the halo density profile, which we assume to match an Navarro-Frenk-White (NFW) model (Navarro et al. 1997), and the HOD parameters, which describe how galaxies populate halos on average. Indeed, according to the three-parameter model by Zehavi et al. (2005), the mean number of galaxies within a halo of mass M is given by

$$\begin{array}{c}\hfill {⟨N_{\mathrm{g}}⟩}_{\mathrm{M}}={⟨N_{c_{\mathrm{g}}}⟩}_{\mathrm{M}}+{⟨N_{c_{\mathrm{s}}}⟩}_{\mathrm{M}}=[1+(\frac{M}{M_1}{)}^{\alpha }]\mathrm{\Theta }(M-M_{\min }),\end{array}$$(4)

where ⟨N_cg⟩_M and ⟨N_sg⟩_M are the mean number of central and satellite galaxies in a halo of mass M, respectively. Parameter M_min describes the minimum mean mass for a halo to host a (central) galaxy, while M₁ represents the mean halo mass at which there is exactly one satellite galaxy. As described in Appendix A, we define halos as spherical overdense regions with an overdensity equal to the virial value, as estimated via the fit by Weinberg & Kamionkowski (2003) and a concentration parameter given by that of Bullock et al. (2001).

As a consequence, for a fixed cosmology, which we assume to match Planck’s (Planck Collaboration VI 2020), the angular cross-correlation function (3) depends on both the HMF and the HOD parameters. Our goal is thus to derive posterior probability distributions for these two sets of parameters.

3. Work methodology

3.1. Galaxy samples

The foreground and background sources have been extracted, respectively, from the GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014) and H-ATLAS (Pilbratt et al. 2010; Eales et al. 2010) surveys, which were coordinated to maximize the common area. Indeed, they both covered the three equatorial regions at 9, 12 and 14.5 h and part of the south galactic pole, amounting to a common area of ∼207 deg².

The foreground sample consists, in principle, of GAMA II sources within the interval 0.1 < z < 0.8, which results in ∼225 000 galaxies with a median spectroscopic redshift of z_med = 0.28. This is the same sample that was used in Cueli et al. (2021), but given the tomographic character of this work, it has been further divided into four redshift bins, namely, [0.1, 0.2], [0.2, 0.3], [0.3, 0.5] and [0.5, 0.8], referred to as bin 1 to bin 4. It should be noted that this particular choice was made to ensure a similar number of sources in each bin. Moreover, we have included additional sources with respect to our previous works (González-Nuevo et al. 2017; Bonavera et al. 2020; Cueli et al. 2021), which excluded galaxies with z < 0.2. The redshift distribution of our foreground sample is shown in Fig. 1 (in orange), where the four different bins are also depicted in different color shades. Given the spectroscopic nature of the measured redshifts, no source is assigned the wrong bin.

Fig. 1. Normalized redshift distribution of the background (in green) and foreground (in orange) samples of galaxies. The shaded regions represent the selected redshift ranges: 1.2 < z < 4.0 for the background (in green) and 0.1 < z < 0.2, 0.2 < z < 0.3, 0.3 < z < 0.5 and 0.5 < z < 0.8 for the four different foreground bins (in orange).

The background sample has been selected from the H-ATLAS sources detected in the aforementioned regions. In order to guarantee that there was no overlap with the foreground sample, a photometric redshift selection in the interval 1.2 < z < 4.0 has been applied, which results in ∼57 930 galaxies. The procedure of redshift estimation is extensively described in González-Nuevo et al. (2017) and Bonavera et al. (2019), but mainly consists in a χ² fit to the SPIRE (and PACS when available) data of the SMM J2135–0102 template spectral energy distribution (“the Cosmic Eyelash”; Ivison et al. 2010; Swinbank et al. 2010), which was found to be the best overall fit. The resulting redshift distribution of our background sample if shown in Fig. 1 (in green), where the random errors in photometric redshift estimation have been taken into account, as explained in González-Nuevo et al. (2017).

3.2. Measurement of the cross-correlation function

The H-ATLAS scanning strategy produced characteristic overlapping rhomboidal shapes or “tiles” in most of the fields. Considering the area of study, we have 16 such tiles, each with an area of ∼16 deg². As thoroughly described in González-Nuevo et al. (2021), the choice of different tiling schemes to measure the cross-correlation function (for example, dividing each tile into four “minitiles” in order to have more regions over which to average; see Fig. 2) brings along different large-scale biases in the measurements that need to be corrected so as to have a consistent method to estimate our observable. In principle, since the present work aims to extend that of Cueli et al. (2021), where the tile scheme was adopted, coherence would require the same approach be used. However, in view of the work by González-Nuevo et al. (2021), we opt to also consider the minitile scheme at first given the minimum surface area corrections that are required and, most importantly, the larger sample size, as will now be emphasized.

Fig. 2. Representation of the tiles (in red) and minitiles (in blue) schemes for the 9 h GAMA field. The pattern is analogous for the rest.

For each tiling scheme (tiles or minitiles) and redshift bin, we estimated the angular cross-correlation function via a modified version of the Landy & Szalay (1993) estimator (Herranz 2001),

$$\begin{array}{c}\hfill {\stackrel{\sim }{w}}^i(\theta )=\frac{D_{\mathrm{f}}^i{D}_{\mathrm{b}}(\theta )-D_{\mathrm{f}}^i{R}_{\mathrm{b}}(\theta )-D_{\mathrm{b}}{R}_{\mathrm{f}}^i(\theta )+R_{\mathrm{f}}^i{R}_{\mathrm{b}}(\theta )}{R_{\mathrm{f}}^i{R}_{\mathrm{b}}(\theta )},\end{array}$$

where the superscript i denotes the ith redshift bin and D_fD_b(θ), D_fR_b(θ), D_bR_f(θ) and R_fR_b(θ) denote the normalized foreground-background, foreground-random, background-random and random-random pair counts for the angular separation θ. The random catalog for each tile or minitile is generated from mock random positions for ten times the number of lenses and sources, respectively. As studied in detail by González-Nuevo et al. (2021), due to the relatively low surface density of the foreground sample and the size of the tiles or minitiles, there is no need of a more complex building procedure because the surface density variations can be considered negligible. In addition, the cross-correlation measurement for each tile or minitile is obtained by averaging over ten different random (Poisson) catalog realizations. In order to minimize cosmic variance, the final value for a given θ corresponds to the mean over tiles or minitiles, with an associated uncertainty of $\sigma /\sqrt{n}$, where σ is the standard deviation of the average measurements and n is the number of tiles or minitiles. Finally, to take into account the finite area of the tiles/minitiles, the integral constraints estimated by González-Nuevo et al. (2021), 5 × 10⁻⁴ and 9 × 10⁻⁴ respectively, are used to correct the estimated cross-correlation functions. This correction is, however, mainly relevant at the largest scales (> 20 arcmin).

The data obtained for each foreground redshift bin are shown in Fig. 3, where the red dots refer to the tiles scheme and the blue to the minitiles scheme. It should be noted that cross-correlation measurements performed with the tiles scheme generally have larger error bars. This is the expected behavior, given the fact that the number of regions over which to average is about four times as big in the minitiles case. In fact, some measurements at smaller scales show significant differences between both tiling schemes that were not present in the non-tomographic analysis by González-Nuevo et al. (2021). This is most likely due to the reduction in sample size due to splitting the foreground galaxies in several bins.

Fig. 3. Cross-correlation measurements at different θ values for the tiles (in red) and minitiles (in blue) schemes. The data for foreground bins 1 to 4 are depicted from left to right and from top to bottom.

It should be noted that, on top of the measurements for each individual redshift bin, ${\stackrel{\sim }{\mathit{w}}}^i(\theta )$, we performed (with the same approach) a single-bin computation of the cross-correlation function, $\stackrel{\sim }{\mathit{w}}(\theta )$, in order to compare the tomographic results with those of a non-tomographic analysis in the redshift bin [0.1 − 0.8].

3.3. Parameter estimation

In order to derive Bayesian posterior probability distributions for the involved parameters, we performed an MCMC analysis using the open source emcee software package (Foreman-Mackey et al. 2013), a Python-based implementation of the Goodman & Weare affine invariant MCMC ensemble sampler (Goodman & Weare 2010). An important comment should be made before going further: emcee is not an optimizer, but a sampler. As such, the maximum of the N-dimensional posterior distribution (that is, the unique tuple of parameters that best fits the data) we obtain from it is not a reliable estimate of the parameters, since there can be a considerable N-dimensional degeneracy. Due to this degeneracy, other tuples with relatively different values can provide as good a fit¹. This single tuple could be plotted to ensure the sampling was correctly carried out, but the full statistical picture is provided by the N-dimensional posterior distribution. When one derives the marginalized distribution for a single parameter, its mode, mean or median need not coincide with the corresponding value from the N-dimensional best fit (it can even result in a poor fit!). The bottom line is that point estimates should be treated with care.

We carried out different MCMC runs assuming that errors are Gaussian and independent and parametrizing the HMF with the ST fit, following the conclusions of Cueli et al. (2021). It is at this very point that the choice of a tiling scheme becomes relevant. If the sample size for a given θ value is not large enough, statistical conclusions may be compromised since the random variable w(θ) may significantly deviate from pure Gaussianity. As a consequence, both schemes will be analyzed in the non-tomographic case, since relevant differences are not expected, but only the results coming from averaging over minitiles will be further taken into account, since these provide the most rigorous statistical conclusions.

Firstly, we performed a non-tomographic (or single-bin) analysis with both tiling schemes. If {P_i}_i denotes the set of parameters to be estimated in the run, the log-likelihood function is expressed as

$$\begin{array}{cc}\hfill log\mathcal{L}({\theta }_1,\dots ,{\theta }_m;{\{P_i\}}_i)=& -\frac{1}{2}\sum _{j=1}^m[log2\pi {\sigma }_j^2\hfill \\ \hfill & +\frac{[w({\theta }_j;{\{P_i\}}_i)-\stackrel{\sim }{w}({\theta }_j)){]}^2}{{\sigma }_j^2}],\hfill \end{array}$$

where σ_j is the error in the jth measurement and w(θ_j) and $\stackrel{\sim }{\mathit{w}}({\theta }_j)$ are the theoretical and measured value of the cross-correlation function at angular scale θ_j. Moreover, {P_i}_i = {A, a, p, α, M_min, M₁}.

Secondly, we performed a tomographic run in all four redshift bins altogether, that is, we considered the joint likelihood for all bins. To analyze the effect that a time evolution of the HOD could have in the universal ST fit, the HMF parameters are not considered to vary among bins. The log-likelihood function thus reads²

$$\begin{array}{cc}\hfill log\mathcal{L}({\theta }_1,\dots ,{\theta }_m;{\{P_i\}}_i)=& -\frac{1}{2}\sum _{k=1}^4\sum _{j=1}^m[log2\pi {({\sigma }_j^k)}^2\hfill \\ \hfill & +\frac{[w^k({\theta }_j;{\{P_i\}}_i)-{\stackrel{\sim }{w}}^k({\theta }_j)){]}^2}{{({\sigma }_j^k)}^2}],\hfill \end{array}$$

where ${\sigma }_j^k$ is the error in the j_th measurement in redshift bin k and w_k(θ_j) and ${\stackrel{\sim }{\mathit{w}}}_k({\theta }_j)$ are the theoretical and measured value of the cross-correlation function in redshift bin k at angular scale θ_j. Moreover, {P_i}_i = {A, a, p, α₁, M_min1, M₁₁,α₂, M_min2, M₁₂, α₃, M_min3, M₁₃, α₄, M_min4, M₁₄}. It should be noted that the assumption of statistical independence of the cross-correlation measurement among different bins is justified: there is no possible bin overlap in redshift as the foreground redshifts are all spectroscopic and the probability of double lensing events or of an alignment of two lenses from different bins and a background galaxy is very low. Considering the typical lensing probability of order 10⁻³, the double lensing probability would be of order 10⁻⁶, thus completely negligible from the statistical point of view (see Bonavera et al. 2021, and references therein).

Lastly, we performed individual runs in each redshift bin separately in order to account for the possible redshift variation of the HMF parameters. The log-likelihood function for redshift bin k is then expressed as

$$\begin{array}{cc}\hfill log{\mathcal{L}}^k({\theta }_1,\dots ,{\theta }_m;{\{P_i^k\}}_i)=& -\frac{1}{2}\sum _{j=1}^m[log2\pi {({\sigma }_j^k)}^2\hfill \\ \hfill & +\frac{[w^k({\theta }_j;{\{P_i^k\}}_i)-{\stackrel{\sim }{w}}^k({\theta }_j)){]}^2}{{({\sigma }_j^k)}^2}],\hfill \end{array}$$

where $\left\{P_i^k\right\}=\{A_k,a_k,p_k,{\alpha }_k,M_{{\text{min}}_{\text{k}}},M_{1_k}\}$.

With regard to the prior probability distributions, they were chosen to be uniform with a wide enough range for all the involved parameters to find a reasonable balance between informational bias and computation efficiency, as shown in the tables throughout the paper. The delicate issue regarding the choice of the allowed portion of parameter space in an MCMC analysis was raised in Cueli et al. (2021) and the discussion therein led to the conclusion that the p parameter in the ST fit should be allowed to take nonpositive values.

Lastly, the full posterior distribution of the set of parameters in each case is sampled to derive observational constraints on the HMF for several mass values at any lens redshift. Indeed, although marginalized distributions are also obtained for all parameters, it is the N-dimensional probability distribution as a whole (which in our case will be summarized by a median value and credible intervals for the HMF at different halo masses) that provides the full picture and, not surprisingly, that proves useful for our purposes.

4. Results

4.1. Non-tomographic analysis

Table 1 shows the summarized statistical results from the non-tomographic MCMC run with the tiles scheme. For visual purposes, the part of the full corner plot showing only the HMF parameters is depicted in Fig. 4 (in red) and the corresponding marginalized posterior distribution of the HOD parameters is shown in Fig. 5 (in red). Figure B.1 (in red) depicts the full corner plot.

Fig. 4. One- and two-dimensional (contour) posterior distributions for the HMF parameters from the non-tomographic run with the tiles (in red) and minitiles (in blue) schemes.

Fig. 5. One-dimensional posterior distributions for the HOD parameters from the non-tomographic run using the tiles (in red) and minitiles (in blue) scheme.

As can be seen in Fig. 4 (in red), all three HMF marginalized distributions display clear peaks and confirm the findings from Cueli et al. (2021) with fixed HOD values: the submillimeter galaxy magnification bias constrains the ST parameters when p is allowed to take negative values and A is left free.

As expected, the uncertainties increase with respect to the fixed HOD case, but there is little to no variation in the p parameter and only a widening of the curve for a (the mode remaining practically unchanged). As for the normalization parameter, A, its distribution is displaced to smaller values because the freedom in the HOD allows it to; for instance, as shown by Fig. B.1, fixed values of 12.4 and 13.6 for log M_min and log M₁ imply lower A values. More precisely, the mean and 68% credible intervals are $A=0.{21}_{-0.18}^{+0.07}$, $a=3.{93}_{-3.08}^{+1.23}$, and $p=-0.{94}_{-1.54}^{+1.13}$, with median values of 0.18, 3.37, and −0.83, respectively. It should be noted that, as in the analysis of Cueli et al. (2021) with fixed HOD parameters, p clearly prefers negative values as opposed to the standard 0.3 from N-body simulations. Moreover, and since the a distribution sources the majority of the uncertainty in the HMF, these results already hint not very restrictive constraints on this quantity, as will be revealed later.

As for the marginalized distribution of the HOD parameters, Fig. 5 shows a well-constrained Gaussian-like posterior for the minimum mean halo mass, $log{M}_{\min }=12.{18}_{-0.22}^{+0.91}$ (with a median of 12.51, whereas low values for M₁ seem to be slightly preferred without ruling out the highest masses due to the restriction that M_min should be less than M₁: $log{M}_1=11.{17}_{-2.17}^{+0.55}$ with a median of 10.90. As already expected given the results from previous works, the α parameter is unconstrained, given its little influence on our observable. The 68% full posterior sampling of the cross-correlation is shown in Fig. 6 (top panel), along with the data corresponding to the tiles scheme, confirming that the MCMC algorithm was correctly carried out.

Fig. 6. Posterior-sampled cross-correlation function (68% band) and measured data for the non-tomographic case with the tiles scheme (top panel) and the minitiles scheme (bottom panel).

Table 2 shows the corresponding statistical results for the non-tomographic run with the minitiles scheme. As in the previous case, the part of the full corner plot showing only the HMF parameters and the marginalized posterior distributions of the HOD parameters are depicted in Fig. 4 (in blue) and in Fig. 5 (in blue), respectively. Figure B.1 (in blue) depicts the full corner plot.

From a qualitative point of view, Figs. 4 and 5 show a similar behavior of all six parameters with respect to the previous case using the tiles scheme. However, as clearly revealed by Fig. B.1, there are some slight changes. As regards the HMF parameters, all three marginalized distributions are somewhat displaced toward lower values and those of A and a are narrower; indeed, the mean and 68% credible intervals are $A=0.{18}_{0.16}^{+0.06}$, $a=3.{34}_{-2.66}^{+0.94}$ and $p=-1.{09}_{-0.65}^{+1.23}$, with median values of 0.16, 2.85 and −0.99, respectively.

With respect to the HOD parameters, and in comparison to the previous case, the main difference is that higher values of M₁ are more strongly suppressed, with an 68% and 95% upper bound of 11.72 and 14.26 for log M₁. As for M_min, it barely changes with respect to using the tiles scheme, with a mean of $log{M}_{\min }=12.{02}_{-0.30}^{+1.01}$ The α parameter is again unconstrained, although can be assigned a highest posterior density 68% lower bound of 0.68. The corresponding 68% posterior sampling of the cross-correlation function is shown in Fig. 6 (bottom panel), along with the data corresponding to the minitiles scheme, confirming once more that the MCMC algorithm was correctly carried out.

Moreover, and most importantly for the purpose of this paper (the observational constraining of the number density of halos), credible intervals for the HMF at any lens redshift can be obtained via the sampling of the full posterior distribution. Figure 7 shows the median values and the 68% and 95% credible intervals at several masses for z = 0.4 in the non-tomographic case for both the tiles and minitiles scheme (in red and blue, respectively), along with the ST best fit found by Despali et al. (2016) under a Planck cosmology (green line). One should keep in mind that the algorithm used by Sheth & Tormen (1999) identified halos as spherically overdense regions with an overdensity equal to 178 times the mean density; our mass definition is thus different since it relies on the virial value, and this is why direct comparisons should be performed under the same halo mass definition, hence the aforementioned choice. While in agreement with the usual values (and between themselves) within the uncertainties, the error bars are relatively sizeable, specially in the tiles scheme and for large masses. It should be noted that the error bars from Cueli et al. (2021) (where the tiles scheme was used) were substantially smaller due to fixing the HOD distribution for the final conclusions. No significant differences are found between both approaches in this non-tomographic case (aside from a slight tendency toward larger values of the HMF and the smaller overall uncertainties for the minitiles scheme), which justifies the use of the tiles scheme in the preliminary work by Cueli et al. (2021), which meant to provide a proof-of-concept method to observationally constrain the HMF.

Fig. 7. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in the non-tomographic setup within the minitiles (blue) and tiles (red) schemes. The ST best fit by Despali et al. (2016) is drawn with a green line.

4.2. Tomographic analysis

Having explored the non-tomographic case, we now turn to the main objective of this work, which is to study the potential improvements that can arise from a tomographic analysis. To do so, and as explained in Sect. 3.3, we assume that all HMF parameters do not vary across redshift bins, that is, that universality of the HMF holds, while we let the HOD evolve. As stressed in Sect. 3.3, only the minitiles scheme will be considered from now on.

Table 3 shows the summarized statistical results from this tomographic run, compared with with the results from the non-tomographic run (in blue). The corresponding part of the full corner plot regarding only the HMF parameters is shown in Fig. 8 (in blue) and the marginalized posterior distributions of the HOD are depicted in Fig. 9 (with a different black line for each redshift bin). Figure B.2 (in black) shows the full corner plot compared with the results from the non-tomographic run (in blue).

Fig. 8. One- and two-dimensional (contour) posterior distributions for the HMF parameters from the tomographic (in black) and non-tomographc (in blue) runs in the minitiles scheme.

The improvement with respect to the non-tomographic case is clear-cut. Indeed, as shown by Fig. 8, the marginalized posterior distributions of a and p are significantly narrowed when compared to the analysis with a single redshift bin. More specifically, $a=2.{51}_{-0.75}^{+0.73}$ and $p=-1.{05}_{-0.59}^{+0.44}$, with median values of 2.45 and −1.09, respectively. It should be mentioned that the non-negligible probability region at around a < 1 (also visible in the blobs on the A − a plane) is not due to a lack of convergence of the MCMC chains, but to a small number of walkers that remain near those values, as has been confirmed by analyzing their path in parameter space³. In fact, these values are in the neighborhood of the traditional ST triple. Concerning the marginalized distribution of the A parameter, it is practically neither narrowed nor widened with respect to the non-tomographic case, but clearly moves toward larger values: $A=0.{40}_{-0.07}^{+0.12}$. The overall behavior of the HMF parameters can already hint at universality as a very valid assumption for our redshifts, given that the new a and p distributions emphasize the goodness of the fit, but, most importantly, this run clearly shows a substantial improvement regarding uncertainties.

As for the behavior of the HOD parameters with respect to redshift, Fig. 9 shows a clear redshift evolution of the mean minimum halo mass, with larger values of M_min corresponding to earlier times. Particularly, the mean and 68% credible intervals for log M_min are $10.{83}_{-0.39}^{+0.86}$, $11.{71}_{-0.17}^{+1.00}$, $12.{42}_{-0.25}^{+0.52}$ and $13.{49}_{-0.04}^{+0.40}$ for bins 1 to 4, respectively. This is the expected behavior (González-Nuevo et al. 2017; Bonavera et al. 2021): given that we are working with a flux-limited sample, there is an observational bias toward intrinsically brighter (and thus more massive) objects. With regard to M₁, the results for bins 1 and 2 do not show the corresponding evolution, while those from bins 3 and 4 do, with mean values for log M₁ of $12.{91}_{-0.72}^{+0.49}$ and $14.{60}_{-0.40}^{+1.40}$, respectively. As for the α parameter, it can only be assigned a proper lower bound in bin 3, being > 1.25 at 68%, while its overall redshift behavior remains unclear.

Fig. 9. One-dimensional posterior distributions for the HOD parameters in the tomographic case under the assumption of universality (in black) and in the non-tomographic case (in blue).

Figure 10 shows the corresponding median values and credible intervals for the z = 0.4 HMF and brings to light the remarkable improvements in the actual observational constraints of the HMF. As expected, the uncertainties are substantially smaller with respect to the non-tomographic case, even for 95% credibility. The error bars for small masses (below ∼10¹² M_⊙ h⁻¹) show that our results are not compatible with the best-fit ST values found by Despali et al. (2016) at 95% credibility. In fact, under the assumption that the HMF is well fit by the ST model, our results predict a larger number density of low-mass halos and (although only at 1σ) a smaller number of massive ones, that is, a steeper fall with a cutoff at ∼10¹⁴ M_⊙ h⁻¹.

Fig. 10. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in both the non-tomographic (blue) and tomographic (black) cases for the minitiles scheme. The ST best fit by Despali et al. (2016) is drawn with a green line.

At this point, it is of interest to analyze the potential of a two-parameter ST fit in the tomographic frame, that is, an ST fit where the A parameter is fixed according to the normalization condition mentioned in Sect. 2.1. As discussed in Cueli et al. (2021), this case is not suitable in a non-tomographic setting due to important degeneracies, a conclusion that has been further confirmed for this work with several tests that yielded nonrestrictive constraints on the a parameter. However, it might be the case that this problem is solved via tomography. Table 4 shows the summarized statistical results from this tomographic run. The corresponding part of the full corner plot regarding only the HMF parameters is shown in Fig. 11. Figure B.3 depicts the full corner plot from this run, comparing the results with those from the three-parameter case.

Fig. 11. One- and two-dimensional (contour) posterior distributions for the HMF parameters in the tomographic setup for the three-parameter (in black) and two-parameter (in purple) fits.

The most remarkable result is the fact that constraints on the HMF parameters are even tighter with respect to the three-parameter case and, surprisingly, that the difficulties present in the two-parameter non-tomographic case are surmounted by the adoption of the tomographic frame. In particular we estimate the following values for the HMF parameters: $a=2.{47}_{0.46}^{+0.43}$ and $p=-1.{29}_{-0.37}^{+0.36}$, with median values of 2.45 and −1.30, respectively. The traditional ST values are rejected at more than 3σ, as will be later discussed.

As for the HOD parameters, Fig. B.3 depicts a qualitatively similar situation to the three-parameter case except for a large discrepancy in the third bin. Indeed, the M_min parameter shows the expected redshift evolution in bins 1, 2 and 4, with mean values for log M_min of $10.{75}_{-0.60}^{+1.08}$, $11.{38}_{-0.59}^{+1.27}$ and $13.{12}_{-0.06}^{+0.97}$, respectively. However, the one-dimensional distribution for bin 3 displays an unusually large probability region for low masses, which breaks this tendency. We believe this effect is caused by the seemingly preferred low values of M₁ in this bin, which in turn forces M_min to be smaller. Regardless of this feature, it should be mentioned that the marginalized posterior distributions of both a and p converge much more quickly than those of the HOD parameters, which proves the robustness of these results with respect to the HOD. Moreover, the main goal of this work is to provide observational constraints on the HMF, not the precise determination of the HOD.

Figure 12 shows the median values and credible intervals for the z = 0.4 HMF obtained for this run (in purple) and compared to those from the three-parameter fit (in black) and the ST best fit found by Despali et al. (2016) (in green). As commented before and stressed by this figure, the disagreement between our results for the two-parameter case and the best-fit ST values is stronger than in the previous scenario, even at large masses and it can safely be said that we predict a higher number density of halos for masses below ∼10¹³ M_⊙ h⁻¹ and lower number density for masses above ∼10¹⁴ M_⊙ h⁻¹, and thus a steeper fall in mass.

Fig. 12. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in the tomographic setup for the three-parameter (in black) and the two-parameter (in purple) fits. The ST best fit by Despali et al. (2016) is drawn with a green line.

4.3. A step toward testing universality

Since we have measured the cross-correlation function for each foreground bin, we can now explore the possibility of allowing the HMF parameters to vary with redshift and perform four different non-tomographic MCMC runs, one for each bin. It should be noticed that it is well established that universality in the HMF strongly depends on the mass definition of halos (Courtin et al. 2011; Despali et al. 2016; Ondaro-Mallea et al. 2021), which means that our current results in this area should be regarded as very preliminary, since a serious study concerning differences with regard to mass definition is beyond the scope of this work.

Taking into account that the two-parameter fit does not perform well in a non-tomographic setting, we carried out these runs for the three-parameter case only. Table 5 shows the summarized statistical results from these runs. For visual purposes, the marginalized posterior distributions of the HOD and the HMF parameters in each redshift bin are shown in Fig. 14, while Fig. B.4 depicts the full corner plot.

If we focus first on the HOD parameters, their marginalized posterior distributions surprisingly show the same distinctive feature as the two-parameter fit in the tomographic case. Indeed, both the M_min and M₁ parameters take unusually low values with very high probability in bin 3. The reason why this happens in this three-parameter case is not clear as of now (since the A₃ parameter is given freedom as opposed to the two-parameter scenario), but we stress that the main goal of this work is to constrain the HMF and not give tight bounds on the HOD; nevertheless, this issue will be studied in the future. As for the rest of bins, the behavior is completely similar to that of the tomographic case, with quantitative differences mainly of the form of a larger dispersion for bins 1 and 2. The situation is more interesting for the HMF parameters. With respect to p, Fig. 14 shows a perfect non-evolving probability distribution, with barely any difference among bins. As for A, the situation differs due again to bin 3. While the posterior distributions for the rest bins 1, 2, and 4 remain practically unchanged, that of bin 3 is displaced to notoriously larger values with a low-probability tail at the mode of the other distributions. This is again indicative of a characteristic behavior of bin 3, which however translates into a smaller spread in the distribution (and consequently, in the HMF). This is patently obvious in the marginalized distributions of the a parameter. Bin 1 displays an ever-increasing probability distribution, a behavior utterly opposed to that of the rest of bins, probably due to the poor cross-correlation data at this redshift. As for the a parameter in bin 3, it displays a relatively narrow probability distribution when compared to the rest, which keep the same tendency but with larger uncertainties.

Fig. 13. Posterior-sampled f(ν, z) function at different redshift bins. The 68% probability contours are depicted in red, cyan, orange, and black for bins 1 to 4, respectively. The ST best fit by Despali et al. (2016) is plotted with a black dotted line.

Fig. 14. One-dimensional posterior distributions for the HMF (top panel) and HOD (bottom panel) parameters in the last four nontomographic runs. The results from bins 1, 2, 3, and 4 are shown with a red, cyan, orange and black line, respectively.

Given the very large overlap in the distributions, universality seems to hold at first glance. However, to test it in a quantitative manner, we sampled the function f(ν, z) from (2) using the posterior distributions from all redshift bins, whose corresponding 68% probability bands are shown in Fig. 13, along with the ST best fit found by Despali et al. (2016) (black dotted line). The most obvious aspect to be noted is that the contour from bin 3 is extremely narrow with respect to the rest, which are patently wider due to the long high-probability tails of the marginalized distribution of the a parameter, which sources most of the HMF uncertainty. This causes that these three contours overlap at every value of ν, but those of bin 3 and bin 1 do not at 68% probability. However, as stressed at the start of this section and given the large uncertainties, this should be more thoroughly studied in the future by varying the halo mass definition and performing parallel tomographic analysis in different sets of bins.

5. Summary and conclusions

This paper has addressed the issue of observationally restricting the dark matter halo mass function by measuring the cross-correlation function between samples of GAMA II (with spectroscopic redshifts between 0.1 < z < 0.8) and H-ATLAS galaxies (with photometric redshifts between 1.2 < z < 4.0) in a tomographic analysis. Due to their physical properties, submillimeter galaxies have been chosen as background sources so as to optimize the measurement of the signal, which is a clear manifestation of the phenomenon of magnification bias. By means of a halo model interpretation of the cross-correlation function and the ST parametrization of the HMF, a Bayesian analysis through an MCMC algorithm has been carried out to obtain posterior probability distributions for both the HMF and HOD parameters under a fixed Planck cosmology and credible intervals for the number density of the dark matter halos associated with the lenses at given mass values. Two different tiling schemes, namely the so-called tiles and minitiles schemes, have been initially used for the sake of coherence, although the statistically rigorous results come from using the latter given the larger sample size.

We began the analysis with a non-tomographic study that considered the overall lens sample so as to assess the information that can be extracted without taking into account the possible redshift variation of the HOD parameters. The results from both tiling schemes are qualitatively similar, with slightly smaller overall uncertainties for the minitile case. The observational constraints on the HMF agree very well with the traditional result from numerical simulations, but the error bars are large, specially for halo masses above ∼10¹⁴ M_⊙ h⁻¹.

In order to explore the potential improvements provided by dividing the foreground sample, we performed a tomographic analysis under the assumption of universality of the HMF and the three-parameter ST fit in the minitiles scheme. The main result of this paper is the remarkable improvement in the observational constraints of the HMF when adopting a tomographic setup, which clearly confirms the forecast by Cueli et al. (2021). Under the universal three-parameter ST fit, and with respect to the traditional values from numerical simulations, we predict a higher number density of dark matter halos at masses below ∼10¹² M_⊙ h⁻¹ (at 95% credibility) and, although at only 1σ, a tendency toward a lower number density at masses larger than 10¹⁴ M_⊙ h⁻¹. The ST values found by Despali et al. (2016) are thus not compatible at 95% credibility.

Although a two-parameter ST fit did not perform well in a non-tomographic setup (as shown in Cueli et al. 2021), we decided to test it in the tomographic case. The results are even more remarkable and further restrict the HMF constraints from the previous case, predicting once more a higher number density of halos for low masses (in this case below ∼10¹³ M_⊙ h⁻¹ and a lower number density above ∼10¹⁴ M_⊙ h⁻¹, all this at more than 3σ credibility, thus disagreeing again the traditional ST triples from numerical simulations, this time more strongly. This is worth noticing and a comment should be made. On the one hand, the inconsistency of traditional ST parameter values with the data may challenge the interpretation of the ST formula in terms of the excursion set theory with an ellipsoidal barrier, since the original values are strictly related to the barrier shape and to the way the mass function is computed in such a theoretical framework. On the other hand, if the ST formula is taken at face value as a simple fit to results from N-body simulations, one might not consider this to be too relevant, especially given the compatibility at 3σ of the three-parameter case.

Lastly, and so as to use the cross-correlation measurements in each bin to test the assumption of universality in a very preliminary test, four independent MCMC runs were carried out. While the results for the HOD parameters are similar to those from the tomographic case (except for some characteristic feature in the third redsfhit bin), there appear to be discrepancies between the lenses at redshifts 0.3 − 0.5 and the rest, although given the large uncertainties, a more in-depth study regarding halo mass definition should be performed in the future to provide a more thorough analysis and robust conclusions.

Following this line of thought, an enlargement of both the foreground and background samples is needed to improve statistics (as discussed in Cueli et al. 2021) due to the nature of a tomographic analysis, where the number of sources used to compute the cross-correlation is decreased with respect to the non-tomographic case. Another interesting approach would involve parametrizing the HMF with a different model to study variations in our observational constraints. While Cueli et al. (2021) already performed such an analysis in the non-tomographic case comparing the ST and Tinker (Tinker et al. 2008) fits and concluding that statistical problems arose regarding prior ranges for the latter, it could very well be that these issues would be solved through tomography. Another possibility could of course be the use of other HMF models derived from alternatives to the standard excursion set formalism, like the stochastic theory by Lapi & Danese (2020, 2021).

Lastly, it should be noted that the halo model used in this work has its shortcomings and appears not to show the accuracy required for very precise measurements regarding cosmology, mainly due to the transition between the 1-halo and 2-halo regimes and/or baryonic physics (Mead et al. 2015, 2020; Philcox et al. 2020; Acuto et al. 2021). Future work will also address the possible improvements that could arise when the model is refined.

Acknowledgments

M.M.C., J.G.N., L.B., D.C., J.M.C. acknowledge the PGC 2018 project PGC2018-101948-B-I00 (MICINN/FEDER). M.M.C. acknowledges PAPI-20-PF-23 and PAPI-21-PF-04 (Universidad de Oviedo). A.L. is supported by the EU H2020-MSCA-ITN-2019 Project 860744 “BiD4BESt: Big Data applications for black hole Evolution STudies” and by the PRIN MUR 2017 prot. 20173ML3WW “Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes”. We deeply acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. In particular the COSMOGAL projects “SIS20_lapi”, “SIS21_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. 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: Ingredients of the model

This section is aimed at giving details about the computation of the theoretical model used in this work. The cross-power spectrum present in the angular cross-correlation function of (3) has been calculated via the halo model prescription, that is,

$$\begin{array}{c}\hfill P_{\mathrm{g}\text{-}\mathrm{dm}}(k,z)=P_{\mathrm{g}\text{-}\mathrm{dm}}^{1\mathrm{h}}(k,z)+P_{\mathrm{g}\text{-}\mathrm{dm}}^{2\mathrm{h}}(k,z),\end{array}$$

where

$$\begin{array}{c}\hfill P_{\mathrm{g}\text{-}\mathrm{dm}}^{1\mathrm{h}}(k,z)={\displaystyle {\int }_0^{\infty }dMM\frac{n(M,z)}{\overline{\rho }(z)}\frac{{⟨N_g⟩}_M}{{\overline{n}}_g(z)}|u_{\mathrm{dm}}(k,z|M)||u_{\mathrm{g}}{(k,z|M)|}^{p-1}}\end{array}$$

and

$$\begin{array}{cc}\hfill P_{\mathrm{g}\text{-}\mathrm{dm}}^{2\mathrm{h}}(k,z)=P(k,z)[& {\displaystyle {\int }_0^{\infty }dMM\frac{n(M,z)}{\overline{\rho }(z)}{b}_1(M,z)u_{\mathrm{dm}}(k,z|M)]}\hfill \\ \hfill ·[& {\displaystyle {\int }_0^{\infty }dMn(M,z)b_1(M,z)\frac{{⟨N_g⟩}_M}{{\overline{n}}_g(z)}{u}_g(k,z|M)]}\hfill \end{array}$$

are the 1-halo and 2-halo term of the cross-power spectrum, respectively. In these expressions, $\overline{\rho }(z)$ is the mean matter density of the Universe at redshift z, n(M, z) is the HMF, computed through (1) under the ST parametrization, b₁(M, z) is the determinsitic bias, computed through the peak background split as in Sheth & Tormen (1999). The mean number of galaxies in a halo of mass M, ⟨N_g⟩_M is parametrized as described in (4) and the mean number density of galaxies at redshit z is thus given by

$$\begin{array}{c}\hfill {\overline{n}}_g(z)={\displaystyle {\int }_0^{M}dMn(M,z){⟨N_g⟩}_M.}\end{array}$$

The linear matter power spectrum P(k, z) is evolved to redshift z through the linear growth factor approximation of Carroll et al. (1992) and the corresponding transfer function is computed via Eisenstein and Hu’s fitting formula Eisenstein & Hu (1998), which takes baryonic effects into account. We have chosen an analytic approach as opposed to a numerical computation using a Boltzmann solver mainly due to computation time, since no significant differences were found between the two.

The normalized Fourier transform of the dark matter distribution within an NFW halo of mass M is given by (Cooray & Sheth 2002)

$$\begin{array}{cc}\hfill u_{\mathrm{dm}}(k,z|M)& =\frac{4\pi {\rho }_s{r}_s^3(M,z)}{M}[sink{r}_s[\mathrm{Si}([1+c]k{r}_s)-\mathrm{Si}(k{r}_s)]\hfill \\ \hfill & -\frac{sinck{r}_s}{[1+c]k{r}_s}+cosk{r}_s[\mathrm{Ci}([1+c]k{r}_s)-\mathrm{Ci}(k{r}_s)]],\hfill \end{array}$$

where

$$\begin{array}{c}\hfill r_s(M,z)\equiv \frac{R_{\mathrm{vir}}}{c(M,z)}\end{array}$$(A.1)

and ρ_s are a scale radius and density that parametrize the profile, concentration parameter of a halo of mass M at redshift z, which satisfies

$$\begin{array}{c}\hfill M=4\pi {\rho }_s{r}_s^3[ln[1+c(M,z)]-\frac{c(M,z)}{1+c(M,z)}]\end{array}$$(A.2)

for an NFW profile. The virial radius R_vir was computed via the virial overdensity at redshift z, through the fit by Weinberg & Kamionkowski (2003). It should be noted that Sheth & Tormen (1999) identified halos using a spherical overdensity halo finder, but with an overdensity threshold of 178 times the background density. In practice, for a halo of mass M, we adopted the concentration parameter by Bullock et al. (2001), computed r_s through (A.1) and, subsequently, computed ρ_s using (A.2).

Lastly, two additional comments should be made. Firstly, it is a reasonable approximation (Sheth & Diaferio 2001) to set the Fourier transform of the galaxy distribution, u_g(k, z) equal to that of dark matter. Secondly, the exponent p should be set to 1 for central galaxies and to 2 for satellites.

Appendix B: Additional plots

Fig. B.1. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the non-tomographic with the tiles (in red) and minitiles (in blue) schemes.

Fig. B.2. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the tomographic run (in black) and the non-tomographic run (in blue) for the minitiles scheme. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.

Fig. B.3. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the tomographic runs in the minitiles scheme for the three-parameter (in black) and the two-parameter (in purple) cases. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.

Fig. B.4. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the last four non-tomographic runs. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.

All Tables

All Figures

	Fig. 1. Normalized redshift distribution of the background (in green) and foreground (in orange) samples of galaxies. The shaded regions represent the selected redshift ranges: 1.2 < z < 4.0 for the background (in green) and 0.1 < z < 0.2, 0.2 < z < 0.3, 0.3 < z < 0.5 and 0.5 < z < 0.8 for the four different foreground bins (in orange).
In the text

	Fig. 2. Representation of the tiles (in red) and minitiles (in blue) schemes for the 9 h GAMA field. The pattern is analogous for the rest.
In the text

	Fig. 3. Cross-correlation measurements at different θ values for the tiles (in red) and minitiles (in blue) schemes. The data for foreground bins 1 to 4 are depicted from left to right and from top to bottom.
In the text

	Fig. 4. One- and two-dimensional (contour) posterior distributions for the HMF parameters from the non-tomographic run with the tiles (in red) and minitiles (in blue) schemes.
In the text

	Fig. 5. One-dimensional posterior distributions for the HOD parameters from the non-tomographic run using the tiles (in red) and minitiles (in blue) scheme.
In the text

	Fig. 6. Posterior-sampled cross-correlation function (68% band) and measured data for the non-tomographic case with the tiles scheme (top panel) and the minitiles scheme (bottom panel).
In the text

	Fig. 7. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in the non-tomographic setup within the minitiles (blue) and tiles (red) schemes. The ST best fit by Despali et al. (2016) is drawn with a green line.
In the text

	Fig. 8. One- and two-dimensional (contour) posterior distributions for the HMF parameters from the tomographic (in black) and non-tomographc (in blue) runs in the minitiles scheme.
In the text

	Fig. 9. One-dimensional posterior distributions for the HOD parameters in the tomographic case under the assumption of universality (in black) and in the non-tomographic case (in blue).
In the text

	Fig. 10. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in both the non-tomographic (blue) and tomographic (black) cases for the minitiles scheme. The ST best fit by Despali et al. (2016) is drawn with a green line.
In the text

	Fig. 11. One- and two-dimensional (contour) posterior distributions for the HMF parameters in the tomographic setup for the three-parameter (in black) and two-parameter (in purple) fits.
In the text

	Fig. 12. Credible intervals (68% in bold and 95% in faint colors) for the z = 0.4 HMF at different mass values when the full posterior distribution is sampled in the tomographic setup for the three-parameter (in black) and the two-parameter (in purple) fits. The ST best fit by Despali et al. (2016) is drawn with a green line.
In the text

	Fig. 13. Posterior-sampled f(ν, z) function at different redshift bins. The 68% probability contours are depicted in red, cyan, orange, and black for bins 1 to 4, respectively. The ST best fit by Despali et al. (2016) is plotted with a black dotted line.
In the text

	Fig. 14. One-dimensional posterior distributions for the HMF (top panel) and HOD (bottom panel) parameters in the last four nontomographic runs. The results from bins 1, 2, 3, and 4 are shown with a red, cyan, orange and black line, respectively.
In the text

	Fig. B.1. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the non-tomographic with the tiles (in red) and minitiles (in blue) schemes.
In the text

	Fig. B.2. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the tomographic run (in black) and the non-tomographic run (in blue) for the minitiles scheme. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.
In the text

	Fig. B.3. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the tomographic runs in the minitiles scheme for the three-parameter (in black) and the two-parameter (in purple) cases. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.
In the text

	Fig. B.4. One- and two-dimensional (contour) posterior distributions for the HMF and HOD parameters from the last four non-tomographic runs. The results from bins 1 to 4 are shown, respectively, from left to right and from top to bottom.
In the text