KiDS+VIKING+GAMA: Halo occupation distributions and correlations of satellite numbers with a new halo model of the galaxy-matter bispectrum for galaxy-galaxy-galaxy lensing

Halo models and halo occupation distributions (HODs) are important tools to model the galaxy and matter distribution. We present and assess a new method for constraining the parameters of HODs using the gravitational lensing shear around galaxy pairs, galaxy-galaxy-galaxy-lensing (G3L). In contrast to galaxy-galaxy-lensing, G3L is sensitive to correlations between the per-halo numbers of galaxies from different populations. We use G3L to probe these correlations and test the default hypothesis that they are negligible. We derive a halo model for G3L and validate it with realistic mock data from the Millennium Simulation and a semi-analytic galaxy model. Then, we analyse public data from the Kilo-Degree Survey (KiDS), the VISTA Infrared Kilo-Degree 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_\odot$ and $10^{11.5}h^{-2} M_\odot$ and two colours (red and blue), as well as correlations between satellite numbers. The analysis recovers the true HODs in the simulated data within the $68\%$ credibility range. The inferred HODs vary significantly with colour and stellar mass. There is also strong evidence ($>3\sigma$) for correlations, increasing with halo mass, between the numbers of red and blue satellites and galaxies with stellar masses below $10^{10} \Msun. Possible causes of these correlations are the selection of similar galaxies in different samples, the survey flux limit, or physical mechanisms like 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 low-occupancy halos. The inferred HODs can be used to complement galaxy-galaxy-lensing or galaxy clustering HOD studies or as input to cosmological analyses and improved mock galaxy catalogues.


Introduction
Accurate models of the distribution of galaxies inside the cosmic large-scale 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 galaxy-galaxy-galaxy lensing (G3L) -the mean gravitational lensing shear around galaxy pairs . We demonstrate that with this higher-order statistic we can obtain the HODs of various galaxy samples A&A proofs: manuscript no. main sured galaxy and matter statistics, such as galaxy clustering (Zehavi et al. 2011;Ishikawa et al. 2021) and the galaxy-matter 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 second-order statistics of galaxies and matter. For example, halo model predictions for the galaxy-galaxy two-point correlation function agree well with observations (Zehavi et al. 2011). Mead et al. (2015) shows that halo models can also describe the non-linear matter power spectrum in N-body simulations with 5% accuracy.
Although halo models are prevalent for two-point statistics, it is unclear whether they are sufficiently accurate for modelling higher-order statistics, such as the galaxy-and matter bispectrum. These higher-order statistics contain complementary information to the two-point statistics (Berlind & Weinberg 2002) so it is worthwhile to expand halo models to include them to improve and cross-validate constraints from second-order statistics.
Here we extend halo models to measurements of G3L (Simon et al. 2008(Simon et al. , 2013Linke et al. 2020b), a third-order 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 galaxy-matter bispectrum integrated over the spread of lenses along the line-of-sight . 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 anti-correlated. While the default assumption in the literature is that satellite numbers are uncorrelated (Scranton 2001(Scranton , 2002Zehavi et al. 2005), some galaxy clustering 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 galaxy-galaxy-matter 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), 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 real-data application, we analyse G3L measurements in the overlap region of KiDS, VIKING, and GAMA, with lenses selected from GAMA and sources from KiDS+VIKING (Linke et al. 2020b). 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 cross-correlations between the per-halo 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, Fig. 1: Geometry of a G3L configuration with one source and two lens galaxies on the sky; adapted from . Lens galaxies are at angular positions θ 1 = θ + ϑ 1 and θ 2 = θ + ϑ 2 on the sky; the source galaxy is at θ. The angle between the source-lens connections is the opening angle φ.

Theory of galaxy-galaxy-galaxy 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 herafter). There are two types of this effect: The lensing of background galaxy pairs by matter around individual foreground galaxies (lens-shear-shear G3L) and the lensing of individual background galaxies by matter around foreground galaxy pairs (lens-lens-shear G3L). We concentrate on the latter. Figure 1 shows the geometric configuration for a lens-lensshear 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 source-lens connections. The main observables for G3L are the three-point correlation functionG ab and the closely related aperture statistics N a N b M , where a and b denote the considered galaxy populations. We describe the observables in Sects. 2.2 and 2.3. To modelG ab and N a N b M , we first introduce the galaxy-galaxy-matter bispectrum and its relation to the matter and galaxy number density distributions. and δ a g (x, z), defined as and δ a g (x, z) = n a g (x, z) − n a g (z) n a g (z) , 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 galaxygalaxy-matter bispectrum B ab ggδ (k 1 , k 2 ; k 3 , z) at comoving modes k 1 , k 2 , and k 3 and redshift z. This bispectrum is defined by δ a g (k 1 , z)δ b g (k 2 , z)δ(k 3 , z) = (2π) 3 δ D (k 1 + k 2 + k 3 ) B ab ggδ (k 1 , k 2 ; k 3 , z) 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 B ab ggδ (k 1 , k 2 ; −k 1 − k 2 , z) =: B ab ggδ (k 1 , k 2 , z) .
As a gravitational lensing effect, G3L depends mainly on the projections κ and κ a g of δ and δ a g along the line-of-sight. 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 N a g as N a g (θ) = ∞ 0 dχ p a (χ) n a g χ(θ, χ), z(χ) .
As shown in Kaiser (1992), under the assumptions of the Limber equation, the projected bispectrum is b ab ggκ ( 1 , 2 ) = 3H 2 0 Ω m 2c 2 dχ q(χ) p a (χ) p b (χ) χ 3 a(χ) × B ab ggδ 1 , 2 , z(χ) , which we use to model the observables of G3L: the correlation functionG ab and the aperture statistics N a N b M . 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 Mpc) of the galaxy distribution (Bartelmann & Schneider 2001).

G3L correlation function
The G3L correlation functionG ab correlates the projected lens galaxy number densities N a g and N b g 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,G ab only depends on ϑ 1 = |ϑ 1 | and ϑ 2 = |ϑ 2 | of the lens-source 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,G ab can be calculated from the projected galaxy-galaxy-matter bispectrum b ab ggκ (c.f. their Equation 40), by integrating b ab ggκ ( 1 , 2 ), multiplied with a complex kernel function containing the second-order Bessel function. Consequently, G ab 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 ofG. We describe them in the following subsection.

G3L aperture statistics
Aperture statistics are moments of aperture masses M θ (ϑ) and aperture number counts N p θ (ϑ), and with a filter function U θ of aperture scale radius θ. The filter function U θ has to be compensated, i.e., dϑ ϑ U θ (ϑ) = 0. For G3L, the relevant aperture statistics are N a N b M , given by Article number, page 3 of 35 A&A proofs: manuscript no. main The aperture statistics are related to the bispectrum b ab ggκ 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 higher-order aperture statistics (e.g. , 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 functionG ab can be connected analytically with the kernel function A NNM (ϑ 1 , ϑ 2 , φ | θ 1 , θ 2 , θ 3 ) given in the appendix of SW05. As alluded to before, there are practical advantages to modelling N a N b M instead ofG ab . First, evaluating Eq. (17) is numerically more stable than calculatingG from the bispectrum because the filter function U θ is localized and non-oscillating. Second, in contrast toG, the aperture statistics do not depend on the galaxy-matter two-point correlation function. Therefore, they do not need a model of the galaxy-galaxy-lensing signal γ t . Third, a N a N b M data vector is a condensed summary statistic, where a few aperture radii (∼ tens) contain a similar amount of information asG over hundreds of bins.
Consequently, we use the aperture statistics N a N b M as the primary observable and model it with Eq. (17), based on a halo-model based bispectrum b ab ggκ . To measure the aperture statistics, we estimateG ab (see Sect. 5.1) and convert it to aperture statistics with Eq. (19). We focus on the equilateral statistics, i.e., θ 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.

G3L halo model
In this section, we derive the galaxy-galaxy-matter 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 second-order moments of the HOD. We will give details on these ingredients next.

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 ρ vir = ∆(z)ρ, wherē ρ 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) Navarro-Frenk-White 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 normalized 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 , i.e., σ(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Ŵ(x), The third group of ingredients are the HMF n(m, z) dm and the halo bias. The HMF describes the comoving number density of dark matter halos with mass between m and m + dm. We use the HMF by Sheth & Tormen (1999), where the parameters A = 0.322, p = 0.3, and q = 0.707 were found in N-body simulations, and ν = δ c /σ(r m , z) for r 3 m = 3m (4πρ) −1 . 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 so-called peak-background split formalism (Mo & White 1996;Scoccimarro et al. 2001) the linear bias b 1 is with the q and p as in the Sheth-Tormen 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), Fourth, for the average number density of satellite galaxies a inside halos of mass m, N a sat u a g (x | 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, e.g., introduced by Cacciato et al. (2012) in galaxy clustering 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, N a cen = 1, situated at the halo centre but can contain several satellite galaxies, N a sat . With this split, the galaxy-galaxy-matter bispectrum requires the specification of the mean numbers N a cen | m and N a sat | m of central and satellite galaxies, as well as the numbers of central and satellite pairs N a cen N b cen | m , N a sat N b sat | m , for a b, and N a sat (N a sat − 1) | m for a = b. The mean number N a cen | m of central galaxies depends only on halo mass. For small halo masses, no galaxy formation occurs, so N a cen | m = 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 , M a th , 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 N a cen | m from 0 to α a . If σ a is small, the transition from N a cen | m = 0 to N a cen | m = α 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 M a th ) 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 N a sat | 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 low-mass halos. Satellite galaxies predominate in massive halos. The galaxy-galaxy matter bispectrum depends not only on the expected numbers of centrals and satellites, but also on the expected number N a sat (N a sat − 1) | m of unmixed satellite pairs, the expected number N a sat N b sat | m of mixed satellite pairs, i.e, a b, and the number N a cen N b sat | m of central-satellite pairs per halo of mass m. In principle, it also depends on the number N a cen N b cen | m of central pairs, but this number is 0 for all m, as each halo contains only one central galaxy. We assume, as, e.g., 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.  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.
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 anti-correlated; and positive if they are positively correlated. If the satellite galaxies are distributed Poissonian, i.e., σ 2 (N a sat | m) = N a sat | 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 cross-correlation between galaxy samples is independent of halo mass if = 0. Ignoring a mass dependence of r ab in a more simplistic approach, as, e.g., 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 halo-mass dependent r ab (m), on the other hand, is easier to interpret, which is why we use the mass-dependent 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 N b sat = p b + p ab for galaxies b. Similarly, the mean number of pairs is N a sat N b sat = p ab , and the variances are σ 2 (N a sat ) = (N a sat ) 2 − N a sat 2 = (p a + p ab )(1 − p a − p ab ) and σ 2 (N b sat ) = (p b + p ab )(1− p b − p ab ). The Pearson cross-correlation coefficient of variations in the galaxy numbers consequently becomes Therefore, at the extreme end, for halos too small to host more than one galaxy, i.e., p ab = 0, the value of r ab pear exactly converges to zero for p a , p b → 0 (that is N a sat , N b sat → 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 r ab pear 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 r ab pear for m → 0 is a natural outcome of halos sparsely populated with galaxies a and b. This also applies to r ab for sub-Poisson variances, expected at low halo masses, since r ab < r ab pear 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 N|m 1, or smaller for N|m 1. Consequently, r ab is at small halo masses systematically lower than the Pearson correlation coefficient (r ab < r ab pear ) and higher for high halo masses (r ab > r ab pear ). 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 galaxy-galaxy-matter bispectrum B ggδ . We derive the bispectrum in the following subsection.

Galaxy-galaxy-matter bispectrum
Using Eqs. (1), (2) and (3), the galaxy-galaxy-matter 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 1-halo term 1 B ab ggδ (k 1 , k 2 , z), the 2-halo term 2 B ab ggδ (k 1 , k 2 , z), and the 3-halo term 3 B ab ggδ (k 1 , k 2 , z), or together The 1-halo term depends on the correlation between galaxies numbers and matter density in the same halo. One part of the 2-halo 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 3-halo term.
To deriveρ andn, 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 normalized density profile u the matter density field is whose Fourier transform iŝ 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 ∆x a i j from the halo centre. Centrals are exactly at the halo centre. The number density of galaxies from sample a is therefore where N a sat,i is the number of satellite galaxies and N a cen,i the number of central galaxies from sample a inside halo i. The Fourier transform of the number density iŝ To derive the bispectrum we insert Eqs. (49) and (51) into Eq. (46), which leads to 1 n a g (z)n b g (z)ρ(z) dm n(m) × mû(k 1 + k 2 | m, z) G ab (k 1 , k 2 | m, z) ; 2 B ab ggδ (k 1 , k 2 , z) = 1 n a g (z)n b g (z)ρ(z) dm 1 dm 2 n(m 1 ) n(m 2 ) and where the Kronecker symbol δ K ab 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.

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, i.e., 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 NNM (θ, p).
The priors are aimed to be uninformative. For σ, log 10 (M a th /M ), β a , and log 10 (M a /M ), they are based on the galaxy-galaxy-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 sanity-check 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 N a N a M for unmixed lens pairs, when varying α a , σ a , M a th , 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 M a th and the satellite mass scale M a are tightly correlated: Increasing M a th while decreasing M a could lead to the same N a N a M .
The threshold halo mass M a th and the slope β a of the satellite HOD affect the signal the strongest on the scales we consider here. Increasing the threshold mass M a th from 3.2 × 10 12 M to 7.5 × 10 12 M roughly doubles the aperture statistics across the whole range of θ from 0. 1 to 100 because M a th is closely connected to the galaxy bias. A higher M a th 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 th 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 N q sat |m has the strongest impact on the statistics at θ ∼ 9 , an up to a 20-fold increase in amplitude. While a larger β a increases the signal on larger scales, at scales below 0. 3, 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 low-mass halos with fewer satellite pairs. This increases N a N a M . 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 lower-mass halos are more relevant.
The concentration parameter f a affects the aperture statistics stronger at small scales than at large scales (57% at 0. 1 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 galaxy-galaxy-matter bispectrum and the aperture statistics increase at small scales. At scales above 11. 5, though, the change is less than 10%.
The sensitivity of the statistics for r ab is visible in Fig. 4, where we show N a N b M for mixed lens pairs for correlated satellite numbers (r ab = 1), uncorrelated numbers (r ab = 0), or anti-correlated 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 N a N b M changes strongest for θ = 1. 1 (72%), where the signal is dominated by the 1-halo term containing N a sat N b sat |m and has a r ab -dependence. The cross-correlation parameter has only a small effect (less than 10% change) at scales above 30 because these scales are dominated by the 3-halo term, independent of the per-halo number of satellite pairs.
In conclusion, the aperture statistics between 0. 1 and 100 are most sensitive to the threshold halo mass M a th and the slope β a . In contrast, we do not expect G3L to improve the σ a constraints beyond the prior range. To probe the cross-correlation parameter r ab , the statistic N a N b M needs to be measured at scales below 30 and preferably around 1 . Both of these measurements are achieved for our survey data.

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. (2020b): 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: Aghanim et al. (2020). This difference in cosmology should not impact our conclusions because we are not interested in a oneto-one comparison between the observation and simulation and merely use the mock data to validate our method.

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 near-infrared 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 r-band. 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.

1-Halo Term 2-Halo Term 3-Halo Term Total
Notes. Lenses are selected either according to their stellar mass M * or to their rest-frame (g − r) 0 colour and absolute r-band magnitude M r and need to have r < 19.8 mag.
Galaxies observed by GAMA are brighter than r = 19.8 mag, rendering our lens galaxy sample flux-limited. 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 r-band, and five samples defined by stellar masses. For this we use the same cuts as in Farrow et al. (2015) and Linke et al. (2020b). An overview of these cuts is given in Table 2.
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 stellar-mass 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.
The differences in the redshift distribution of galaxies from different samples are also visible in Fig. 6, which shows the stel- lar 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.

Millennium Simulation
The simulation data are constructed from the MS (Springel et al. 2005), a dark-matter only cosmological N-body 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 multiplelens-plane 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 lines-of-sights (total area of 1024 deg 2 ) at each redshift slice of the MS. 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, i.e., 291 arcmin −2 .
Since all 64 lines-of-sight originate from the same simulation, they are not independent at the largest scales. However, the lines-of-sight correspond to different observers placed in the simulation box such that the overlaps between the lines-of-sights 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 1-halo term, originating from small scales of approximately 1 h −1 Mpc. Therefore, we treat the lines-of-sight 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. (2020b), 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).

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. (2020a) 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. For an assessment of 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 2 .

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 estimatorG ab est 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 ofG ab (B) is the real part of where ω ab (θ) is the angular two-point correlation function of lens galaxies from samples a and b with lag θ, and The angles ϕ ik and ϕ jk are the polar angles of the lens-source separation vectors θ i − θ k and θ j − θ k (corresponding to ϕ 1 and ϕ 2 in Fig. 1) and φ i jk = ϕ ik − ϕ jk is the opening angle between the lens-source 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 ofG ab est (B) is pure noise in the absence of systematic errors in the shear data (Linke et al. 2020b). We give in the estimator equal weight to all lens pairs, which is in contrast to the recent analysis in Linke et al. (2020b), where we weighted galaxy pairs based on the redshift difference between their members. The lack of weighting dilutes the signalto-noise ratio of the G3L signal since non-physical pairs (widely separated along the line-of-sight) 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 line-of-sight 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 down-weighting 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 two-point correlation function of the lens clustering, ω ab (θ), is that by Szapudi & Szalay (1998), for two lens samples a and b with N a d and N b d galaxies, and two 'random samples'. These random samples contain N a r and N b r 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)G ab 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 theG ab est with jackknife resampling on a tile by tile basis. We estimateG ab i 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 ϑ 1 , ϑ 2 ∈ [0. 15, 85 ] and φ ∈ [0, 2π]. For the totalG ab est , individual tile estimatesG ab i 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 lens-lens-source 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. (2020a). 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. (2020a) that estimatingG in 128 × 128 × 128 bins and then applying the tesselation leads to a measurement accuracy within 5%.
We convert a binnedG ab est 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 A NNM (b i | θ 1 , θ 2 , θ 3 ) is the kernel function of Eq. (19), evaluated at the center of bin b i . We estimate the covariance matrix of the estimates with jackknife resampling. For this, we combine theG ab estimates for all but the k-th tile to the k-th jackknife sampleG ab k,jn , which is then converted to the aperture statistics, Eq. (19), leading to N sam- 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 ofG are statistically independent. However, as all observed tiles originate from the same observation and are adjacent to each other, this assumption is not exactly valid. It is still probably a good approximation on scales smaller than the tile sizes.
Article number, page 11 of 35 A&A proofs: manuscript no. main A possible bias in the empirical covariance is less pronounced for the simulation, as the tiles originate from 64 (mostly) independent line-of-sights.

Fitting the halo model
We constrain the parameters of the G3L halo model by fitting it to measurements of both the auto-correlation aperture statistics, N a N a M (θ) and N b N b M (θ), and cross-correlation statistics, N a N b M (θ), of two galaxy samples a and b for N θ = 30 scale radii θ between 0. 1 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 Eqs. (46) and the bispectrum in Sect. 3.2, Optimal parameters p opt are those that minimise the goodness-of-fit determined by the Nelder-Mead algorithm (Nelder & Mead 1965) as implemented in the GNU Scientific Library (Gough 2009). This algorithm is well-suited to multi-dimensional 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 best-fitting 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 L(d ab | p) is the likelihood of the data given the parameters p. We assume Gaussian statistical errors in the data, i.e., the likelihood function is with the χ 2 as defined in Eq. (66); the Bayesian evidence and thus the normalisation of L 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 L(d ab | p) in the proximity of the optimal parameters p opt by the Gaussian probability densitỹ 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 normalized such that their sum is unity. For each parameter p i we give the α credibility interval (CI) on the marginalized posterior P(p i |d ab ), defined by We take the mode of the posterior as optimal parameter value p i opt and define the CI I α of a single parameter as the interval around the optimal parameter value including α of the marginalized posterior, i.e., α = I α dp i P(p i |d ab ) .
To find this interval, we sort the sampling points by the distance |p i − p i opt | 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.

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 friends-of-friends 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 = N a cen | m + N a sat | 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 lines-of-sights of the simulation.

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.

Colour-selected 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 3-halo terms. The parameter values corresponding to the best fit are given in the first two columns in Table 3. The goodness-of-fit is χ 2 = 93.63 for 90 − 14 = 76 degrees-of-freedom (dof), or χ 2 /dof = 1.28. This χ 2 corresponds to a p-value of 0.083, indicating no significant deviation between the fit and the simulated G3L signal within the 95% confidence level (CL). For our science verification, we compare in Fig. 8a the HODs inferred by the best-fitting G3L halo model to the directly estimated HODs of the simulated galaxies. 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 red-blue and red-blue , 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.
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 3-halo terms. The model fit has χ 2 /dof = 0.977, corresponding to a p-value 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: M red th = 1.5 +6.7 −1.2 × 10 12 M is roughly ten times larger than M blue th = 1.4 +7.2 −0.9 × 10 11 M (68% credibility intervals, CI herafter). The per-halo number of blue satellites increases slower with halo mass than for red galaxies: the mass scale M blue = 2.0 +3.0 −1.2 × 10 14 M of blue satellites is more than five times larger than M red = 3.6 +3.5 −0.7 × 10 13 M (68% CI). For central galaxies, the sum of α red and α blue is 0.47 +0.51 −0.31 (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 non-zero for halo masses above 10 12 M (5 × 10 11 M ), whereas, at lower halo masses, the constraints become essentially upper limits for N|m . 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 N a N b M 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 red-blue and red-blue 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, ab = 0.99 +0.22 −0.12 versus ab = 0.69 +0.08 −0.04 at 68% CI), while the correlation amplitude at 10 12 M is lower (A ab = 1.62 +0.62 −0.51 × 10 −2 versus A ab = 5.31 +0.87 −0.92 × 10 −2 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 r red−blue 0.16 +0.06 −0.05 for red−blue = 1 (68% CI). 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 red-red and blue-blue 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 red-red pairs.

Stellar mass-selected lens samples
We repeat the science verification and G3L analysis for the stellar mass-selected galaxies. We consider five stellar-mass bins m1 to m5, so there are ten distinct combinations of two samples a and b. For each combination, we again use the statistics N a N a M , N a N b M , and N b N b M 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 Article number, page 13 of 35 A&A proofs: manuscript no. main  Fig. 7: G3L measurement (points) and best-fitting 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 1-halo, dotted lines the 2-halo, and dash-dotted lines the 3-halo term of the fit. The left panels show the result for red-red galaxy pairs, the central panels for blue-blue galaxy pairs, and the right panels for red-blue mixed pairs. realistic model and successful fits, these four HODs should be consistent with each other. Fitting the halo model to two samples out of five individually 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 stellar-mass 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 p-values 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 realizations of the model. A Kolmogorov-Smirnov (KS) test on the distribution of p-values for the simulation yields a KS distance of 0.118, which for the 11 samples and dof = 76 in the distribution indicates no significant deviation from a uniform distribution at 95% CL. Additionally, the four parameter sets for each stellar-mass bin agree within the 68% CI. The distribution of p-values 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.
To verify the reconstruction of the HODs for stellar-mass samples, we compare the inferred HODs (lines and shaded areas) to the true HODs (data points) in Figure 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 N|m 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, N|m 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.
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 Fig. 8: Mean per-halo 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 line-of-sights. The lines indicate the per-halo 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.
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 1.20 +0.45 −0.36 × 10 −2 , whereas A m1m5 is 0.029 +0.024 −0.020 × 10 −2 (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).
After the science verification, we fit the model to the observed G3L signal of KV450×GAMA. The best-fitting model parameters, χ 2 -values and plots are reported in Table C.3. For the observations, all p-values exceed 0.05, indicating no significant deviation between fit and measurement at a 95% CL. Again, we perform a KS-test on the cumulative distribution of p-values 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  For a perfect description of G3L signal and data noise, the distributions would be consistent with a uniform distribution (black, dotted). 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 stellar-mass 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 M a th 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 M a th against the average stellar mass of the samples. For each stellar mass bin, there are four estimates of M a th , from the four possible sample combinations, slightly displaced along the y-axis. As noted above, these four estimates are consistent with each other for all samples.
The central galaxies parameter α a also increases with stellar mass: The central galaxy of a halo with m max{M a th , M b th } 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ᾱ a for each sample a, yieldingᾱ m1 = 0.11 +0.07 −0.04 ,ᾱ m2 = 0.14 +0.08 −0.06 ,ᾱ m3 = 0.22 +0.10 −0.08 , α m4 = 0.47 +0.11 −0.13 ,ᾱ m5 = 0.54 +0.19 −0.19 . The sum of these α a is 1.48 +0.26 −0.25 , 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 aᾱ a requires a simultaneous fit of all 10 combinations of the five mass bins. This is, unfortunately, not feasible in our case because our jackknife-derived 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.
We remind here that r ab > 1 are possible because the r ab are a Pearson correlation coefficient, r ab pear hereafter, only for a Poissonian variance of satellite numbers. In particular, we have the relation so that values of r ab > 1 indicate a super-Poissonian variance. For the strongly non-Poissonian m1 and m2 SAM galaxies at m = 10 14 M , r m1 m2 = 1.2, which corresponds to r m1 m2 pear = 0.18 (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 r m4 m5 pear = 0.04. For red and blue SAM galaxies, the measured r red blue = 3.2 at m = 10 14 M corresponds to r red blue pear = 0.44. 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 m2-galaxies. 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 the 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, 10 12 and observation (bottom) for lens galaxies from each stellar mass bin as function of halo mass. Crosses indicate the directly estimated per-halo numbers of simulated galaxies, while lines show the predictions from the G3L fits. The shaded areas indicate the 68% confidence areas. The left panels show the mean per-halo 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. irrespective of their magnitude, and divide them by the same stellar-mass cuts as the flux-limited 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 flux-limited 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.

Discussion
We presented a new method to measure galaxy HODs and the correlation of per-halo satellite galaxy numbers using G3L, a weak gravitational lensing effect measuring the projected galaxy-galaxy-matter bispectrum. To this end, we constructed a new halo model for the galaxy-galaxy-matter 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 Halo mass [M ]  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 best-constrained HOD parameter for the KV450 × GAMA is the threshold halo mass M a th , 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 ex-pectation 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 N|m , 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.

Mass bins
Fig. 14: Threshold mass M a th 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 M a th for each sample a, slightly displaced along the y-axis for visibility. Horizontal errors correspond to 68% CI of M a th ; vertical errors show the standard deviation of the stellar masses of galaxies within a sample.
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 N a cen |m by a step function.
The 1-halo term of the G3L signal of red-red lens pairs stretches to larger scales than for blue-blue 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 1-halo term extends to larger scales for red-red lens pairs than for blue-blue ones. However, mixed red-blue pairs exist in intermediate halos, large enough to host red galaxies but small enough to contain a significant fraction of blue galaxies. Consequently, the cross-over between the domination of the 1-halo and the 3-halo term occurs at larger scales than for blue-blue 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 (r 0.16 +0.06 −0.05 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 (e.g., 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 low-mass halos, m 10 13 M , are considered.
The obvious presence of correlations, especially for halos of the group-and cluster-mass scale, questions the assumption of no correlation in halo models for the cross-correlation statistics between galaxy samples in galaxy clustering studies (e.g., Scranton 2001Scranton , 2002Zehavi 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 cluster-sized halos (m 10 14 M ): Galaxy models commonly assume that star-forming 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.
Using G3L to infer the correlations between the HODs of different galaxy populations could be interesting for other sample selections, for example, Luminous Red Galaxies (LRGs) or Emission Line Galaxies (ELGs). It would also be interesting to investigate the evolution of the cross-correlation 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 sub-halos, known to be relevant for galaxy-galaxy lensing at sub-arcmin scales (e.g., Velander et al. 2014). 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 N a cen |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 N a cen |m above a certain halo mass. 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 sub-Poissonian variance for low-mass halos and super-Poissonian variance for high-mass 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 N 2 sat |m 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 non-Poissonian 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, better HOD constraints could be achieved by incorporating the redshift-weighting scheme of lens pairs by Linke et al. (2020a). In this scheme, lens galaxy pairs that appear close in projection but are physically distant are down-weighted compared to physically close lens pairs in estimating the G3L correlation function. This weighting increases the signal-to-noise 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 line-of-sight 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 in the future.
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., Hildebrandt et al. 2021) or inference pipelines (e.g., Ferrero et al. 2021;DeRose et al. 2021). Including the cross-correlation 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 stellar-to-halo mass ratio, traditionally obtained from second-order statistics, such as galaxy-galaxy-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.

Appendix A: Calculation of galaxy-galaxy-matter bispectrum
Here we derive the galaxy-galaxy-matter 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 , 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 depend only 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 sub-halos and that galaxy positions inside halos are statistically independent from each other. These approximations greatly simplify the model formalism, while still being suffienciently 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 P a N (N a cen,h , N a sat,h | m h ) is the joint probability for halo h to have N a cen,h central and N a sat,h satellite galaxies from sample a, and wheren h = H/V is the mean halo number density 3 . The PDF P ab g is a product of the normalized spatial distributions u a g (∆x|m) of satellites in each halo, under the assumption of independent satellite positions inside a halo. The joint probability P N of per-halo galaxy numbers has the moments where p, q, r, s ∈ N 0 . After inserting Eqs. (49) and (51) into (46) and using Eq. (A.1), the bispectrum B ab ggδ for mixed pairs, a b, is given by where unconnected terms do not contain δ D (k 1 + k 2 + k 2 ),n a g 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 δ K i j is the Kronecker symbol. We divide the sum over i, j and k into the 1-halo term 1 B ab ggδ for i = j = k, the 2-halo term 2 B ab ggδ for i = j k, i = k j, and j = k i, and the 3-halo term 3 B ab ggδ for i j k. In the following, we calculate these terms individually.

Appendix A.1: 1-halo term
To calculate the 1-halo term 1 B ab ggδ with i = j = k, we distinguish between mixed (a b) and unmixed (a = b) pairs. For a = b, the 1-halo term is given by Evaluating all m-integrals aside from the one over m i (H − 1 integrals in total) leads to Using Eqs. (A.6), d 3 x u a g (x | m) = 1, and d 3 x u a g (x | m) exp(−ik · x) =û a g (k | m), we find explicitly using the second-order moments of P N in the last line. We consider large volumes V here, rendering the V-integrals 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 Dirac-distribution (e.g., Weisstein 2022), as (A.14) In this limit (assumed henceforth), For mixed pairs, a b, an analogous calculation gives In the main text, we denote 1 B ab ggδ (k 1 , k 2 , k 3 ) as 1 B ab ggδ (k 1 , k 2 ) for brevity.

Appendix A.2: 2-halo term
The 2-halo term 2 B ab ggδ in Eq. (A.7) contains the contributions from i = j k, i = k j, and k = j i. For unmixed pairs, a = b, the 2-halo term is given by × m jû (k 3 | m j ) e −ik 3 ·x j N a cen,i (N a cen,i − 1) e −i(k 1 +k 2 )·x i + N a .
In the main text we just denote 2 B ab ggδ (k 1 , k 2 ).
We evaluate all m-integrals 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 3 B aa ggδ (k 1 , k 2 , k 3 ) + unconnected terms Now, again in the limit of an infinite volume V, Eq. (A.14), we use the two-and three-point 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 Again, we just denote 3 B ab ggδ (k 1 , k 2 ) in the main text. Left: Ratio for all simulated galaxies (black crosses), red simulated galaxies (red crosses), and blue simulated galaxies (blue dots). Right: Ratio for stellar mass-selected samples.

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 non-Poisson satellites. To some degree, deviations from Poisson satellites are indeed visible in our SAM galaxies: Figure B.1 shows the ratio of σ(N a sat |m) to the Poisson variance N a sat |m for different samples of the simulated galaxies. Poisson satellites, where this ratio is unity, are present in the intermediate halo-mass range, such as between 2 × 10 11 M m 10 12 M for galaxies from stellar-mass bin m2. On the low-mass end, satellites are distributed sub-Poissonian (ratio is below unity), and on the high-mass end they are super-Poissonian (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.
We estimate the bias induced by the mismatch between the true satellite variance and the model assumption by computing the fractional change of N a N b M 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 N a N b M mod , using the HOD parameters from the best-fit but updating the number of satellite pairs inside a halo to where σ 2 (N a sat |m) is the true satellite variance in the simulation. Clearly, for σ(N a sat |m) = N a sat |m , i.e., 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 strongest 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 1-halo term, containing N a sat (N a sat − 1)|m , dominates here. Smaller scales are dominated by the halo term containing N a cen N a sat |m , whereas larger scales are dominated by the 3-halo 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 non-Poisson satellites affects the model prediction, a Poisson model is accurate enough for the G3L analysis of measurements in this work.
We have repeated this test also for the satellite galaxies from the stellar mass-selected 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, i.e., r ab = 1 for perfectly correlated N a sat and N b sat , r ab = −1 for perfectly anti-correlated satellite numbers, and r ab = 0 for uncorrelated samples. However, for a super-Poisson variance in the high-mass regime, our r ab is larger than the Pearson coefficient. Conversely, for low-mass halos with sub-Poisson variance, r ab is smaller than the Pearson coefficient. This in combination 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 mass-selected galaxies
In Fig C.1 we show the model fits to the G3L aperture statistics from stellar mass-selected lenses in the MS and the χ 2 of the fits. All fits have 76 degrees-of-freedom (the same as the fit to the colour-selected samples in Sect. 6.1), so a χ 2 < 97.4 signifies agreement between fit and data at the 95% CL. The HOD parameters of the best-fitting models are given in Table C.1; the parameters A ab and ab that determine the cross-correlation between satellite numbers together with the χ 2 and p-values of the fits are given in Table C.2. Figure C.2 shows the model fits to the G3L aperture statistics from stellar mass-selected lenses in KV450× GAMA and the χ 2 of the fits. The parameters of the best-fitting models are given in Tables C.3 and C.4.   Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1-halo, dotted lines the 2-halo, and dash-dotted lines the 3-halo 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 Table 2. Solid lines indicate the total aperture statistics, dashed lines the 1-halo, dotted lines the 2-halo, and dash-dotted lines the 3-halo 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 Table C.3: Best-fit values for halo model parameters for stellar-mass-selected lenses in KV450 × GAMA for each stellar mass sample a.

From correlation with m1
From correlation with m2 From correlation with m3 From correlation with m4 From correlation with m5