Issue 
A&A
Volume 665, September 2022



Article Number  A38  
Number of page(s)  35  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202243711  
Published online  08 September 2022 
KiDS+VIKING+GAMA: Halo occupation distributions and correlations of satellite numbers with a new halo model of the galaxymatter bispectrum for galaxygalaxygalaxy lensing
^{1}
ArgelanderInstitut für Astronomie, Rheinische FriedrichWilhems Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
email: llinke@astro.unibonn.de
^{2}
MaxPlanckInstitut für extraterrestriche Physik, Giessenbachstrasse 1, 85748 Garching, Germany
^{3}
UniversitätsSternwarte, Fakultät für Physik, LudwigMaximiliansUniversität München, Scheinerstr. 1, 81679 München, Germany
^{4}
Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany
Received:
5
April
2022
Accepted:
28
June
2022
Context. Halo models and halo occupation distributions (HODs) are important tools to model the distribution of galaxies and matter.
Aims. We present and assess a new method for constraining the parameters of HODs using the mean gravitational lensing shear around galaxy pairs, socalled galaxygalaxygalaxy lensing (G3L). In contrast to galaxygalaxy lensing, G3L is also sensitive to the correlations between the perhalo numbers of galaxies from different populations. We employed our G3L halo model to probe these correlations and test the default hypothesis that they are negligible.
Methods. We derived a halo model for G3L and validated it with realistic mock data from the Millennium Simulation and a semianalytic galaxy model. Then, we analysed public data from the KiloDegree Survey (KiDS), the VISTA Infrared KiloDegree Galaxy Survey (VIKING) and data from the Galaxy And Mass Assembly Survey (GAMA) to infer the HODs of galaxies at z < 0.5 in five different stellar mass bins between 10^{8.5}h^{−2} M_{⊙} and 10^{11.5}h^{−2} M_{⊙} and two colours (red and blue), as well as correlations between satellite numbers.
Results. The analysis accurately recovers the true HODs in the simulated data for all galaxy samples within the 68% credibility range. The model best fits agree with the observed G3L signal on the 95% confidence level. The inferred HODs vary significantly with colour and stellar mass. In particular, red galaxies prefer more massive halos ≳10^{12} M_{⊙}, while blue galaxies are present in halos ≳10^{11} M_{⊙}. There is strong evidence (> 3σ) for a high correlation, increasing with halo mass, between the numbers of red and blue satellites and between galaxies with stellar masses below 10^{10} M_{⊙}.
Conclusions. Our G3L halo model accurately constrains galaxy HODs for lensing surveys of up to 10^{3} deg^{2} and redshift below 0.5 probed here. Analyses of future surveys may need to include nonPoisson variances of satellite numbers or a revised model for central galaxies. Correlations between satellite numbers are ubiquitous between various galaxy samples and are relevant for halos with masses ≳10^{13} M_{⊙}, that is, of galaxygroup scale and more massive. Possible causes of these correlations are the selection of similar galaxies in different samples, the survey flux limit, or physical mechanisms such as a fixed ratio between the satellite numbers of distinct populations. The decorrelation for halos with smaller masses is probably an effect of shot noise by lowoccupancy halos. The inferred HODs can be used to complement galaxygalaxy lensing or galaxyclustering HOD studies or as input to cosmological analyses and improved mock galaxy catalogues.
Key words: gravitational lensing: weak / galaxies: halos / cosmology: observations / largescale structure of Universe
© L. Linke et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Accurate models of the distribution of galaxies inside the cosmic largescale structure are crucial to understanding the physics of galaxy evolution and inferring cosmological parameters from galaxy surveys. Popular frameworks for analytically describing the galaxy and matter distribution are halo models (e.g. Cooray & Sheth 2002; Scoccimarro et al. 2001; Kravtsov et al. 2004; Zheng et al. 2005). Their key ingredient is the halo occupation distribution (HOD), which gives the expected number of galaxies inside a halo of a given mass (Berlind & Weinberg 2002). Here we present a new method to accurately infer galaxy HODs with a halo model for galaxygalaxygalaxy lensing (G3L) – the mean gravitational lensing shear around galaxy pairs (Schneider & Watts 2005). We demonstrate that with this higherorder statistic, we can obtain the HODs of various galaxy samples selected, for example, by their colour, in current galaxy surveys. Additionally, in contrast to other lensing methods, G3L can probe the correlation of perhalo numbers of satellite galaxies from different populations. In a first application, we infer the HODs and correlations of various galaxy selections in the KiloDegree Survey (KiDS; Kuijken et al. 2015), the VISTA Infrared Kilodegree Galaxy survey (VIKING; Edge et al. 2013; Venemans et al. 2015), and the Galaxy And Mass Assembly survey (GAMA; Driver et al. 2009; Liske et al. 2015).
Halo models are the backbone of many analytical expressions for the statistics of the largescale structure. They postulate that matter is distributed in virialised halos, and galaxies only exist within these halos (White & Rees 1978). For known halo density profiles and HODs, they can predict all statistics of the galaxy and matter distributions (Cooray & Sheth 2002). Therefore, with halo models, we can infer HODs from the measured galaxy and matter statistics, such as galaxy clustering (Zehavi et al. 2011; Ishikawa et al. 2021) and the galaxymatter power spectrum (Mandelbaum et al. 2006; Clampitt et al. 2016; Dvornik et al. 2018).
Since halo models are analytic, calculating their predictions is faster and easier than more complex approaches such as semianalytic galaxy models (SAMs; e.g. Henriques et al. 2015) or hydrodynamical simulations (Vogelsberger et al. 2020). Even though they rely on simple assumptions, they can accurately describe secondorder statistics of galaxies and matter. For example, halo model predictions for the galaxygalaxy twopoint correlation function agree well with observations (Zehavi et al. 2011). Mead et al. (2015) shows that halo models can also describe the nonlinear matter power spectrum in Nbody simulations with 5% accuracy.
Although halo models are prevalent for twopoint statistics, it is unclear whether they are sufficiently accurate for modelling higherorder statistics, such as the galaxy and matter bispectrum. These higherorder statistics contain complementary information to the twopoint statistics (Berlind & Weinberg 2002) so it is worthwhile to expand halo models to include them to improve and crossvalidate constraints from secondorder statistics. Here we extend halo models to measurements of G3L (Simon et al. 2008, 2013; Linke et al. 2020a), a thirdorder statistic.
The G3L signal is induced on a background galaxy (the ‘source’) by the weak gravitational lensing of matter around a pair of foreground galaxies (the ‘lenses’). It directly depends on the galaxymatter bispectrum integrated over the spread of lenses along the lineofsight (Schneider & Watts 2005). Unlike galaxygalaxy lensing, G3L is sensitive to the mean number of galaxy pairs inside halos and, therefore, the correlation of halo satellite numbers. If the galaxies in a lens pair belong to different samples, the G3L signal is higher for positively correlated satellite numbers and lower if they are anticorrelated. While the default assumption in the literature is that satellite numbers are uncorrelated (Scranton 2001, 2002; Zehavi et al. 2005), some galaxyclustering studies suggest a correlation between galaxy populations, such as red and blue galaxies (Zehavi et al. 2011; Ross & Brunner 2009; Wang et al. 2007; Simon et al. 2009).
We challenge the default assumption by constructing a halo model for the galaxygalaxymatter bispectrum and the G3L signal, inferring the HODs and satellite correlations for various galaxy samples. Our model is based on the approaches by Zheng et al. (2007), Zehavi et al. (2011) and Clampitt et al. (2017), but goes further by including the correlation between the satellite numbers of different galaxy samples. We validate the model and our inference procedure with G3L estimates in a simulated lensing survey based on the Millennium Simulation (MS; Springel et al. 2005) populated with SAM galaxies by Henriques et al. (2015). As a first realdata application, we analyze G3L measurements in the overlap region of KiDS, VIKING, and GAMA, with lenses selected from GAMA and sources from KiDS+VIKING (Linke et al. 2020a). We infer the HODs of galaxies in five different stellar mass bins between 10^{8.5} to 10^{11.5} h^{−2} M_{⊙} and two colours (red and blue), as well as the crosscorrelations between the perhalo satellite numbers for these galaxy samples.
This paper is structured as follows: In Sect. 2, we give an overview of the basics of G3L. We present our halo model for G3L in Sect. 3. The simulated and observed data sets are described in Sect. 4. Section 5 describes the estimator for the G3L signal and our analysis of the measurements. We present our results in Sect. 6 and discuss them in Sect. 7.
In the analysis of the MS, we use the cosmological parameters of this simulation, namely Ω_{m} = 0.25, Ω_{b} = 0.045, H_{0} = 73 km s^{−1} Mpc^{−1}, and σ_{8} = 0.9. For the analysis of the observations, we use the parameters from the Planck Collaboration I (2020), namely Ω_{m} = 0.315, Ω_{b} = 0.049, H_{0} = 67.4 km s^{−1} Mpc^{−1}, and σ_{8} = 0.811. Throughout we assume a flat cosmology with Ω_{Λ} = 1 − Ω_{m}.
2. Theory of galaxygalaxygalaxy lensing
G3L is a weak gravitational lensing effect (see, e.g. Bartelmann & Schneider 2001 for a review on weak lensing), first described by Schneider & Watts (2005, SW05 hereafter). There are two types of this effect: The lensing of background galaxy pairs by matter around individual foreground galaxies (lensshearshear G3L) and the lensing of individual background galaxies by matter around foreground galaxy pairs (lenslensshear G3L). We concentrate on the latter.
Figure 1 shows the geometric configuration for a lenslensshear G3L system on the sky. A background galaxy (‘source’), located at angular position θ, is gravitationally lensed by matter around two foreground galaxies (‘lenses’), located at θ + ϑ_{1} and θ + ϑ_{2}. Due to the lensing, the source experiences a tangential shear γ_{t}, which is measured with respect to the bisector of the angle ϕ between the sourcelens connections. The main observables for G3L are the threepoint correlation function and the closely related aperture statistics ⟨𝒩^{a}𝒩^{b}ℳ⟩, where a and b denote the considered galaxy populations. We describe the observables in Sects. 2.2 and 2.3. To model and ⟨𝒩^{a}𝒩^{b}ℳ⟩, we first introduce the galaxygalaxymatter bispectrum and its relation to the matter and galaxy number density distributions.
Fig. 1. Geometry of a G3L configuration with one source and two lens galaxies on the sky; adapted from Schneider & Watts (2005). Lens galaxies are at angular positions θ_{1} = θ + ϑ_{1} and θ_{2} = θ + ϑ_{2} on the sky; the source galaxy is at θ. The angle between the sourcelens connections is the opening angle ϕ. 
2.1. Projected fields and the galaxygalaxymatter bispectrum
The distribution of matter and galaxies is defined by the density ρ(x, z) and discrete galaxy number density at comoving position x and redshift z. The subscript a refers to the ‘sample’ of the galaxies, chosen with the same selection function, for example early or latetype galaxies. Fluctuations in the densities ρ and are the matter and galaxy number density contrast, δ(x, z) and , defined as
and
where bars denote ensemble averages.
Since G3L involves the correlation of galaxy pairs with the matter field, it is natural to derive its model from the galaxygalaxymatter bispectrum at comoving modes k_{1}, k_{2}, and k_{3} and redshift z. This bispectrum is defined by
where hats denote Fourier transforms for which we use the convention
The bispectrum is defined only for k_{3} = −k_{1} − k_{2}, so we abbreviate
As a gravitational lensing effect, G3L depends mainly on the projections κ and of δ and along the lineofsight. In a flat universe the lensing convergence is
where χ is the comoving distance, z(χ) is the redshift at χ, x(θ, χ) = (χ θ, χ), and q(χ) is an integral over the probability distribution function p_{s}(χ) of source galaxies with comoving distance χ, given as
The projected galaxy number density contrast is given by
where p^{a}(χ) is the probability distribution function of lens galaxies from sample a in comoving distance. This distribution strongly depends on the selection criteria for sample a. With p^{a}(χ) we also define the projected galaxy number density as
To arrive at a G3L signal, we convert the galaxygalaxymatter bispectrum to its projected counterpart , defined by
where hats again denote Fourier transforms. Similar as for , we abbreviate
As shown in Kaiser (1992), under the assumptions of the Limber equation, the projected bispectrum is
which we use to model the observables of G3L: the correlation function and the aperture statistics ⟨𝒩^{a}𝒩^{b}ℳ⟩. Equation (12) requires that the product q(χ) p^{a}(χ) p^{b}(χ) in the integrand does not vary strongly on scales smaller than the typical correlation length (i.e. a few megaparsecs) of the galaxy distribution (Bartelmann & Schneider 2001).
2.2. G3L correlation function
The G3L correlation function correlates the projected lens galaxy number densities and of two lens galaxy samples a and b with the tangential shear γ_{t}(θ) of the source galaxies,
Due to the statistical homogeneity and isotropy of the matter and galaxy fields, only depends on ϑ_{1} = ϑ_{1} and ϑ_{2} = ϑ_{2} of the lenssource separations, and the angle ϕ between ϑ_{1} and ϑ_{2}. For a ≠ b, the two lenses in each G3L configuration are from two different samples (mixed lens pairs); otherwise, they are from the same sample (unmixed lens pairs).
As shown in SW05, can be calculated from the projected galaxygalaxymatter bispectrum (cf. their Eq. (40)), by integrating (ℓ_{1}, ℓ_{2}), multiplied with a complex kernel function containing the secondorder Bessel function. Consequently, can be modelled from the bispectrum and directly compared to measurements. However, it is numerically preferable to consider the G3L aperture statistics, a linear transform of . We describe them in the following subsection.
2.3. G3L aperture statistics
Aperture statistics are moments of aperture masses ℳ_{θ}(ϑ) and aperture number counts ,
and
with a filter function U_{θ} of aperture scale radius θ. The function U_{θ} has to be a compensated filter, that is ∫ⅆϑ ϑ U_{θ}(ϑ) = 0. For G3L, the relevant aperture statistics are ⟨𝒩^{a}𝒩^{b}ℳ⟩, given by
The aperture statistics are related to the bispectrum by
where is the Fourier transform of U_{θ}. We choose the exponential filter function by Crittenden et al. (2002),
which is commonly used for studies of higherorder aperture statistics (e.g. Schneider et al. 2005; Jarvis et al. 2004; Simon et al. 2009; Saghiha et al. 2017) due to its favourable analytical properties. In particular, SW05 show that for this filter function the correlation function can be connected analytically to ⟨𝒩^{a}𝒩^{b}ℳ⟩ through
with the kernel function 𝒜_{𝒩𝒩ℳ}(ϑ_{1}, ϑ_{2}, ϕ ∣ θ_{1}, θ_{2}, θ_{3}) given in the appendix of SW05.
As alluded to before, there are practical advantages to modelling ⟨𝒩^{a}𝒩^{b}ℳ⟩ instead of . First, evaluating Eq. (17) is numerically more stable than calculating from the bispectrum because the filter function U_{θ} is localised and nonoscillating. Second, in contrast to , the aperture statistics do not depend on the galaxymatter twopoint correlation function. Therefore, they do not need a model of the galaxygalaxy lensing signal ⟨γ_{t}⟩. Third, a ⟨𝒩^{a}𝒩^{b}ℳ⟩ data vector is a condensed summary statistic, where a few aperture radii (∼tens) contain a similar amount of information as over hundreds of bins.
Consequently, we use the aperture statistics ⟨𝒩^{a}𝒩^{b}ℳ⟩ as the primary observable and model it with Eq. (17), based on a halomodel based bispectrum . To measure the aperture statistics, we estimate (see Sect. 5.1) and convert it to aperture statistics with Eq. (19). We focus on the equilateral statistics, that is, θ_{1} = θ_{2} = θ_{3}, to reduce the size of the data vector and use the shorthand
for convenience. We describe the ‘ingredients’ of our halo model and the derivation of the model bispectrum in the next section.
3. G3L halo model
In this section, we derive the galaxygalaxymatter bispectrum based on the halo model assumption, from which we can predict the G3L aperture statistics with Eq. (17). For this derivation, we require several ingredients: the linear matter power spectrum and critical density contrast, the halo density profile, the halo mass function (HMF), the halo bias, the spatial distribution of galaxies within a halo, and the first and secondorder moments of the HOD. We give details on these ingredients next.
3.1. Ingredients
Halo models assume that all matter and galaxies are within selfbound, virialized halos. Halos form in regions at redshift z where the linear density contrast of matter exceeds the critical density contrast δ_{c}(z). They reach the virial density , where is the cosmic mean density of matter and Δ(z) is the fractional overdensity within the virialized region.
The first ingredients are the linear matter power spectrum P_{lin} and the critical density contrast δ_{c}. We use the P_{lin} by Eisenstein & Hu (1998) that includes baryonic effects, and the fitting formula by Nakamura & Suto (1997) for the critical density contrast,
which was derived for flat ΛCDM universes. For Ω_{m} = 1, it reduces to the value for an Einstein–de Sitter Universe, δ_{c} ≃ 1.69.
Second, we require the density profile ρ(r − r_{0}  m, z) of a halo with mass m and redshift z, centred at r_{0}. We assume halos follow (truncated) NavarroFrenkWhite profiles (NFW; Navarro et al. 1996), given as
where r_{200} is the radius of a sphere around the halo centre, in which the mean density^{1} is 200 times the Universes mean density, m is the mass inside this sphere, c is the halo concentration parameter and u is the normalised density profile.
We model the mass and redshift dependence of c with the fitting formula by Bullock et al. (2001),
with c_{0} = 9 and α = 0.13. The mass m_{⋆} is that enclosed by a sphere of radius r_{⋆},
where r_{⋆} is the scale at which the standard deviation σ(r, z) of linear density fluctuations is equal to the critical overdensity δ_{c}, that is σ(r_{⋆}, z) = δ_{c}(z). The σ(r, z) is given by the convolution of the linear matter power spectrum P_{lin}(k, z) with the Fourier transformed tophat filter ,
The third group of ingredients are the HMF n(m, z)ⅆm and the halo bias. The HMF describes the comoving number density of dark matter halos with mass between m and m + ⅆm. We use the HMF by Sheth & Tormen (1999),
where the parameters A = 0.322, p = 0.3, and q = 0.707 were found in Nbody simulations, and ν = δ_{c}/σ(r_{m}, z) for .
The halo bias quantifies the clustering of halos by the ratio of the halo density contrast, δ_{h}, and the matter density contrast, δ. We assume a linear halo bias,
neglecting terms higher than linear in δ. Using the socalled peakbackground split formalism (Mo & White 1996; Scoccimarro et al. 2001) the linear bias b_{1} is
with the q and p as in the ShethTormen HMF. Assuming a deterministic halo bias, the halo power spectrum is then
and the halo bispectrum is
for the linear matter bispectrum by Bernardeau et al. (2002),
where k_{3} = −k_{1} − k_{2} and
Fourth, for the average number density of satellite galaxies a inside halos of mass m, , we assume an NFW profile with concentration
This concentration may differ from that of the halo matter if f^{a} ≠ 1. A similar parameter was, for example, introduced by Cacciato et al. (2012) in galaxyclustering halo models, who show that f^{a} affects galaxy clustering at scales below 1 h^{−1} Mpc.
Finally, we express the first– and second–order moments of the HODs by the model of Zheng et al. (2007), where galaxies are split into centrals and satellites, each with their own expected number per halo, such that the expected number ⟨N^{a}  m⟩ of galaxies from sample a in a halo of mass m is
Each halo hosts at most one central galaxy, , situated at the halo centre but can contain several satellite galaxies, . With this split, the galaxygalaxymatter bispectrum requires the specification of the mean numbers ⟨⟩ and ⟨  m⟩ of central and satellite galaxies, as well as the numbers of central and satellite pairs ⟨⟩, ⟨⟩, for a ≠ b, and ⟨( − 1)  m⟩ for a = b.
The mean number ⟨⟩ of central galaxies depends only on halo mass. For small halo masses, no galaxy formation occurs, so ⟨⟩ = 0. Halos with masses above a certain threshold, though, will contain central galaxies, at most one per halo. Similar to Zheng et al. (2007), we assume
with the free parameters α^{a}, , and σ^{a}. The mass M_{th} is the halo mass below which we do not expect halos to contain galaxies. The parameter σ^{a} determines the transition of ⟨⟩ from 0 to α^{a}. If σ^{a} is small, the transition from ⟨⟩ = 0 to ⟨⟩=α^{a} occurs quickly, whereas the transition is slower for larger σ^{a}. The parameter 0 ≤ α^{a} ≤ 1, for the maximum of ⟨N^{a}m⟩, gives the fraction of massive halos (m ≫ ) with a central galaxy from population a. Its inclusion is necessary because we are splitting galaxies into samples: only one sample can contain the central galaxy in a halo at a time. Since no more than one central galaxy per halo is allowed, the sum of all α^{a} from disjunct samples can never exceed unity.
For the mean number ⟨  m⟩ of satellites, we assume, based on Zehavi et al. (2005),
with the free parameters M^{′a} and β^{a}. The satellite number therefore follows the central galaxy number for small halo masses and becomes a power law at high halo masses. To illustrate the dependence of the HOD terms on halo mass, Fig. 2 shows the expected number of satellite and central galaxies and the parameters in Table 1. The total number of galaxies depends strongly on the central galaxy distribution for lowmass halos. Satellite galaxies predominate in massive halos.
Fig. 2. Mean perhalo numbers of galaxies for fiducial HOD parameters in Table 1. The solid black line shows the total galaxy number per halo, the dashed red line shows the fraction of halos with central galaxies, and the dotted blue line shows the number of satellite galaxies per halo. 
Fiducial values and flat priors of the halo model parameters.
The galaxygalaxy matter bispectrum depends not only on the expected numbers of centrals and satellites, but also on the expected number ⟨ ( − 1)  m⟩ of unmixed satellite pairs, the expected number ⟨  m⟩ of mixed satellite pairs, i.e, a ≠ b, and the number ⟨  m⟩ of centralsatellite pairs per halo of mass m. In principle, it also depends on the number ⟨⟩ of central pairs, but this number is 0 for all m, as each halo contains only one central galaxy. We assume, as, for example in Cacciato et al. (2012), that the occupation numbers of centrals and satellites are statistically independent,
This common assumption entails that the mean number of type b satellite galaxies in halos of mass m is independent of whether a central galaxy a is present or not. Furthermore, following Kravtsov et al. (2004), we assume that satellite occupation numbers vary according to a Poisson statistic,
We test this assumption and its impact on the model accuracy with galaxies inserted via SAM in a cosmological simulation in Appendix B.
Finally, we introduce as a new parameter the crosscorrelation coefficient r^{ab} of satellites, defined for mixed galaxy samples a ≠ b by
The coefficient r^{ab} is negative if the numbers of satellite galaxies a and b are anticorrelated; and positive if they are positively correlated. If the satellite galaxies are distributed Poissonian, that is, σ^{2}(  m) = ⟨  m⟩, then r^{ab} is a Pearson correlation coefficient. As we show in Sect. 6, mock galaxies from a SAM imply a mass dependence of r^{ab}(m), which we model by
with the free parameters A^{ab} and ϵ^{ab}. A value of A^{ab} = 0 corresponds to uncorrelated galaxy samples. The crosscorrelation between galaxy samples is independent of halo mass if ϵ = 0. Ignoring a mass dependence of r^{ab} in a more simplistic approach, as, for example in Simon et al. (2009), would result in a weighted average over the true r^{ab}(m). The comparison of such an average value to the true r^{ab}(m) in the simulation would not be straightforward, as the weighting per halo mass, given the survey characteristics, needs to be determined. A halomass dependent r^{ab}(m), on the other hand, is easier to interpret, which is why we use the massdependent fitting formula above.
A decrease of r^{ab}(m) towards small halo masses can be understood qualitatively as follows. The mean satellite number of galaxies from sample a inside a halo decreases with halo mass until a regime is reached where we find at most one. Likewise, at small enough halo masses, we find at most one galaxy from a second sample b (other galaxy populations may be present as well but are not of interest here). Then there are basically only four possibilities to populate a halo: (i) one galaxy from sample a, (ii) one galaxy from sample b, (iii) a pair of galaxies from a and b, and (iv) no galaxies from a or b. We shall denote the probabilities of these cases (i) to (iv) by p_{a}, p_{b}, p_{ab}, and 1 − p_{a} − p_{b} − p_{ab}, respectively. For this Poisson process with four outcomes, we find for the mean number of galaxies a in a halo
and ⟨⟩=p_{b} + p_{ab} for galaxies b. Similarly, the mean number of pairs is ⟨ ⟩=p_{ab}, and the variances are σ^{2}() = ⟨()^{2}⟩−⟨⟩^{2} = (p_{a} + p_{ab})(1 − p_{a} − p_{ab}) and σ^{2}() = (p_{b} + p_{ab})(1 − p_{b} − p_{ab}). The Pearson crosscorrelation coefficient of variations in the galaxy numbers consequently becomes
Therefore, at the extreme end, for halos too small to host more than one galaxy, that is, p_{ab} = 0, the value of exactly converges to zero for p_{a}, p_{b} → 0 (that is ⟨⟩, ⟨⟩→0). But, already in the intermediate regime, where 0 ≤ p_{ab} ≤ 1 − p_{a} − p_{b} and p_{a}, p_{b} ∼ 0.1, the upper limit of in the last line confines the correlation factor to smaller values: for example, r^{ab} ≤ 0.6 for p_{a} = p_{b} = 0.1, or r^{ab} ≤ 0.81 for p_{a} = 0 and p_{b} = 0.1. Therefore, a decreasing for m → 0 is a natural outcome of halos sparsely populated with galaxies a and b. This also applies to r^{ab} for subPoisson variances, expected at low halo masses, since r^{ab} < in this case.
As stated above, r^{ab} is a Pearson correlation coefficient only for Poisson variances of the satellite numbers, as assumed in our halo model. We show in Appendix B, however, that this assumption is inaccurate for the SAMs, and possibly real galaxies, at high or low halo masses, with specifics depending on the galaxy selection. In extreme cases, variances may be an order of magnitude larger for ⟨Nm⟩≫1 or smaller for ⟨Nm⟩≪1. Consequently, r^{ab} is at small halo masses systematically lower than the Pearson correlation coefficient (r^{ab} < ) and higher for high halo masses (r^{ab} > ). Nevertheless, true and reconstructed correlation parameters in this paper are comparable because we use the same r^{ab} definition in both cases (and in the science verification). In addition, our following verification results show that the strict assumption of Poisson satellites still allows an accurate HOD reconstruction within the statistical errors of KV450×GAMA data.
After describing our choices for the halo model parameters, we can model the galaxygalaxymatter bispectrum B_{ggδ}. We derive the bispectrum in the following subsection.
3.2. Galaxygalaxymatter bispectrum
Using Eqs. (1)–(3), the galaxygalaxymatter bispectrum is given by
where a and b denote the galaxy samples, and ‘unconnected terms’ are those proportional to δ_{D}(k_{1}), δ_{D}(k_{2}), or δ_{D}(k_{3}), which do not affect the bispectrum. The bispectrum can be divided into three terms: the 1halo term _{1}(k_{1}, k_{2}, z), the 2halo term _{2}(k_{1}, k_{2}, z), and the 3halo term _{3}(k_{1}, k_{2}, z), or together
The 1halo term depends on the correlation between galaxies numbers and matter density in the same halo. One part of the 2halo term is caused by the correlation of galaxies in one halo with the matter in a different halo. Its second part is due to the correlation of galaxies and matter in one halo with galaxies in another halo. Correlations between matter and galaxies in three distinct halos cause the 3halo term.
To derive and , we use that in the halo model the cosmic density field consists of H halos with masses {m_{1}, …, m_{H}} at positions {x_{1}, …, x_{H}}. With the normalised density profile u the matter density field is
whose Fourier transform is
Galaxies are treated as discrete point particles. Each satellite galaxy j from sample a belongs to a halo centred at x_{i} and is at separation from the halo centre. Centrals are exactly at the halo centre. The number density of galaxies from sample a is therefore
where is the number of satellite galaxies and the number of central galaxies from sample a inside halo i. The Fourier transform of the number density is
To derive the bispectrum we insert Eqs. (49) and (51) into Eq. (46), which leads to
and
with
and
where the Kronecker symbol is 1 for a = b, and 0 otherwise. We derive these equations in full in Appendix A.
With the ingredients of the halo model from Sect. 3.1, the bispectrum is fully specified. By inserting it into Eq. (12) and (17) we obtain the aperture statistics. Before fitting the model to simulated and actual G3L measurements, we discuss the parameter sensitivity of the model in the next section.
3.3. Discussion of model parameters
We qualitatively study the impact of the parameters on the theoretical G3L signal. For this, we vary each parameter inside the prior range, adopted for the likelihood analysis below (Table 1), while keeping the other parameters fixed to ‘fiducial’ values, that is, to the centres of the parameter prior range. There are six free parameters for each galaxy sample and two additional parameters for the correlation of satellite numbers, hence in total 14 parameters for the bispectrum of a combination of galaxy samples. We arrange these parameters inside the vector p,
and write the modelled aperture statistics for a set p of parameters and scale radius θ as ⟨𝒩𝒩ℳ⟩(θ, p).
The priors are aimed to be uninformative. For σ, log_{10}(/M_{⊙}), β^{a}, and log_{10}(M^{′a}/M_{⊙}), they are based on the galaxygalaxy lensing halo model by Clampitt et al. (2017), whose prior ranges we doubled for the G3L analysis. For α, the prior range contains all possible values between 0 and 1. For the parameters f, A, and ϵ, which are unique to our G3L halo model, there are no reference values. Therefore, we chose arbitrary broad ranges centred on the values f = 0 (indicating galaxy distributions perfectly following dark matter halos), A = 1, and ϵ = 1 (indicating no correlations between galaxy HODs). To sanitycheck the dependence on priors, we performed the analysis on the validation data also with prior ranges two times larger and found no difference in the final parameter constraints.
Figure 3 shows ⟨𝒩^{a}𝒩^{a}ℳ⟩ for unmixed lens pairs, when varying α^{a}, σ^{a}, , M^{′a}, and f^{a}. All parameters, except for σ, visibly impact the aperture statistics. However, we note that the trends seen in Fig. 3 are for varying each parameter individually and do not show the correlation between the parameters. For example, the threshold mass and the satellite mass scale M^{′a} are tightly correlated: Increasing while decreasing M^{′a} could lead to the same ⟨𝒩^{a}𝒩^{a}ℳ⟩.
Fig. 3. Impact of halo model parameters on ⟨𝒩^{a}𝒩^{a}ℳ⟩ for unmixed lens pairs. In each panel, only one parameter is varied at a time. Solid lines indicate the total ⟨𝒩^{a}𝒩^{a}ℳ⟩, while dashed lines show the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term. 
The threshold hal o mass and the slope β^{a} of the satellite HOD affect the signal the strongest on the scales we consider here. Increasing the threshold mass from 3.2 × 10^{12} M_{⊙} to 7.5 × 10^{12} M_{⊙} roughly doubles the aperture statistics across the whole range of θ from to 100′ because is closely connected to the galaxy bias. A higher causes galaxies to reside in more massive halos, so the galaxy bias and the shear amplitude increase, which increases the amplitude of the G3L aperture statistics. The satellite parameter M^{′a}, which describes the mass scale where ⟨N_{sat}m⟩≈1 if ≪ M^{′a}, changes the signal amplitude for θ ≳ 1′. A smaller M^{′a} produces a higher amplitude by scaling up the number of satellites in a halo of a given mass.
The slope β^{a} of ⟨⟩ has the strongest impact on the statistics at θ ∼ 9′, an up to a 20fold increase in amplitude. While a larger β^{a} increases the signal on larger scales, at scales below , it decreases the aperture statistics by up to 50%. These opposite trends have two reasons. First, a larger β^{a} causes more massive halos to contain more satellite galaxies and pairs, giving the halos more weight in the second term in G^{a} (Eq. (A.27)) and the last term in G^{ab} (Eq. (A.18)), compared to lowmass halos with fewer satellite pairs. This increases ⟨𝒩^{a}𝒩^{a}ℳ⟩. Second, more massive halos are more extended and contribute predominantly at larger scales. Therefore, the signal increases at larger angular scales but decreases at scales where lowermass halos are more relevant.
The concentration parameter f^{a} affects the aperture statistics stronger at small scales than at large scales (57% at and 40% at 1′). A more concentrated galaxy profile (f^{a} > 1) leads to more galaxies in the inner regions of each halo, which increases the galaxy density contrast δ_{g} at small scales. Consequently, the galaxygalaxymatter bispectrum and the aperture statistics increase at small scales. At scales above , though, the change is less than 10%.
The sensitivity of the statistics for r^{ab} is visible in Fig. 4, where we show ⟨𝒩^{a}𝒩^{b}ℳ⟩ for mixed lens pairs for correlated satellite numbers (r^{ab} = 1), uncorrelated numbers (r^{ab} = 0), or anticorrelated numbers (r^{ab} = −1). For this example, we keep r^{ab} constant with halo mass, so ϵ^{ab} = 0, while A^{ab} is varied between A^{ab} ∈ {0, ±1}. The amplitude of ⟨𝒩^{a}𝒩^{b}ℳ⟩ changes strongest for (72%), where the signal is dominated by the 1halo term containing ⟨m⟩ and has a r^{ab}dependence. The crosscorrelation parameter has only a small effect (less than 10% change) at scales above 30′ because these scales are dominated by the 3halo term, independent of the perhalo number of satellite pairs.
Fig. 4. Impact of the correlation of satellite numbers on the aperture statistics of mixed lens pairs. Satellite numbers are either fully correlated (r^{ab} = 1, blue lines), uncorrelated (r^{ab} = 0, black lines) or anticorrelated (r^{ab} = −1, red lines). Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term. 
In conclusion, the aperture statistics between and 100′are most sensitive to the threshold halo mass and the slope β^{a}. In contrast, we do not expect G3L to improve the σ^{a} constraints beyond the prior range. To probe the crosscorrelation parameter r^{ab}, the statistic ⟨𝒩^{a}𝒩^{b}ℳ⟩ needs to be measured at scales below 30′and preferably around 1′. Both of these measurements are achieved for our survey data.
4. Data
Before describing the G3L measurement procedure and halo model fitting procedures, we briefly explain the simulated and observed data. These are the same data sets as in Linke et al. (2020a): simulated (verification) data based on the MS and the SAM by Henriques et al. (2015), and one composed of the overlap of KiDS, VIKING, and GAMA. As mentioned in Sect. 1, the cosmology used to analyze the observation differs from the fiducial cosmology of the MS and is based on Planck Collaboration I (2020). This difference in cosmology should not impact our conclusions because we are not interested in a onetoone comparison between the observation and simulation and merely use the mock data to validate our method.
4.1. KV450 × GAMA
We use observational data from the overlap of KiDS, VIKING, and GAMA, KV450 × GAMA for short, which is an area of approximately 180 deg^{2}. KiDS (Kuijken et al. 2015; de Jong et al. 2015) and VIKING (Edge et al. 2013; Venemans et al. 2015) are two photometric surveys covering approximately the same 1350 deg^{2} area, with KiDS observed in the optical with the VLT survey telescope, and VIKING observed in the nearinfrared with the VISTA telescope. We use the public combined data release KV450, described in detail in Wright et al. (2019). Galaxies were observed in the u, g, r, i band for KiDS and Z, Y, J, H, K_{s} bands for VIKING, with shape measurements performed in the rband. This data set was processed with AstroWISE (de Jong et al. 2015) and THELI (Erben et al. 2005; Schirmer 2013) to give a catalogue of observed galaxies. The shapes of these galaxies were measured with lensfit (Miller et al. 2013; Kannawadi et al. 2019). We use the KV450 galaxies with photometric redshifts between 0.5 and 1.2 (obtained by Hildebrandt et al. 2020) as sources for the G3L measurements. The cosmological analysis of KV450 found that the n(z) for individual tomographic bins of the KV450 galaxies are offset by values between −0.007 and 0.013 (Hildebrandt et al. 2020).
Our lens galaxies are from GAMA (Driver et al. 2009; Liske et al. 2015), a spectroscopic survey conducted at the Anglo Australian Telescope. We use all galaxies listed in the table distanceFramesv14 from the data management unit (DMU) LocalFlowCorrection with a redshift quality parameter N_Q > 2 and with a spectroscopic redshift less than 0.5 to avoid overlap between lenses and sources. Each lens galaxy is assigned absolute magnitudes, restframe photometry, and stellar masses according to the table stellarMassesLambdarv20 from the DMU StellarMasses (Taylor et al. 2011). These were obtained with matched aperture photometry and the LAMBDAR code (Wright et al. 2016), assuming the initial mass function by Chabrier (2003), stellar population synthesis according to Bruzual & Charlot (2003), and dust extinction according to Calzetti et al. (2000). To calculate the angular correlation function ω^{ab}, we use the randoms in the table randomsv02 from the DMU Randoms, created by Farrow et al. (2015).
Galaxies observed by GAMA are brighter than r = 19.8 mag, rendering our lens galaxy sample fluxlimited. We divide the GAMA galaxies into different samples: a ‘red’ and a ‘blue’ sample, using their restframe (g − r)_{0} colour and their absolute magnitude M_{r} in the rband, and five samples defined by stellar masses. For this we use the same cuts as in Farrow et al. (2015) and Linke et al. (2020a). An overview of these cuts is given in Table 2.
Selection criteria for lens samples.
These different samples have differing redshift distributions n(z), shown in Fig. 5. While the red and blue galaxies are distributed similarly, the n(z) of the stellarmass selected lenses differ strongly. Galaxies with lower stellar mass are predominantly observed at smaller redshifts; galaxies with higher stellar mass are found up to the limiting redshift 0.5. These differences are caused by the flux limit of the survey; galaxies with lower stellar masses, typically fainter than more massive galaxies, are only visible at smaller redshifts.
Fig. 5. Normalised redshift distributions n(z) of GAMA galaxies, selected by colour (left) and stellar mass (right). 
The differences in the redshift distribution of galaxies from different samples are also visible in Fig. 6, which shows the stellar mass M_{*} and redshift z of the GAMA galaxies. Observed galaxies at higher redshift have higher stellar masses on average, which is a direct consequence of the flux limit. Red and blue galaxies are distributed similarly with redshift but differ in stellar mass: red galaxies tend to have higher stellar masses than blue galaxies.
Fig. 6. Stellarmass and redshift of GAMA galaxies, divided by colour (left) and stellar mass (right). 
4.2. Millennium Simulation
The simulation data are constructed from the MS (Springel et al. 2005), a darkmatter only cosmological Nbody simulation. Its simulation box has a comoving side length of 500 h^{−1}Mpc and contains 2160^{3} particles with mass 8.6 × 10^{8} h^{−1} M_{⊙} (h = 0.73).
Maps of the lensing shear, γ, were created with the multiplelensplane raytracing algorithm by Hilbert et al. (2009). This algorithm generated shear maps of size 4 × 4 deg^{2} on a regular grid of 4096 × 4096 deg^{2} for 64 linesofsights (total area of 1024 deg^{2}) at each redshift slice of the MS^{2}. We combine the shear redshift slices in a weighted average according to the observed KV450 galaxy redshift distribution, resulting in a combined shear without shear noise and intrinsic source alignment. The effective source density for the science verification corresponds to the pixel density, that is, 291 arcmin^{−2}.
Since all 64 linesofsight originate from the same simulation, they are not independent at the largest scales. However, the linesofsight correspond to different observers placed in the simulation box such that the overlaps between the linesofsights are minimal, so we expect correlations predominantly at large scales, close to the simulation box size of 500 h^{−1} Mpc. In contrast, the G3L signal is dominated by the 1halo term, originating from small scales of approximately 1 h^{−1}Mpc. Therefore, we treat the linesofsight as independent for our purposes.
The lens galaxies for the simulated lenses Henriques et al. (2015) use the same initial mass function (Chabrier 2003) as assumed for the stellar mass estimates from GAMA, but a different stellar population model (Maraston 2005). Nevertheless, as shown in Linke et al. (2020a), the predictions of the MS combined with this SAM for the G3L signal agree with observations in KV450×GAMA.
To mimic the GAMA selection function, we use all lens galaxies with redshift less than 0.5 and brighter than r = 19.8 mag. We also divide the simulated lens galaxies by their colour and stellar masses by the same cuts as for the observed lens galaxies (Table 2).
5. Methods
This section outlines our G3L estimators and the model fit to the data. The measurement procedure largely follows Simon et al. (2008) and Linke et al. (2020b) and is summarised in Sect. 5.1. The statistical analysis for the inference of HOD parameters from the G3L data is detailed in Sect. 5.2. To assess the accuracy of the inferred HOD in a science verification with mock data, later on, we describe in Sect. 5.3 the determination of the true HODs of the simulated galaxies. We make our codes publicly available at github.com^{3}.
5.1. Measuring G3L
We estimate the G3L aperture statistics from a catalogue of galaxy positions and source shapes for bins B of similar triangles (ϑ_{1}, ϑ_{2}, ϕ) with the estimator in Simon et al. (2008). For N^{a} lens galaxies from sample a, N^{b} lens galaxies from sample b, and N_{s} source galaxies with complex ellipticities ϵ_{k}, the estimate of is the real part of
where ω^{ab}(θ) is the angular twopoint correlation function of lens galaxies from samples a and b with lag θ, and
The angles φ_{ik} and φ_{jk} are the polar angles of the lenssource separation vectors θ_{i} − θ_{k} and θ_{j} − θ_{k} (corresponding to φ_{1} and φ_{2} in Fig. 1) and ϕ_{ijk} = φ_{ik} − φ_{jk} is the opening angle between the lenssource separation vectors (corresponding to ϕ in Fig. 1). The weight w_{k} of the source k (set to w_{k} ≡ 1 for our simulated data) measures the confidence in the shape measurements for the source, with higher weights indicating more precise ellipticities. The imaginary part of is pure noise in the absence of systematic errors in the shear data (Linke et al. 2020a).
We give in the estimator equal weight to all lens pairs, which is in contrast to the recent analysis in Linke et al. (2020a), where we weighted galaxy pairs based on the redshift difference between their members. The lack of weighting dilutes the signaltonoise ratio of the G3L signal since nonphysical pairs (widely separated along the lineofsight) carry no signal but increase the noise. However, a model with weights requires abandoning the Limber approximation for the projection of the bispectrum in Eq. (12). The weighting introduces an additional factor into the integral, which depends on the lineofsight distance between the lens galaxies in a pair. This factor, by definition, varies strongly on scales corresponding to the correlation length between galaxies. It is designed to give high weights to galaxy pairs within a correlation length while downweighting galaxies outside this range. Consequently, the assumption that the integrand in Eq. (12) is only slowly varying is no longer valid. With an alternative formalism for Eq. (12) unclear at this point, we apply equal weights to our lens pairs for the scope of this work.
Our estimator for the twopoint correlation function of the lens clustering, ω^{ab}(θ), is that by Szapudi & Szalay (1998),
for two lens samples a and b with and galaxies, and two ‘random samples’. These random samples contain and unclustered galaxies subject to the same selection functions as the lens samples. The D_{a}D_{b}, D_{a}R_{b}, D_{b}R_{a}, and R_{a}R_{b} are the pair counts of observed (D) and random galaxies (R).
Having obtained ω^{ab}, we measure (on an approximate flat sky) with (58) individually for tiles of size 1° ×1°. For this, we divide the observational data into N = 189 and the simulation data into N = 4 × 64 = 256 tiles. The tiles allow us to project the galaxy positions and shear to Cartesian coordinates via an orthographic transformation and to quantify the uncertainty of the with jackknife resampling on a tile by tile basis. We estimate for each tile i on a regular grid of 128 × 128 × 128 bins. These bins are spaced logarithmically for ϑ_{1} and ϑ_{2}, and linearly for ϕ. We use and ϕ ∈ [0, 2π]. For the total , individual tile estimates are averaged, weighted by the effective number of triplets per bin.
Due to the finite number of galaxies, some of the bins will remain ‘empty’, meaning that the data contains no lenslenssource triplet fitting the configuration of the bin. Setting the correlation function in these bins to an arbitrary value biases the measurement, so we use the adaptive binning from Linke et al. (2020b). This scheme uses Voronoi tesselation to redefine the bins, now b_{i}, such that no empty bins occur, effectively merging empty and ‘filled’ bins. We found in Linke et al. (2020b) that estimating in 128 × 128 × 128 bins and then applying the tesselation leads to a measurement accuracy within 5%.
We convert a binned to aperture statistics with a numerical approximation to Eq. (19),
where N_{bin} is the number of bins, V(b_{i}) is the tesselation volume of bin b_{i}, and 𝒜_{𝒩𝒩ℳ}(b_{i}  θ_{1}, θ_{2}, θ_{3}) is the kernel function of Eq. (19), evaluated at the centre of bin b_{i}.
We estimate the covariance matrix of the estimates with jackknife resampling. For this, we combine the estimates for all but the kth tile to the kth jackknife sample , which is then converted to the aperture statistics, Eq. (19), leading to N samples ⟨𝒩^{a}𝒩^{b}ℳ⟩_{k}. The estimate of the ⟨𝒩^{a}𝒩^{b}ℳ⟩ covariance matrix is then
where is the average of all aperture statistics jackknife samples. The (meansquare) statistical uncertainty of the aperture statistics ⟨𝒩^{a}𝒩^{b}ℳ⟩(θ_{i}) is . The inverse of this covariance, needed for the likelihood analysis, is estimated with
where N_{θ} is the number of aperture radii bins (Hartlap et al. 2007; Anderson 2003). In our case N_{θ} = 3 × 30 = 90.
Formally, jackknife resampling assumes that all individual estimates of are statistically independent. However, as all observed tiles originate from the same observation and are adjacent to each other, this assumption is not strictly valid. It is still probably a good approximation on scales smaller than the tile sizes. A possible bias in the empirical covariance is less pronounced for the simulation, as the tiles originate from 64 (mostly) independent lineofsights.
5.2. Fitting the halo model
We constrain the parameters of the G3L halo model by fitting it to measurements of both the autocorrelation aperture statistics, ⟨𝒩^{a}𝒩^{a}ℳ⟩(θ) and ⟨𝒩^{b}𝒩^{b}ℳ⟩(θ), and crosscorrelation statistics, ⟨𝒩^{a}𝒩^{b}ℳ⟩(θ), of two galaxy samples a and b for N_{θ} = 30 scale radii θ between and 30′. To this end, we combine the measurements into a data vector of 3N_{θ} elements,
Likewise, we define the halo model vector m^{ab}(p) for each parameter set p, which is obtained from Eq. (46) and the bispectrum in Sect. 3.2,
Optimal parameters p_{opt} are those that minimise the goodnessoffit
determined by the NelderMead algorithm (Nelder & Mead 1965) as implemented in the GNU Scientific Library (Gough 2009). This algorithm is wellsuited to multidimensional minimisation problems and requires only few (typically 1 or 2) function evaluations per iteration step. To avoid local minima, the algorithm is restarted multiple times at different, randomly chosen, initial parameter values.
To estimate the uncertainties on the bestfitting parameters p_{opt} and to quickly sample the posterior distribution of the parameters P(p  d^{ab}), we approximate the posterior distribution with the importance function q(p  d^{ab}) following the importance sampling scheme (e.g. Liu 2004). According to the Bayes theorem, the posterior density of p given the data is
where P_{prior}(p) is the prior density of the parameters, as given in Table 1, and ℒ(d^{ab}p) is the likelihood of the data given the parameters p. We assume Gaussian statistical errors in the data, that is, the likelihood function is
with the χ^{2} as defined in Eq. (66); the Bayesian evidence and thus the normalisation of ℒ is not of interest here and will be ignored in what follows. The importance sampling function q should be close to the posterior P to achieve an efficient sampling. To find an appropriate q, we approximate ℒ(d^{ab}  p) in the proximity of the optimal parameters p_{opt} by the Gaussian probability density
where the matrix F is the Fisher information,
(e.g. Tegmark et al. 1997), and choose q as
We now draw N_{p} parameter sets p_{i} from q which are then, by importance sampling, weighted to sample P. The allocated weights are
normalised such that their sum is unity.
For each parameter p^{i} we give the α credibility interval (CI) on the marginalised posterior P(p^{i}d^{ab}), defined by
We take the mode of the posterior as optimal parameter value and define the CI I_{α} of a single parameter as the interval around the optimal parameter value including α of the marginalised posterior, that is,
To find this interval, we sort the sampling points by the distance p^{i} −  from the optimal parameter value and take the first N_{p} points, for which the sum of their weights equals α. The interval defined by these points is our estimate for the α credibility region of the optimal parameter.
5.3. Science verification
We assess the accuracy of the model, estimators, and our statistical setup by comparing the inferred HODs to the true ones in mock data. In these mock data, each galaxy is associated with a dark matter halo, identified by its halo ID and virial halo mass from a friendsoffriends halo finder. To extract the true HOD, we count the number of galaxies of each sample a in each halo and divide the halos into 50 mass bins in the range between 10^{11} and 10^{15} h^{−2} M_{⊙}. For each mass bin, the number of satellite and central galaxies per halo is averaged, yielding the true average ⟨N^{a}  m⟩=⟨⟩+⟨  m⟩ in the data. We also verify the inference estimate of the correlation coefficient r^{ab}(m), Eq. (40), by computing
from the halo catalogue and galaxies in the mock data. The uncertainties on ⟨N^{a}  m⟩ and r^{ab}(m) are the standard errors on the mean over the 64 linesofsights of the simulation.
6. Results
In this section, we give the results of the science verification and the G3L analysis for KV450 × GAMA. We first present the results for lens samples defined by their colour in Sect. 6.1 and then for lens samples defined by their stellar mass in Sect. 6.2.
6.1. Colourselected lens samples
The G3L aperture statistics of simulated red and blue galaxies are shown in Fig. 7a, along with the best fit of the halo model and its decomposition into the 1, 2, and 3halo terms. The parameter values corresponding to the best fit are given in the first two columns in Table 3. The goodnessoffit is χ^{2} = 93.63 for 90 − 14 = 76 degrees of freedom (d.o.f.), or χ^{2}/d.o.f. = 1.28. This χ^{2} corresponds to a pvalue of 0.083, indicating no significant deviation between the fit and the simulated G3L signal within the 95% confidence level (CL).
Fig. 7. G3L measurement (points) and bestfitting halo model (lines) for red and blue galaxies in the MS (upper plot) and the KV450 × GAMA (lower plot). Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. The left panels show the result for redred galaxy pairs, the central panels for blueblue galaxy pairs, and the right panels for redblue mixed pairs. Error bars correspond to the standard deviation estimated from the jackknife resampling detailed in Sect. 5.1. (a) MS. (b) KV450 × GAMA. 
Bestfitting values of halo model parameters for colourselected lenses and 68% credibility intervals (d.o.f. = 76).
For our science verification, we compare in Fig. 8a the HODs inferred by the bestfitting G3L halo model to the directly estimated HODs of the simulated galaxies. The model prediction and direct estimate agree within the 68% credibility band (shaded areas) for red and blue galaxies. Likewise, the correlation of numbers of red and blue satellites, A^{redblue} and ϵ^{redblue}, is detected at 3σ significance in the verification, and r^{ab}(m) agrees in Fig. 9 with the true r^{ab} in the galaxy SAM within the 68% credibility band. Therefore, the G3L fit recovers the galaxy HODs and r^{ab}(m) sufficiently accurate within the statistical errors for a survey similar to our verification data (∼10^{3} deg^{2}, a high source number density without shape noise) and surveys with higher estimator noise, such as KV450 × GAMA. As an additional test, we explore the impact of a systematically wrong source n(z) in the HOD inference by shifting the full n(z)↦n(z − δz) by δz = 0.02, roughly twice the expected bias on the mean redshift of n(z) (see Sect. 4.1). The resulting systematic shift of the best fit HOD parameters stays within the 68% CI of the verification data, for example, δM_{th} = −1.8 × 10^{11} M_{⊙} (−0.05 × 10^{11} M_{⊙}) for the red (blue) galaxies. Systematic errors in the HOD due to bias in n(z) are hence negligible for our analysis of KV450×GAMA.
Fig. 8. Mean perhalo numbers of simulated galaxies (top) and observed galaxies (bottom) as function of halo mass. Red crosses (blue points) indicate the true HOD of simulated red (blue) galaxies, where the error bars are the standard deviation of the mean over the 64 lineofsights. The lines indicate the perhalo numbers inferred from the fit to the G3L signal for red (solid red) and blue galaxies (dashed blue). The shaded areas are the 68% credibility areas of the halo model fit. (a) MS. (b) KV450 × GAMA. 
We determine the HOD parameters of real red and blue galaxies in KV450×GAMA in Fig. 7b, where we show the aperture statistics and the best fits of our halo model, together with a decomposition into the 1, 2, and 3halo terms. The model fit has χ^{2}/d.o.f. = 0.977, corresponding to a pvalue of 0.53 and an agreement with the model within the 95% CL. The best fitting parameters are reported in the second pair of columns in Table 3. They show that red galaxies clearly populate more massive halos than blue galaxies: is roughly ten times larger than (68% credibility intervals, CI herafter). The perhalo number of blue satellites increases slower with halo mass than for red galaxies: the mass scale of blue satellites is more than five times larger than (68% CI). For central galaxies, the sum of α^{red} and α^{blue} is (68% CI), which is consistent with unity (68% CI). The concentration of halo satellites is consistent with that of matter (f^{a} ∼ 1) in the 68% CI. As expected from the qualitative analysis (Sect. 3.3) σ cannot be constrained better than the prior.
The HODs corresponding to the best fitting parameter values are shown in Fig. 8b. The HOD of red (blue) galaxies is nonzero for halo masses above 10^{12} M_{⊙} (5 × 10^{11} M_{⊙}), whereas, at lower halo masses, the constraints become essentially upper limits for ⟨Nm⟩. The inferred HODs for GAMA galaxies match those obtained from the fit to the simulated galaxies in Fig. 8b. This agreement reflects the similar G3L aperture statistics of mock data and observations – the SAM predictions for ⟨𝒩^{a}𝒩^{b}ℳ⟩ agree with the measurements in KV450 × GAMA within the errors (Fig. 7). However, the uncertainties on the HODs and their parameters are considerably larger for the observation since our simulated data has less noise in the shear signal.
As for r^{ab} of red and blue satellites in KV450×GAMA, we report a 2σ to 3σ detection of a positive correlation and an amplitude increase towards more massive halos: both A^{redblue} and ϵ^{redblue} are positive. This trend is similar to that of the simulated galaxies. However, the increase with halo mass is steeper for the observed galaxies (Fig. 9 versus at 68% CI), while the correlation amplitude at 10^{12} M_{⊙} is lower ( versus at 68% CI). Consequently, numbers of red and blue satellites are correlated both in the SAMs and for true galaxies, especially beyond the mass scale of galaxy groups m ≳ 10^{13} M_{⊙} where for ϵ^{red − blue} = 1 (68% CI).
Fig. 9. Correlation parameter r^{ab} for red and blue galaxies in the simulation and observation as a function of halo mass. Black crosses show the direct estimate for the simulation, where the error bars are the standard deviation of the mean over the 64 lineofsights. The solid brown line shows the r^{ab} inferred by the halo model fit to the G3L signal of the MS, and the green dashed line is the result of the fit to the KV450 × GAMA G3L signal. The shaded areas show the 68% credibility bands of the fits. 
The correlation matrix of the G3L estimate for the KV450 × GAMA red and blue galaxies is shown in Fig. 10. We see that the signal for similar aperture radii is strongly correlated. The signals for redred and blueblue lens pairs are almost independent, while the signal for mixed galaxy pairs has correlation coefficients of up to 0.3 with the signal for unmixed redred pairs.
Fig. 10. Correlation matrix for aperture statistics measurement in KV450×GAMA for red and blue galaxy samples. The data vector is ordered as given in Eq. (64), starting with the smallest aperture radius. 
6.2. Stellar massselected lens samples
We repeat the science verification and G3L analysis for the stellar massselected galaxies. We consider five stellarmass bins m1 to m5, so there are ten distinct combinations of two samples a and b. For each combination, we again use the statistics ⟨𝒩^{a}𝒩^{a}ℳ⟩, ⟨𝒩^{a}𝒩^{b}ℳ⟩, and ⟨𝒩^{b}𝒩^{b}ℳ⟩ to infer HOD parameters. Since each galaxy sample is used in four combinations, we obtain four versions of the same HOD for each sample. For a realistic model and successful fits, these four HODs should be consistent with each other.
Fitting the halo model individually to two samples out of five neglects the correlations between the ten combinations. Consequently, better constraints on the model parameters could be obtained by fitting the model to the signal of all five stellarmass selected samples simultaneously. However, this would increase the data vector from 90 entries to 450 entries. As the covariance estimated is obtained from only 180 (quasi)independent data realisations, it cannot be used with such a large data vector, and a simultaneous fit of all five samples is unfeasible.
The HOD parameters, χ^{2}, and plots of the best fitting models for the science verification with the simulated data are given in Appendix C. The pvalues exceed 0.05 for all fits, indicating no significant deviation between the fits and the measurements within the 95% CL. To evaluate the overall agreement of all fits to the mock data, we consider the cumulative distribution of pvalues, shown in Fig. 11 (solid line). This distribution should correspond to a uniform distribution between 0 and 1 (dotted line) if the measurements are unbiased realisations of the model. A KolmogorovSmirnov (KS) test on the distribution of pvalues for the simulation yields a KS distance of 0.118, which for the 11 samples and d.o.f. = 76 in the distribution indicates no significant deviation from a uniform distribution at 95% CL. Additionally, the four parameter sets for each stellarmass bin agree within the 68% CI. The distribution of pvalues and the consistency of the HOD parameters for the four model fits supports the view that the model coherently and accurately reproduces the G3L signal of the verification data within the statistical errors.
Fig. 11. Cumulative distribution of pvalues of G3L halo model fits for MS (orange, solid) and KV450 × GAMA (green, dashed). For a perfect description of G3L signal and data noise, the distributions would be consistent with a uniform distribution (black, dotted). 
To verify the reconstruction of the HODs for stellarmass samples, we compare the inferred HODs (lines and shaded areas) to the true HODs (data points) in Fig. 12a, which only shows the reconstructions for the combinations m1–m5, m2–m5, m3–m5, and m4–m5; the other reconstructions for the same sample but in a different combination agree with those within the 68% CI. The inferred HODs agree with the true HODs within the 68% credibility band. However, for the lower stellar mass galaxies from samples m1, m2, and m3, there is a ‘dip’ for the true ⟨Nm⟩ with a local minimum near 6 × 10^{11} M_{⊙}, 1 × 10^{12} M_{⊙}, and 2 × 10^{12} M_{⊙}, respectively. Our halo model cannot trace this feature because, by construction, ⟨Nm⟩ increases monotonically with halo mass m; the halo model fits a smooth profile across the dip. However, this smoothing only increases the HOD reconstruction uncertainty without introducing a significant bias within the 68% CI. Therefore, at the level of the precision of the verification data and for the noisier observational data, the missing HOD model feature is acceptable for KV450×GAMA.
Fig. 12. Mean perhalo galaxy numbers in the simulation (top) and observation (bottom) for lens galaxies from each stellar mass bin as function of halo mass. Crosses indicate the directly estimated perhalo numbers of simulated galaxies, while lines show the predictions from the G3L fits. The shaded areas indicate the 68% confidence areas. Left panels: the mean perhalo numbers for galaxies from stellar mass bins m1, m3, and m5 obtained from the fits to the G3L signal for m1–m5, m3–m5, and m4–m5. The right panels show the same for galaxies from stellar mass bins m2 and m5, obtained from the fits to the signal for m2–m5 and m4–m5. The corresponding HOD parameters are listed in Table C.1. (a) MS. (b) KV450 × GAMA. 
To verify the inference of r^{ab}, we compare the model fits (solid lines and 68% credibility regions) to the true correlation of satellite numbers (data points) in the simulation in Fig. 13. The corresponding values of A^{ab} and ϵ^{ab} are listed in Table C.1 in the Appendix. For all combinations of stellarmass bins, r^{ab} is positive and scales approximately linearly with halo mass (ϵ^{ab} ≈ 1). The amplitude A^{ab}, though, drops if one of the samples is m4 or m5 in the combination. For example, A^{m1m2} is , whereas A^{m1m5} is (68% CI). The r^{ab} from the halo model fits agree with the true r^{ab} for all sample combinations at the 68% CI (cf. data points to lines and dark green areas).
Fig. 13. Correlation parameter r^{ab} for stellar massselected galaxies in the MS and in KV450×GAMA. Black crosses show the true correlation for the MS, where the error bars are the standard deviation of the mean over the 64 simulated lineofsights. The solid brown line shows r^{ab} inferred from a simulated G3L analysis of MS, and the green dashed line the inference for KV450×GAMA. The shaded areas show the 68% credibility bands of the inferences. Blue points show the correlation parameter of galaxies in the MS selected without assuming a flux limit. 
After the science verification, we fit the model to the observed G3L signal of KV450×GAMA. The bestfitting model parameters, χ^{2}values and plots are reported in Table C.3. For the observations, all pvalues exceed 0.05, indicating no significant deviation between fit and measurement at a 95% CL. Again, we perform a KStest on the cumulative distribution of pvalues to evaluate the overall agreement of the model and measurements (dashed lines in Fig. 11). The KS distance of this distribution to a uniform distribution is now 0.31, larger than for the science verification but still consistent with a uniform distribution at the 95% CL.
The HODs inferred for KV450×GAMA are shown in Fig. 12b by solid lines and 68% credibility regions, again using the combinations m1–m5, m2–m5, m3–m5, and m4–m5; other combinations for the same sample yield consistent results within the 68% CI. The HODs vary with the stellar mass of the lenses. In contrast to massive halos, low mass halos are mainly populated by galaxies of low stellar mass. For example, in halos with masses below 2 × 10^{11} M_{⊙}, m1 galaxies are the most numerous, whereas m4 galaxies dominate halos with masses above 4 × 10^{14} M_{⊙}. The most massive m5 galaxies, however, are never the most numerous sample between 10^{11} to 10^{15} M_{⊙}. Instead, they are outnumbered by satellites of lower stellar mass. Compared to the SAM HOD (data points), the HODs of KV450×GAMA are consistent with the model.
Each sample HOD is an average over the sample redshift distribution, shown in Fig. 5. Therefore, the inferred HODs are affected by the survey flux limit. Most affected are the faint m1galaxies. Although they are observed up to z ∼ 0.34, for z ≳ 0.15 only the most massive m1 galaxies are seen (M_{*} > 2 × 10^{9} h^{−2} M_{⊙}). Due to this flux selection, the inferred average HOD is skewed towards the HOD of the more massive galaxies in this stellarmass bin. Likewise, the HODS of the other samples also give more weight to the most massive galaxies in the samples but with less bias than for the faint m1 sample.
The threshold masses increase with stellar mass, indicating that galaxies with higher stellar mass prefer to inhabit more massive halos. We show this trend in Fig. 14, which plots against the average stellar mass of the samples. For each stellar mass bin, there are four estimates of , from the four possible sample combinations, slightly displaced along the yaxis. As noted above, these four estimates are consistent with each other for all samples.
Fig. 14. Threshold mass measured for the GAMA galaxies as a function of the average stellar mass of each stellar mass bin (Green crosses). We show all four estimates for for each sample a, slightly displaced along the yaxis for visibility. Horizontal errors correspond to 68% CI of ; vertical errors show the standard deviation of the stellar masses of galaxies within a sample. 
The central galaxies parameter α^{a} also increases with stellar mass: The central galaxy of a halo with is more likely to be a than b, if a is a sample with higher stellar mass than b. As the four α^{a}s obtained for each stellar mass bin are consistent, we can average them to a for each sample a, yielding , , , , . The sum of these α^{a} is , which is consistent with unity within a 2σ CI. To derive the uncertainties, we assumed that each individual estimate of α^{a} is an independent measurement. This assumption is not true in our case (the data vectors for the ten fits contain multiple times the same measurements), so the derived uncertainties are probably larger in reality. A better estimate of the uncertainty on requires a simultaneous fit of all 10 combinations of the five mass bins. This is, unfortunately, not feasible in our case because our jackknifederived covariance becomes singular for such a large data vector.
The satellite parameters σ^{a}, f^{a}, and β^{a} do not change significantly with stellar mass (68% CI), which is not surprising for σ^{a} because we cannot constrain σ^{a} better than the prior range. The concentration parameter is consistent with f^{a} ≈ 1 for all samples (68% CI), so there is generally no detectable deviation of the satellite distribution from the matter distribution inside halos. The slope of the satellite numbers is β^{a} ≈ 1 for all mass samples.
Lastly, we show the correlation r^{ab} of GAMA satellite numbers in Fig. 13 (green, dotted lines). There is a similar trend to the simulated galaxies (black data points), with a positive r^{ab}, approximately scaling linearly with halo mass (ϵ^{ab} ≈ 1). For the sample combinations m1–m2, m1–m3, m2–m3, and m2–m4, we find A^{ab} > 0 (95% CI): the satellite numbers of galaxies below stellar mass of 10^{10.5} h^{−2} M_{⊙} are positively correlated. Since r^{ab} increases with halo mass, the correlations between low stellar masses become relevant beyond galaxygroup size halos, where , , and for m = 10^{13} M_{⊙} and ϵ = 1 (68% CI).
We remind here that r^{ab} > 1 are possible because the r^{ab} are a Pearson correlation coefficient, hereafter, only for a Poissonian variance of satellite numbers. In particular, we have the relation
so that values of r^{ab} > 1 indicate a superPoissonian variance. For the strongly nonPoissonian m1 and m2 SAM galaxies at m = 10^{14} M_{⊙}, r^{m1 m2} = 1.2, which corresponds to (see Fig. B.1 in Appendix B for the deviation from a Poisson statistic). For the samples m4 and m5, which are closer to Poissonian, the r^{m4 m5} = 0.07 at the same halo mass corresponds to . For red and blue SAM galaxies, the measured r^{red blue} = 3.2 at m = 10^{14} M_{⊙} corresponds to .
The positive correlation coefficients are partly caused by the flux limit of the survey. To see this effect, consider as an example the correlation between m1 and m2 galaxies. At low z ≲ 0.1, all galaxies with stellar masses in the m1 range are observed, so the correlation parameter at these redshifts measures the true correlation between all m1 and m2 galaxies. However, at higher redshifts, only m1 galaxies at the high end of the stellar mass bin are observed due to the flux limit. These galaxies have stellar masses closer to m2galaxies. They are therefore more similar to m2 galaxies and stronger correlated than the overall less massive m1 galaxies. This systematically increases the inferred value of r^{ab}, which is an average over the correlations at different redshifts, compared to the correlation without flux limit.
To test whether the flux limit is the only cause of th e correlation, we consider samples of galaxies without a flux limit in the MS. For this purpose, we select all galaxies up to z = 0.5, irrespective of their magnitude, and divide them by the same stellarmass cuts as the fluxlimited samples. We then directly estimate the parameter r^{ab} for these samples, shown in grey in Fig. 13. The r^{ab} of the unlimited samples are generally lower than for the fluxlimited samples, indicating that the flux limit indeed causes a higher correlation. Nevertheless, the r^{ab} are still clearly positive for the sample combinations m1–m2, m1–m3, m2–m3, m2–m4, and m3–m4 without flux limit. Consequently, the measured positive correlation is not purely an effect of the survey incompleteness.
7. Discussion
We presented a new method to measure galaxy HODs and the correlation of perhalo satellite galaxy numbers using G3L, a weak gravitational lensing effect measuring the projected galaxygalaxymatter bispectrum. To this end, we constructed a new halo model for the galaxygalaxymatter bispectrum and the G3L aperture statistics. We validated the science analysis for different selections of lens galaxies and demonstrated that it accurately recovers the true HODs within 68% errors for a survey with an area of ≃10^{3}deg^{2} and lens galaxies between 0 ≤ z ≤ 0.5 and brighter than r = 19.8 mag. Therefore, our G3L analysis is accurate enough to infer HODs from KV450×GAMA, which has a smaller footprint (180 deg^{2}) and similar selections for the lens galaxies.
We inferred HODs for GAMA for galaxies in two colour samples (red and blue) and five stellar mass bins between 10^{8.5}h^{−2} M_{⊙} and 10^{11.5}h^{−2} M_{⊙} (m1 to m5). The bestconstrained HOD parameter for the KV450 × GAMA is the threshold halo mass , which is M_{th} ≈ 10^{11} M_{⊙} for our blue galaxies and M_{th} ≈ 10^{12} M_{⊙} for our red galaxies. These values match the expectation that red galaxies are typically group or cluster galaxies, whereas blue galaxies tend to be field galaxies. At lower halo masses, we find tight upper limits on ⟨Nm⟩, indicating that such halos rarely host galaxies satisfying our selection criteria. The transition region between this regime and halos typically containing galaxies is only poorly constrained. These poor constraints are reflected by the large uncertainty on the transition parameter σ, even for the mock data in our science verification. Therefore, G3L alone cannot constrain σ better than our prior, and future studies may fix σ to a fiducial value or dispense with this parameter altogether by modelling ⟨m⟩ by a step function.
The 1halo term of the G3L signal of redred lens pairs stretches to larger scales than for blueblue lens pairs (Fig. 7). This finding can be explained by the tendency of red galaxies to populate more massive halos than blue galaxies. Since massive halos are larger than less massive ones, pairs of red galaxies inside the same halo can have wider separations than pairs of blue galaxies. Accordingly, the 1halo term extends to larger scales for redred lens pairs than for blueblue ones. However, mixed redblue pairs exist in intermediate halos, large enough to host red galaxies but small enough to contain a significant fraction of blue galaxies. Consequently, the crossover between the domination of the 1halo and the 3halo term occurs at larger scales than for blueblue lens pairs.
In contrast to previous HOD studies, we measured the correlation of numbers of satellites inside halos between galaxy samples. We report a > 3σ detection of a positive correlation for red and blue GAMA galaxies ( for m ≳ 10^{13} M_{⊙}, rising towards galaxy cluster scales). Similar positive correlations are present between the samples m1 (stellar masses below 10^{9.5} M_{⊙}) and m2 or m3 (stellar masses between 10^{9.5} M_{⊙} and 10^{10.5} M_{⊙}), as well as between sample m2 (stellar masses between 10^{9.5} M_{⊙} and 10^{10} M_{⊙}) and samples m3 or m4 (stellar masses between 10^{10} M_{⊙} and 10^{11} M_{⊙}). In particular, galaxies with similar stellar masses (for example from neighbouring samples m1 and m2) show positive correlations. Towards smaller halos, correlations become irrelevant due to the almost linear decrease of r^{ab}, visible both in the SAMs and in our observational data (Figs. 9 and 13). This trend fits our toy model consideration in Sect. 3.1, where a decreasing correlation is the consequence of Poisson shot noise inside lowoccupancy halos. This finding implies that the assumption of uncorrelated satellite distributions is probably still appropriate if only galaxies in lowmass halos, m ≲ 10^{13} M_{⊙}, are considered.
The obvious presence of correlations, especially for halos of the group and clustermass scale, questions the assumption of no correlation in halo models for the crosscorrelation statistics between galaxy samples in galaxyclustering studies (e.g. Scranton 2001, 2002; Zehavi et al. 2005; Simon et al. 2008). The correlations also raise questions about their origin and impact on galaxy models. To address the first, we note that finding strong correlations for similar galaxy samples is not unexpected. If two galaxy samples are drawn randomly from a pool of galaxies, they are statistically identical and their satellite numbers inside halos differ only by Poisson shot noise. Consequently, the satellite numbers are tightly correlated, which could explain the correlations between satellite galaxies in neighbouring stellar mass bins (Fig. 13). Stellar mass varies continuously and does not uniquely define a category of galaxies. Therefore, the division of galaxies by their stellar mass is ultimately arbitrary and satellite galaxies from neighbouring stellar mass bins have similar statistical properties. In addition, errors in the estimators for stellar mass blur the stellar mass bin on the edges, similar to randomly drawing galaxies from the same pool. This effect, in combination with depth variations between the samples due to the survey flux limit (Sect. 6.2), systematically increases r^{ab}. Moreover, strong correlations also emerge if two samples have a (close to) fixed ratio of satellite numbers inside a halo. This effect could occur for red and blue galaxies for clustersized halos (m ≳ 10^{14} M_{⊙}): Galaxy models commonly assume that starforming blue galaxies falling into a galaxy cluster are quenched and turn into red galaxies with a certain probability. If this probability is roughly constant for halos of similar mass, satellite numbers of red and blue galaxies would be strongly correlated. Whatever the cause for the correlations for the GAMA galaxies, the consistency with the mock galaxies shows that galaxy SAMs already account for it at z < 0.5. Therefore, our findings do not hint at a need to improve these models at low redshift. However, as the galaxy population in clusters evolves with time, it will be interesting to study the correlations at higher redshift in the future.
It could be interesting to use G3L to investigate the evolution of the crosscorrelation parameter r^{ab} with cosmic time, by selecting similar galaxy samples at different redshifts. Such a selection is difficult with the GAMA sample as the detection limit of r = 19.8 mag restricts our analysis to z ≲ 0.5. An analysis of the redshift evolution of the HOD parameters would consequently require a deeper lens galaxy sample, so lenses could be divided into different redshift bins. It would also require a large number of source galaxies with higher redshifts than the lens samples so that the source number density is sufficient for a significant G3L signal.
Our G3L halo model makes several unrealistic assumptions, although they are sufficient for analysing KV450 × GAMA. Concerning the matter distribution, we ignore halo exclusion (Smith et al. 2007), the dependency of halo properties on environment and assembly history (assembly bias, Gao & White 2007), or galactic subhalos, known to be relevant for galaxygalaxy lensing at subarcmin scales (e.g. Velander et al. 2014). Studies of simulated galaxies from SAMs (e.g. Zehavi et al. 2018) and hydrodynamical simulations (e.g. Hadzhiyska et al. 2020) have shown that HODs also vary with halo formation time and environment. As we only include a dependence on halo mass in our HOD model, our parameter values can be considered an average over all secondary halo properties. Furthermore, concerning the lens galaxy distribution, there is a known feature of the HOD of the simulated galaxies our model does not include: a local minimum in the HODs of galaxies with stellar masses below 10^{10.5} h^{−2} M_{⊙} (see Fig. 12b). This minimum occurs at the transition region between the domination of the HODs by central and satellite galaxies. There, contrary to our model for ⟨m⟩ the mean number of centrals decreases again with halo mass for galaxies from samples m1, m2, and m3 (stellar masses below 10^{10.5} h^{−2} M_{⊙}). A similar observation was reported by Berlind et al. (2003) and Zheng et al. (2005) for galaxies with young stellar populations: Due to quenching, massive halos are unlikely to have central galaxies with young stars suppressing ⟨m⟩ above a certain halo mass. Another case are emission line galaxies (Guo et al. 2019) where both the central and satellite number decreases with mass for massive halos. Therefore, our specific model implementation is not directly applicable to other lens sample selections, and the HOD model needs to be adjusted on a casebycase basis. Another issue is the assumption of a Poissonian variance for satellite numbers. Contrary to this assumption, the SAM galaxies by Henriques et al. (2015) have a subPoissonian variance for lowmass halos and superPoissonian variance for highmass halos, as we show in Fig. B.1. We investigated the impact of this effect by comparing the bestfitting G3L model to the prediction by a modified model using the actual variance ⟨⟩ of the SAM galaxies. The difference between the two G3L signals is smaller than the statistical errors in the verification data, so we can neglect nonPoissonian variances here (see Appendix B). However, a G3L analysis in future surveys may require a more sophisticated model with additional parameters for the satellite number variance, such as the models in Dvornik et al. (2018). Relaxing the Poisson assumption is a new test for galaxy models with G3L. While some simulations assume Poisson satellites (Kravtsov et al. 2004; Zheng et al. 2005), recent studies suggest otherwise (Dvornik et al. 2018; Gruen et al. 2018). Despite the above model approximations, our inferred HODs agree with the true HODs of the validation data within the 68% CI. Therefore, these approximations are acceptable for surveys at the same level (or worse) of statistical estimator noise as our validation data.
Aside from increasing the survey area or the source density in a new survey, better HOD constraints are achievable by tapping into two more sources of information in the observational data – after a future revision of the model and the statistical analysis. First, the lens redshifts can be exploited in the redshiftweighting scheme of lens pairs by Linke et al. (2020b). In this scheme, lens galaxy pairs that appear close in projection but are physically distant are downweighted compared to physically close lens pairs in estimating the G3L correlation function. This weighting increases the signaltonoise by up to 35% for lenses with spectroscopic redshifts. However, incorporating the scheme into the halo model requires abandoning the Limber approximation and revising Eq. (12). Without the Limber approximation, the projected bispectrum includes three lineofsight integrals over Bessel functions, whose evaluation is computationally challenging. There are several efforts to obtain projected matter power and bispectra without relying on the Limber approximation (Assassi et al. 2017; Campagne et al. 2017; Deshpande & Kitching 2020), which might be adopted for a weighted G3L. Second, we expect the parameter constraints to improve by combining G3L with other statistics, breaking degeneracies between HOD parameters (Sect. 3.3). For example, secondorder galaxygalaxy lensing or the mean number density of the lens galaxies are easily available from the same observational data. In particular, the mean number density can be measured to high precision, and it depends directly on ⟨m⟩+⟨m⟩, while G3L alone depends on the secondorder moments of the central and satellite galaxy numbers. The overall benefit of this combination is unclear because it requires a refinement of the statistical analysis in order to account for the (expected) covariance between the different probes. Focusing here purely on constraints with G3L, we also leave this promising prospect to a future analysis.
In conclusion, G3L with our halo model is a viable new method to infer the HODs of galaxies and the correlation between the halo distributions of galaxies from different populations. This information is useful for several other analyses. First, HODs offer a fast method to obtain realistic mock galaxy catalogues from halo catalogues (e.g. Carretero et al. 2015; Avila et al. 2018), which is faster than using SAMs or hydrodynamical simulations. Realistic mock galaxies are vital to validate the redshift calibration (e.g. van den Busch et al. 2020; Hildebrandt et al. 2021) or inference pipelines (e.g. Ferrero et al. 2021; DeRose et al. 2022). Including the crosscorrelation of galaxy samples constrained with G3L will increase the realism of the mock galaxy distribution and, therefore, the robustness of the analyses based on these simulated data. Second, HODs can be used to provide a physical interpretation of the galaxy bias, which tells us how galaxies trace the cosmic largescale structure (Cacciato et al. 2012; Simon & Hilbert 2018; Dvornik et al. 2018). Cosmological studies often use simple linear or quadratic models for galaxy biasing (Joachimi et al. 2021; Krause et al. 2021). They, therefore, need to exclude small scales from their analysis, where the simple models are not sufficiently accurate. HODs may provide a more accurate model on small scales so that cosmological parameters could be inferred from a larger scale range. However, the HODs obtained in this work depend on the assumed fiducial cosmology. Accordingly, using the G3L halo model for cosmological inference requires simultaneously constraining the HOD and cosmological parameters. This inference might be unfeasible in practice due to the highdimensional parameter space. Third, HODs are important tools to constrain the stellartohalo mass ratio, traditionally obtained from secondorder statistics, such as galaxygalaxy lensing (Velander et al. 2014; van Uitert et al. 2018; Dvornik et al. 2020). Combining these measurements with the HOD constraints from G3L could lead to tighter constraints on the relationship between baryonic and dark matter.
In their original derivation, Navarro et al. (1996) defined the mass by a sphere of enclosed mean density 200 ρ_{crit}. However, the fitting function for the concentration by Bullock et al. (2001), employed here, defines m in terms of the mean cosmic density .
The raytracing algorithm introduces a smoothing, which lowers the lensing power spectrum for large ℓ; for ℓ < 2 × 10^{4}, however, the bias is below 5% (Hilbert et al. 2009). The smoothing bias beyond ℓ ∼ 10^{4} is almost irrelevant for the G3L aperture statistics, as the filter function in Eq. (17) decreases sharply for high ℓ: For a small aperture of , we have .
Acknowledgments
We thank Elisa Chisari for her helpful comments and for acting as internal referee for the KiDS Collaboration. This work has been supported by the Deutsche Forschungsgemeinschaft through the project SCHN 342/151. LL received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. AHW is supported by a European Research Council Consolidator Grant (No. 770935). Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A3016, 177.A3017, 177.A3018, 179.A2004, 298.A5015. We also use products from the GAMA survey. GAMA is a joint EuropeanAustralasian project based around a spectroscopic campaign using the AngloAustralian 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 several independent survey programmes, including GALEX MIS, VST KiDS, VISTA VIKING, WISE, HerschelATLAS, 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.gamasurvey.org/. Author contributions. All authors contributed to the development and writing of this paper. The authorship list is given in two groups: The lead authors (LL, PSi, PS), followed by an alphabetical list of contributors to either the scientific analysis or the data products.
References
 Anderson, T. W. 2003, An Introduction to Multivariate Statistical Analysis (WileyInterscience) [Google Scholar]
 Assassi, V., Simonović, M., & Zaldarriaga, M. 2017, J. Cosmol. Astropart. Phys., 2017, 054 [CrossRef] [Google Scholar]
 Avila, S., Crocce, M., Ross, A. J., et al. 2018, MNRAS, 479, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [Google Scholar]
 Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587 [Google Scholar]
 Berlind, A. A., Weinberg, D. H., Benson, A. J., et al. 2003, ApJ, 593, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1 [Google Scholar]
 Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
 Cacciato, M., Lahav, O., van den Bosch, F. C., Hoekstra, H., & Dekel, A. 2012, MNRAS, 426, 566 [Google Scholar]
 Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Campagne, J. E., Neveu, J., & Plaszczynski, S. 2017, A&A, 602, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carretero, J., Castander, F. J., Gaztañaga, E., Crocce, M., & Fosalba, P. 2015, MNRAS, 447, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
 Clampitt, J., Miyatake, H., Jain, B., & Takada, M. 2016, MNRAS, 457, 2391 [NASA ADS] [CrossRef] [Google Scholar]
 Clampitt, J., Sánchez, C., Kwan, J., et al. 2017, MNRAS, 465, 4204 [NASA ADS] [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [Google Scholar]
 Crittenden, R. G., Natarajan, P., Pen, U.L., & Theuns, T. 2002, ApJ, 568, 20 [NASA ADS] [CrossRef] [Google Scholar]
 de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2022, Phys Rev. D, 105, 123520 [NASA ADS] [CrossRef] [Google Scholar]
 Deshpande, A. C., & Kitching, T. D. 2020, Phys Rev. D, 101, 103531 [NASA ADS] [CrossRef] [Google Scholar]
 Driver, S. P., Norberg, P., Baldry, I. K., et al. 2009, Astron. Geophys., 50, 12 [Google Scholar]
 Dvornik, A., Hoekstra, H., Kuijken, K., et al. 2018, MNRAS, 479, 1240 [Google Scholar]
 Dvornik, A., Hoekstra, H., Kuijken, K., et al. 2020, A&A, 642, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32 [NASA ADS] [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [Google Scholar]
 Erben, T., Schirmer, M., Dietrich, J. P., et al. 2005, Astron. Nachr., 326, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Farrow, D. J., Cole, S., Norberg, P., et al. 2015, MNRAS, 454, 2120 [Google Scholar]
 Ferrero, I., Crocce, M., Tutusaus, I., et al. 2021, A&A, 656, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gao, L., & White, S. D. M. 2007, MNRAS, 377, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, B. 2009, GNU Scientific Library Reference Manual– Third Edition (Network Theory Ltd.) [Google Scholar]
 Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys Rev. D, 98, 023507 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, H., Yang, X., Raichoor, A., et al. 2019, ApJ, 871, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Hadzhiyska, B., Bose, S., Eisenstein, D., Hernquist, L., & Spergel, D. N. 2020, MNRAS, 493, 5506 [NASA ADS] [CrossRef] [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [Google Scholar]
 Hilbert, S., Hartlap, J., White, S. D. M., & Schneider, P. 2009, A&A, 499, 31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124 [EDP Sciences] [Google Scholar]
 Ishikawa, S., Okumura, T., Oguri, M., & Lin, S.C. 2021, ApJ, 922, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338 [Google Scholar]
 Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaiser, N. 1992, ApJ, 388, 272 [Google Scholar]
 Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krause, E., Fang, X., Pandey, S., et al. 2021, ArXiv eprints [arXiv:2105.13548]. [Google Scholar]
 Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35 [Google Scholar]
 Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500 [Google Scholar]
 Linke, L., Simon, P., Schneider, P., et al. 2020a, A&A, 640, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Linke, L., Simon, P., Schneider, P., & Hilbert, S. 2020b, A&A, 634, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liske, J., Baldry, I. K., Driver, S. P., et al. 2015, MNRAS, 452, 2087 [Google Scholar]
 Liu, J. S. 2004, Monte Carlo Strategies in Scientific Computing, 1st edn. (New York, NY: Springer), 31 [CrossRef] [Google Scholar]
 Mandelbaum, R., Hirata, C. M., Broderick, T., Seljak, U., & Brinkmann, J. 2006, MNRAS, 370, 1008 [Google Scholar]
 Maraston, C. 2005, MNRAS, 362, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, S. M. 2019, Ph.D. Thesis, University of Bonn, Germany [Google Scholar]
 Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858 [Google Scholar]
 Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [Google Scholar]
 Nakamura, T. T., & Suto, Y. 1997, Prog. Theor. Phys., 97, 49 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rödiger, J. 2009, Ph.D. Thesis, University of Bonn, Germany [Google Scholar]
 Ross, A. J., & Brunner, R. J. 2009, MNRAS, 399, 878 [NASA ADS] [CrossRef] [Google Scholar]
 Saghiha, H., Simon, P., Schneider, P., & Hilbert, S. 2017, A&A, 601, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schirmer, M. 2013, ApJS, 209, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Watts, P. 2005, A&A, 432, 783 [CrossRef] [EDP Sciences] [Google Scholar]
 Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Scranton, R. 2001, MNRAS, 332, 697 [Google Scholar]
 Scranton, R. 2002, MNRAS, 339, 410 [Google Scholar]
 Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]
 Simon, P., & Hilbert, S. 2018, A&A, 613, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, P., Watts, P., Schneider, P., et al. 2008, A&A, 479, 655 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, P., Hetterscheidt, M., Wolf, C., et al. 2009, MNRAS, 398, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, P., Erben, T., Schneider, P., et al. 2013, MNRAS, 430, 2476 [Google Scholar]
 Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys Rev. D, 75 [CrossRef] [Google Scholar]
 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]
 Szapudi, I., & Szalay, A. S. 1998, ApJ, 494, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587 [Google Scholar]
 Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22 [NASA ADS] [CrossRef] [Google Scholar]
 van den Busch, J. L., Hildebrandt, H., Wright, A. H., et al. 2020, A&A, 642, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662 [NASA ADS] [CrossRef] [Google Scholar]
 Velander, M., van Uitert, E., Hoekstra, H., et al. 2014, MNRAS, 437, 2111 [Google Scholar]
 Venemans, B. P., Verdoes Kleijn, G. A., Mwebaze, J., et al. 2015, MNRAS, 453, 2259 [Google Scholar]
 Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nat. Rev. Phys., 2, 42 [Google Scholar]
 Wang, Y., Yang, X., Mo, H. J., & van den Bosch, F. C. 2007, ApJ, 664, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Watts, P., & Schneider, P. 2005, in Gravitational Lensing Impact on Cosmology, eds. Y. Mellier, & G. Meylan, IAU Symp., 225, 243 [NASA ADS] [Google Scholar]
 Weisstein, E. W. 2022, Delta Function. From MathWorld–A Wolfram Web Resource, http://mathworld.wolfram.com/DeltaFunction.html, Last visited on 02/3/2022 [Google Scholar]
 White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
 Wright, A. H., Robotham, A. S. G., Bourne, N., et al. 2016, MNRAS, 460, 765 [Google Scholar]
 Wright, A. H., Hildebrandt, H., Kuijken, K., et al. 2019, A&A, 632, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2005, ApJ, 630, 1 [Google Scholar]
 Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Zehavi, I., Contreras, S., Padilla, N., et al. 2018, ApJ, 853, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760 [Google Scholar]
Appendix A: Calculation of galaxygalaxymatter bispectrum
Here we derive the galaxygalaxymatter bispectrum for mixed (a ≠ b) and unmixed (a = b) galaxy pairs in a halo model. While the bispectrum for unmixed pairs has already been derived partially in Watts et al. (2005), Rödiger (2009), and Martin (2019), a model for mixed pairs is, to our knowledge, still pending and therefore shown in detail here.
Our model assumes that all statistics of matter and galaxies inside halos only depend on halo mass, halo positions are clustered according to a deterministic linear bias, matter and galaxy number density profiles are spherically symmetric, satellite galaxies have no subhalos, and galaxy positions inside halos are statistically independent of each other. These approximations greatly simplify the model formalism while still being sufficiently accurate to describe the G3L aperture statistics in the science verification data. For brevity, we drop all explicit redshift dependencies, and all quantities are understood to be evaluated at the same redshift.
The halo model averages over the probability density functions (PDFs) of the halos masses, positions, galaxy numbers and galaxy positions. The average of an arbitrary quantity f(…) depending on the halo variables is, for H halos distributed in the volume V and containing two distinct galaxy samples a ≠ b,
For two identical galaxy samples a = b, this expression simplifies to
where is the joint probability for halo h to have central and satellite galaxies from sample a, and is the PDF for the positions of these satellites relative to the halo centre. The joint PDF P_{m} of halo masses, assuming that halo masses are independent of each other, is a product of the HMF,
where is the mean halo number density^{4}. The PDF is a product of the normalised spatial distributions of satellites in each halo,
under the assumption of independent satellite positions inside a halo. The joint probability P_{N} of perhalo galaxy numbers has the moments
where p, q, r, s ∈ ℕ_{0}.
After inserting Eqs. (49) and (51) into (46) and using Eq. (A.1), the bispectrum for mixed pairs, a ≠ b, is given by
where unconnected terms do not contain δ_{D}(k_{1} + k_{2} + k_{2}), is the mean galaxy number density and is the mean matter density (as in Eqs. 1 and 2). For unmixed pairs, a = b, the bispectrum is given by
where is the Kronecker symbol.
We divide the sum over i, j and k into the 1halo term _{1} for i = j = k, the 2halo term _{2} for i = j ≠ k, i = k ≠ j, and j = k ≠ i, and the 3halo term _{3} for i ≠ j ≠ k. In the following, we calculate these terms individually.
A.1. 1halo term
To calculate the 1halo term _{1} with i = j = k, we distinguish between mixed (a ≠ b) and unmixed (a = b) pairs. For a = b, the 1halo term is given by
Evaluating all mintegrals aside from the one over m_{i} (H − 1 integrals in total) leads to
Using Eqs. (A.6), ∫ⅆ[3]x ug^{a}(x  m) = 1, and , we find
explicitly using the secondorder moments of P_{N} in the last line.
We consider large volumes V here, rendering the Vintegrals over the halo positions x with exp−ik ⋅ x in the integrand essentially Fourier transforms. This can, for example, be seen for a cubic volume V_{s} with side length s,
where the k_{i} are the components of k. For large V_{s} (s → ∞), this expression approximates a Diracdistribution (e.g. Weisstein 2022), as
and hence
In this limit (assumed henceforth),
which leads to
with k_{3} = −k_{1} − k_{2} and
For mixed pairs, a ≠ b, an analogous calculation gives
In the main text, we denote _{1}(k_{1}, k_{2}, k_{3}) as _{1}(k_{1}, k_{2}) for brevity.
A.2. 2halo term
The 2halo term _{2} in Eq. (A.7) contains the contributions from i = j ≠ k, i = k ≠ j, and k = j ≠ i. For unmixed pairs, a = b, the 2halo term is given by
We now evaluate all mintegrals independent of m_{i} and m_{j} (these are H(H − 1) terms in total), use Eq. (A.6), and evaluate the Δx integrals over the halo profiles . This leads to
In the limit of an infinite volume V, Eq. (A.14), the halo twopoint correlation function ξ_{h} and its Fourier transform P_{h}, gives
With this expression, neglecting unconnected terms proportional to δ_{D}(k_{1}) and δ_{D}(k_{2}) and the approximation H(H − 1)/H^{2} ≃ 1, we find
where again k_{3} = −k_{1} − k_{2} and
For mixed pairs, a ≠ b, a similar calculation yields
In the main text we just denote _{2}(k_{1}, k_{2}).
A.3. 3halo term
The 3halo term for unmixed pairs a = b is given by
We evaluate all mintegrals independent of m_{i}, m_{j}, and m_{k} (H(H − 1)(H − 2) terms in total), and evaluate the Δx integrals. This leads to
Now, again in the limit of an infinite volume V, Eq. (A.14), we use the two and threepoint functions ξ_{h} and ζ_{h} of halo clustering and their Fourier transforms P_{h} and B_{h} to find
As before, all terms not proportional to δ_{D}(k_{1} + k_{2} + k_{3}) are unconnected and therefore neglected. Using this expression and the approximation H(H − 1)(H − 2)/H^{3} ≃ 1, we obtain
For mixed pairs a ≠ b a similar calculation yields
Again, we just denote _{3}(k_{1}, k_{2}) in the main text.
Appendix B: Poissonianity of satellite galaxies
Our HOD model explicitly assumes Poisson satellites – satellite numbers inside halos that vary according to a Poisson statistic. In this section, we explore the bias due to this assumption in the presence of nonPoisson satellites. To some degree, deviations from Poisson satellites are indeed visible in our SAM galaxies: Figure B.1 shows the ratio of σ(m) to the Poisson variance for different samples of the simulated galaxies. Poisson satellites, where this ratio is unity, are present in the intermediate halomass range, such as between 2 × 10^{11} M_{⊙} ≲ m ≲ 10^{12} M_{⊙} for galaxies from stellarmass bin m2. On the lowmass end, satellites are distributed subPoissonian (ratio is below unity), and on the highmass end, they are superPoissonian (ratio exceeds unity). The trend is similar across all galaxy populations, although with a shift of the intermediate range depending on the stellar mass of the sample.
Fig. B.1. Ratio of the standard deviation of the galaxy satellite number per halo and the Poisson variance. For Poissonian satellites this ratio is unity. Left: Ratio for all simulated galaxies (black crosses), red simulated galaxies (red crosses), and blue simulated galaxies (blue dots). Right: Ratio for stellar massselected samples. 
We estimate the bias induced by the mismatch between the true satellite variance and the model assumption by computing the fractional change of ⟨𝒩^{a}𝒩^{b}ℳ⟩_{mod} when switching from Poisson satellites to the true variance in the SAM. This means: For each best fit to the simulated G3L signal, we calculate new model predictions ⟨𝒩^{a}𝒩^{b}ℳ⟩_{mod}, using the HOD parameters from the best fit but updating the number of satellite pairs inside a halo to
where σ^{2}(m) is the true satellite variance in the simulation. Clearly, for , that is, Poisson satellites, this reduces to Eq. (39).
Figure B.2 shows the fractional change between the updated model and the original best fit with Poisson satellites for the red and blue galaxies. Also shown is the fractional difference of the measured aperture statistics to the best fit. For red galaxies, the updated model differs from the original strongest at 4′ by roughly 50%. For blue galaxies, the difference is greatest at 2′ and has a similar magnitude. It is to be expected that the aperture statistics at scales between 1′and 10′are most affected because the 1halo term, containing ⟨( − 1)m⟩, dominates here. Smaller scales are dominated by the halo term containing ⟨m⟩, whereas larger scales are dominated by the 3halo term. Although a bias is visible by the solid line in Fig. B.2, it is of the order of or smaller than the uncertainties on the aperture statistics measurements (error bars). In particular, for blue galaxies, the measurements cannot discriminate between the two models. Therefore, while the presence of nonPoisson satellites affects the model prediction, a Poisson model is accurate enough for the G3L analysis of measurements in this work.
Fig. B.2. Fractional difference of aperture statistics measurement ⟨𝒩^{a}𝒩^{b}ℳ⟩_{meas} (points) and modified aperture statistics model ⟨𝒩^{a}𝒩^{b}ℳ⟩_{mod}(lines) to the original model ⟨𝒩^{a}𝒩^{b}ℳ⟩ for redred lens pairs (left) and blueblue lens pairs (right) in the MS. The original model assumes a Poissonian satellite distribution, while the modified model uses the directly measured distribution of satellite galaxies. 
We have repeated this test also for the satellite galaxies from the stellar massselected samples. Galaxies from samples m4 and m5 are closer to a Poisson satellite model than blue galaxies and therefore show less bias (below 10%). For the other samples, the bias is larger (up to 61% for m1, 57% for m2, 19% for m3 satellites), but since the measured aperture statistics have larger uncertainties, the bias is of even less relevance than for red and blue galaxies.
Finally, a deviation of the variance in satellite numbers from a Poisson statistic also affects the interpretation of the correlation parameter r^{ab}, as defined in Eq. (40). For Poisson satellites, r^{ab} corresponds to a Pearson coefficient, that is, r^{ab} = 1 for perfectly correlated and , r^{ab} = −1 for perfectly anticorrelated satellite numbers, and r^{ab} = 0 for uncorrelated samples. However, for a superPoisson variance in the highmass regime, our r^{ab} is larger than the Pearson coefficient. Conversely, for lowmass halos with subPoisson variance, r^{ab} is smaller than the Pearson coefficient. The combination of these trends increases the slope ϵ^{ab} of our r^{ab}(m) compared to a Pearson definition of correlation. Nevertheless, the sign of r^{ab} remains unchanged, and if the samples a and b are uncorrelated, both r^{ab} and the Pearson coefficient vanish.
Appendix C: Results of model fit to stellar massselected galaxies
In Fig C.1 we show the model fits to the G3L aperture statistics from stellar massselected lenses in the MS and the χ^{2} of the fits. All fits have 76 d.o.f. (the same as the fit to the colourselected samples in Sect. 6.1), so a χ^{2} < 97.4 signifies an agreement between fit and data at the 95% CL. The HOD parameters of the bestfitting models are given in Table C.1; the parameters A^{ab} and ϵ^{ab} that determine the crosscorrelation between satellite numbers together with the χ^{2} and pvalues of the fits are given in Table C.2. Figure C.4 shows the model fits to the G3L aperture statistics from stellar massselected lenses in KV450× GAMA and the χ^{2} of the fits. The parameters of the bestfitting models are given in Tables C.3 and C.4.
Fig. C.1. G3L measurement in MS (points) and bestfitting halo model (lines) for stellar massselected lens samples, as defined in Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. Each row was fitted individually, leading to the χ^{2} values in the last panel. The corresponding halo model parameters are given in Table C.1. 
Fig. C.4. G3L measurement in KV450 × GAMA (points) and bestfitting halo model (lines) for stellar massselected lens samples, as defined in Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. Each row was fitted individually, leading to the χ^{2} values in the last panel. The corresponding halo model parameters are given in Table C.3 
best fit values for halo model parameters for stellarmassselected lenses in MS for each stellar mass sample a.
best fit values of HOD parameters describing satellite number crosscorrelation, and χ^{2} and pvalues for the G3L halo model fit to MS
best fit values for halo model parameters for stellarmassselected lenses in KV450 × GAMA for each stellar mass sample a.
best fit values of HOD parameters describing satellite number crosscorrelation, and χ^{2} and pvalues for the G3L halo model fit to KV450 × GAMA
All Tables
Bestfitting values of halo model parameters for colourselected lenses and 68% credibility intervals (d.o.f. = 76).
best fit values for halo model parameters for stellarmassselected lenses in MS for each stellar mass sample a.
best fit values of HOD parameters describing satellite number crosscorrelation, and χ^{2} and pvalues for the G3L halo model fit to MS
best fit values for halo model parameters for stellarmassselected lenses in KV450 × GAMA for each stellar mass sample a.
best fit values of HOD parameters describing satellite number crosscorrelation, and χ^{2} and pvalues for the G3L halo model fit to KV450 × GAMA
All Figures
Fig. 1. Geometry of a G3L configuration with one source and two lens galaxies on the sky; adapted from Schneider & Watts (2005). Lens galaxies are at angular positions θ_{1} = θ + ϑ_{1} and θ_{2} = θ + ϑ_{2} on the sky; the source galaxy is at θ. The angle between the sourcelens connections is the opening angle ϕ. 

In the text 
Fig. 2. Mean perhalo numbers of galaxies for fiducial HOD parameters in Table 1. The solid black line shows the total galaxy number per halo, the dashed red line shows the fraction of halos with central galaxies, and the dotted blue line shows the number of satellite galaxies per halo. 

In the text 
Fig. 3. Impact of halo model parameters on ⟨𝒩^{a}𝒩^{a}ℳ⟩ for unmixed lens pairs. In each panel, only one parameter is varied at a time. Solid lines indicate the total ⟨𝒩^{a}𝒩^{a}ℳ⟩, while dashed lines show the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term. 

In the text 
Fig. 4. Impact of the correlation of satellite numbers on the aperture statistics of mixed lens pairs. Satellite numbers are either fully correlated (r^{ab} = 1, blue lines), uncorrelated (r^{ab} = 0, black lines) or anticorrelated (r^{ab} = −1, red lines). Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term. 

In the text 
Fig. 5. Normalised redshift distributions n(z) of GAMA galaxies, selected by colour (left) and stellar mass (right). 

In the text 
Fig. 6. Stellarmass and redshift of GAMA galaxies, divided by colour (left) and stellar mass (right). 

In the text 
Fig. 7. G3L measurement (points) and bestfitting halo model (lines) for red and blue galaxies in the MS (upper plot) and the KV450 × GAMA (lower plot). Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. The left panels show the result for redred galaxy pairs, the central panels for blueblue galaxy pairs, and the right panels for redblue mixed pairs. Error bars correspond to the standard deviation estimated from the jackknife resampling detailed in Sect. 5.1. (a) MS. (b) KV450 × GAMA. 

In the text 
Fig. 8. Mean perhalo numbers of simulated galaxies (top) and observed galaxies (bottom) as function of halo mass. Red crosses (blue points) indicate the true HOD of simulated red (blue) galaxies, where the error bars are the standard deviation of the mean over the 64 lineofsights. The lines indicate the perhalo numbers inferred from the fit to the G3L signal for red (solid red) and blue galaxies (dashed blue). The shaded areas are the 68% credibility areas of the halo model fit. (a) MS. (b) KV450 × GAMA. 

In the text 
Fig. 9. Correlation parameter r^{ab} for red and blue galaxies in the simulation and observation as a function of halo mass. Black crosses show the direct estimate for the simulation, where the error bars are the standard deviation of the mean over the 64 lineofsights. The solid brown line shows the r^{ab} inferred by the halo model fit to the G3L signal of the MS, and the green dashed line is the result of the fit to the KV450 × GAMA G3L signal. The shaded areas show the 68% credibility bands of the fits. 

In the text 
Fig. 10. Correlation matrix for aperture statistics measurement in KV450×GAMA for red and blue galaxy samples. The data vector is ordered as given in Eq. (64), starting with the smallest aperture radius. 

In the text 
Fig. 11. Cumulative distribution of pvalues of G3L halo model fits for MS (orange, solid) and KV450 × GAMA (green, dashed). For a perfect description of G3L signal and data noise, the distributions would be consistent with a uniform distribution (black, dotted). 

In the text 
Fig. 12. Mean perhalo galaxy numbers in the simulation (top) and observation (bottom) for lens galaxies from each stellar mass bin as function of halo mass. Crosses indicate the directly estimated perhalo numbers of simulated galaxies, while lines show the predictions from the G3L fits. The shaded areas indicate the 68% confidence areas. Left panels: the mean perhalo numbers for galaxies from stellar mass bins m1, m3, and m5 obtained from the fits to the G3L signal for m1–m5, m3–m5, and m4–m5. The right panels show the same for galaxies from stellar mass bins m2 and m5, obtained from the fits to the signal for m2–m5 and m4–m5. The corresponding HOD parameters are listed in Table C.1. (a) MS. (b) KV450 × GAMA. 

In the text 
Fig. 13. Correlation parameter r^{ab} for stellar massselected galaxies in the MS and in KV450×GAMA. Black crosses show the true correlation for the MS, where the error bars are the standard deviation of the mean over the 64 simulated lineofsights. The solid brown line shows r^{ab} inferred from a simulated G3L analysis of MS, and the green dashed line the inference for KV450×GAMA. The shaded areas show the 68% credibility bands of the inferences. Blue points show the correlation parameter of galaxies in the MS selected without assuming a flux limit. 

In the text 
Fig. 14. Threshold mass measured for the GAMA galaxies as a function of the average stellar mass of each stellar mass bin (Green crosses). We show all four estimates for for each sample a, slightly displaced along the yaxis for visibility. Horizontal errors correspond to 68% CI of ; vertical errors show the standard deviation of the stellar masses of galaxies within a sample. 

In the text 
Fig. B.1. Ratio of the standard deviation of the galaxy satellite number per halo and the Poisson variance. For Poissonian satellites this ratio is unity. Left: Ratio for all simulated galaxies (black crosses), red simulated galaxies (red crosses), and blue simulated galaxies (blue dots). Right: Ratio for stellar massselected samples. 

In the text 
Fig. B.2. Fractional difference of aperture statistics measurement ⟨𝒩^{a}𝒩^{b}ℳ⟩_{meas} (points) and modified aperture statistics model ⟨𝒩^{a}𝒩^{b}ℳ⟩_{mod}(lines) to the original model ⟨𝒩^{a}𝒩^{b}ℳ⟩ for redred lens pairs (left) and blueblue lens pairs (right) in the MS. The original model assumes a Poissonian satellite distribution, while the modified model uses the directly measured distribution of satellite galaxies. 

In the text 
Fig. C.1. G3L measurement in MS (points) and bestfitting halo model (lines) for stellar massselected lens samples, as defined in Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. Each row was fitted individually, leading to the χ^{2} values in the last panel. The corresponding halo model parameters are given in Table C.1. 

In the text 
Fig. C.2. Fig. C.1 continued 

In the text 
Fig. C.3. Fig. C.2 continued 

In the text 
Fig. C.4. G3L measurement in KV450 × GAMA (points) and bestfitting halo model (lines) for stellar massselected lens samples, as defined in Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1halo, dotted lines the 2halo, and dashdotted lines the 3halo term of the fit. Each row was fitted individually, leading to the χ^{2} values in the last panel. The corresponding halo model parameters are given in Table C.3 

In the text 
Fig. C.5. Fig. C.4 continued 

In the text 
Fig. C.6. Fig. C.5 continued 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.