© L. Linke et al. 2022

1. 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 (Schneider & Watts 2005). We demonstrate that with this higher-order 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 per-halo numbers of satellite galaxies from different populations. In a first application, we infer the HODs and correlations of various galaxy selections in the Kilo-Degree Survey (KiDS; Kuijken et al. 2015), the VISTA Infrared Kilo-degree 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 large-scale 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 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 semi-analytic 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, 2013; Linke et al. 2020a), 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 (Schneider & Watts 2005). Unlike galaxy-galaxy 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, 2002; Zehavi 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) 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 real-data 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⁻² 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, H₀ = 73 km s⁻¹ Mpc⁻¹, and σ₈ = 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₀ = 67.4 km s⁻¹ Mpc⁻¹, and σ₈ = 0.811. Throughout we assume a flat cosmology with Ω_Λ = 1 − Ω_m.

2. 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 hereafter). 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-lens-shear 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 θ + ϑ₁ and θ + ϑ₂. 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 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 galaxy-galaxy-matter bispectrum and its relation to the matter and galaxy number density distributions.

2.1. Projected fields and the galaxy-galaxy-matter 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 late-type galaxies. Fluctuations in the densities ρ and are the matter- and galaxy number density contrast, δ(x, z) and , defined as

(1)

and

(2)

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 galaxy-galaxy-matter bispectrum at comoving modes k₁, k₂, and k₃ and redshift z. This bispectrum is defined by

(3)

where hats denote Fourier transforms for which we use the convention

(4)

The bispectrum is defined only for k₃ = −k₁ − k₂, so we abbreviate

(5)

As a gravitational lensing effect, G3L depends mainly on the projections κ and of δ and along the line-of-sight. In a flat universe the lensing convergence is

(6)

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

(7)

The projected galaxy number density contrast is given by

(8)

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

(9)

To arrive at a G3L signal, we convert the galaxy-galaxy-matter bispectrum to its projected counterpart , defined by

(10)

where hats again denote Fourier transforms. Similar as for , we abbreviate

(11)

As shown in Kaiser (1992), under the assumptions of the Limber equation, the projected bispectrum is

(12)

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,

(13)

Due to the statistical homogeneity and isotropy of the matter and galaxy fields, only depends on ϑ₁ = |ϑ₁| and ϑ₂ = |ϑ₂| of the lens-source separations, and the angle ϕ between ϑ₁ and ϑ₂. 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 galaxy-galaxy-matter bispectrum (cf. their Eq. (40)), by integrating (ℓ₁, ℓ₂), multiplied with a complex kernel function containing the second-order 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 ,

(14)

and

(15)

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

(16)

The aperture statistics are related to the bispectrum by

(17)

where is the Fourier transform of U_θ. We choose the exponential filter function by Crittenden et al. (2002),

(18)

which is commonly used for studies of higher-order 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

(19)

with the kernel function 𝒜_𝒩𝒩ℳ(ϑ₁, ϑ₂, ϕ ∣ θ₁, θ₂, θ₃) 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 non-oscillating. Second, in contrast to , 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 ⟨𝒩^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 halo-model 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, θ₁ = θ₂ = θ₃, to reduce the size of the data vector and use the shorthand

(20)

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 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 give details on these ingredients next.

3.1. Ingredients

Halo models assume that all matter and galaxies are within self-bound, 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,

(21)

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₀| | m, z) of a halo with mass m and redshift z, centred at r₀. We assume halos follow (truncated) Navarro-Frenk-White profiles (NFW; Navarro et al. 1996), given as

(22)

(23)

where r₂₀₀ is the radius of a sphere around the halo centre, in which the mean density¹ 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),

(24)

with c₀ = 9 and α = 0.13. The mass m_⋆ is that enclosed by a sphere of radius r_⋆,

(25)

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 ,

(26)

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),

(27)

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 .

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,

(28)

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₁ is

(29)

with the q and p as in the Sheth-Tormen HMF. Assuming a deterministic halo bias, the halo power spectrum is then

(30)

and the halo bispectrum is

(31)

for the linear matter bispectrum by Bernardeau et al. (2002),

(32)

where k₃ = −k₁ − k₂ and

(33)

Fourth, for the average number density of satellite galaxies a inside halos of mass m, , we assume an NFW profile with concentration

(34)

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 galaxy-clustering halo models, who show that f^a affects galaxy clustering at scales below 1 h⁻¹ 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

(35)

Each halo hosts at most one central galaxy, , situated at the halo centre but can contain several satellite galaxies, . With this split, the galaxy-galaxy-matter 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

(36)

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),

(37)

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 ⟨ ( − 1) | m⟩ of unmixed satellite pairs, the expected number ⟨  | m⟩ of mixed satellite pairs, i.e, a ≠ b, and the number ⟨ | m⟩ of central-satellite 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,

(38)

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,

(39)

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 cross-correlation coefficient r^ab of satellites, defined for mixed galaxy samples a ≠ b by

(40)

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, that is, σ²( | 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

(41)

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, 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 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

(42)

(43)

(44)

and ⟨⟩=p_b + p_ab for galaxies b. Similarly, the mean number of pairs is ⟨ ⟩=p_ab, and the variances are σ²() = ⟨()²⟩−⟨⟩² = (p_a + p_ab)(1 − p_a − p_ab) and σ²() = (p_b + p_ab)(1 − p_b − p_ab). The Pearson cross-correlation coefficient of variations in the galaxy numbers consequently becomes

(45)

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 sub-Poisson 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 ⟨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 < ) 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 galaxy-galaxy-matter bispectrum B_ggδ. We derive the bispectrum in the following subsection.

3.2. Galaxy-galaxy-matter bispectrum

Using Eqs. (1)–(3), the galaxy-galaxy-matter bispectrum is given by

(46)

where a and b denote the galaxy samples, and ‘unconnected terms’ are those proportional to δ_D(k₁), δ_D(k₂), or δ_D(k₃), which do not affect the bispectrum. The bispectrum can be divided into three terms: the 1-halo term ₁(k₁, k₂, z), the 2-halo term ₂(k₁, k₂, z), and the 3-halo term ₃(k₁, k₂, z), or together

(47)

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 and , we use that in the halo model the cosmic density field consists of H halos with masses {m₁, …, m_H} at positions {x₁, …, x_H}. With the normalised density profile u the matter density field is

(48)

whose Fourier transform is

(49)

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

(50)

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

(51)

To derive the bispectrum we insert Eqs. (49) and (51) into Eq. (46), which leads to

(52)

(53)

and

(54)

with

(55)

and

(56)

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,

(57)

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₁₀(/M_⊙), β^a, and log₁₀(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 ⟨𝒩^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ℳ⟩.

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¹² M_⊙ to 7.5 × 10¹² 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 20-fold 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 low-mass 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 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 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 , 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 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 ⟨𝒩^a𝒩^bℳ⟩ changes strongest for (72%), where the signal is dominated by the 1-halo term containing ⟨|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 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 cross-correlation 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 one-to-one 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². 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² 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. 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 flux-limited. We divide the GAMA galaxies into different samples: a ‘red’ and a ‘blue’ sample, using their restframe (g − r)₀ 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. (2020a). 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 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.

4.2. 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⁻¹Mpc and contains 2160³ particles with mass 8.6 × 10⁸ h⁻¹ M_⊙ (h = 0.73).

Maps of the lensing shear, γ, were created with the multiple-lens-plane raytracing algorithm by Hilbert et al. (2009). This algorithm generated shear maps of size 4 × 4 deg² on a regular grid of 4096 × 4096 deg² for 64 lines-of-sights (total area of 1024 deg²) 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, that is, 291 arcmin⁻².

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⁻¹ Mpc. In contrast, the G3L signal is dominated by the 1-halo term, originating from small scales of approximately 1 h⁻¹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. (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³.

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 (ϑ₁, ϑ₂, ϕ) 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

(58)

where ω^ab(θ) is the angular two-point correlation function of lens galaxies from samples a and b with lag θ, and

(59)

The angles φ_ik and φ_jk are the polar angles of the lens-source separation vectors θ_i − θ_k and θ_j − θ_k (corresponding to φ₁ and φ₂ in Fig. 1) and ϕ_ijk = φ_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 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 signal-to-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),

(60)

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_aD_b, D_aR_b, D_bR_a, and R_aR_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 ϑ₁ and ϑ₂, 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 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. (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),

(61)

where N_bin is the number of bins, V(b_i) is the tesselation volume of bin b_i, and 𝒜_𝒩𝒩ℳ(b_i | θ₁, θ₂, θ₃) 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 k-th tile to the k-th 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

(62)

where is the average of all aperture statistics jackknife samples. The (mean-square) statistical uncertainty of the aperture statistics ⟨𝒩^a𝒩^bℳ⟩(θ_i) is . The inverse of this covariance, needed for the likelihood analysis, is estimated with

(63)

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 line-of-sights.

5.2. 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, ⟨𝒩^a𝒩^aℳ⟩(θ) and ⟨𝒩^b𝒩^bℳ⟩(θ), and cross-correlation 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,

(64)

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,

(65)

Optimal parameters p_opt are those that minimise the goodness-of-fit

(66)

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

(67)

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

(68)

with the χ² 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

(69)

where the matrix F is the Fisher information,

(70)

(e.g. Tegmark et al. 1997), and choose q as

(71)

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

(72)

normalised such that their sum is unity.

For each parameter pⁱ we give the α credibility interval (CI) on the marginalised posterior P(pⁱ|d^ab), defined by

(73)

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,

(74)

To find this interval, we sort the sampling points by the distance |pⁱ − | 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 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¹¹ and 10¹⁵ h⁻² 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

(75)

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.

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. 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 χ² = 93.63 for 90 − 14 = 76 degrees of freedom (d.o.f.), or χ²/d.o.f. = 1.28. This χ² 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. 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^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³ deg², 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¹¹ M_⊙ (−0.05 × 10¹¹ 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.

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 χ²/d.o.f. = 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: is roughly ten times larger than (68% credibility intervals, CI herafter). The per-halo 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 non-zero for halo masses above 10¹² M_⊙ (5 × 10¹¹ 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 ⟨𝒩^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^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 versus at 68% CI), while the correlation amplitude at 10¹² 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¹³ M_⊙ where 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.

6.2. 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 ⟨𝒩^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 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, χ², 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 p-values, 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 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 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 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 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 ⟨N|m⟩ with a local minimum near 6 × 10¹¹ M_⊙, 1 × 10¹² M_⊙, and 2 × 10¹² 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.1 in the Appendix. For all combinations of stellar-mass 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).

After the science verification, we fit the model to the observed G3L signal of KV450×GAMA. The best-fitting model parameters, χ²-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 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¹¹ M_⊙, m1 galaxies are the most numerous, whereas m4 galaxies dominate halos with masses above 4 × 10¹⁴ M_⊙. The most massive m5 galaxies, however, are never the most numerous sample between 10¹¹ to 10¹⁵ 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 m1-galaxies. Although they are observed up to z ∼ 0.34, for z ≳ 0.15 only the most massive m1 galaxies are seen (M_* > 2 × 10⁹ h⁻² 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 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 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 is more likely to be a than b, if a is a sample with higher stellar mass than b. As the four α^as 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 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.

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⁻² M_⊙ are positively correlated. Since r^ab increases with halo mass, the correlations between low stellar masses become relevant beyond galaxy-group size halos, where , , and for m = 10¹³ 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

(76)

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¹⁴ 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¹⁴ 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 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 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 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.

7. 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³deg² 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²) 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.5h⁻² M_⊙ and 10^11.5h⁻² M_⊙ (m1 to m5). The best-constrained HOD parameter for the KV450 × GAMA is the threshold halo mass , which is M_th ≈ 10¹¹ M_⊙ for our blue galaxies and M_th ≈ 10¹² 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 ⟨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. 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 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 ( for m ≳ 10¹³ 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¹⁰ M_⊙) and samples m3 or m4 (stellar masses between 10¹⁰ M_⊙ and 10¹¹ 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 low-occupancy halos. This finding implies that the assumption of uncorrelated satellite distributions is probably still appropriate if only galaxies in low-mass halos, m ≲ 10¹³ 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 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 cluster-sized halos (m ≳ 10¹⁴ 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.

It could be interesting to use G3L 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). 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⁻² 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⁻² 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 case-by-case 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 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 best-fitting 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 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 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 redshift-weighting 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 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. Second, we expect the parameter constraints to improve by combining G3L with other statistics, breaking degeneracies between HOD parameters (Sect. 3.3). For example, second-order galaxy-galaxy 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 second-order 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 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 large-scale 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 high-dimensional 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.

---

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/15-1. 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.A-3016, 177.A-3017, 177.A-3018, 179.A-2004, 298.A-5015. We also use products from the GAMA survey. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by several independent survey programmes, including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP, providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/. 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

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 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 sub-halos, 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,

(A.1)

For two identical galaxy samples a = b, this expression simplifies to

(A.2)

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,

(A.3)

where is the mean halo number density⁴. The PDF is a product of the normalised spatial distributions of satellites in each halo,

(A.4)

(A.5)

under the assumption of independent satellite positions inside a halo. The joint probability P_N of per-halo galaxy numbers has the moments

(A.6)

where p, q, r, s ∈ ℕ₀.

After inserting Eqs. (49) and (51) into (46) and using Eq. (A.1), the bispectrum for mixed pairs, a ≠ b, is given by

(A.7)

where unconnected terms do not contain δ_D(k₁ + k₂ + k₂), 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

(A.8)

where is the Kronecker symbol.

We divide the sum over i, j and k into the 1-halo term ₁ for i = j = k, the 2-halo term ₂ for i = j ≠ k, i = k ≠ j, and j = k ≠ i, and the 3-halo term ₃ for i ≠ j ≠ k. In the following, we calculate these terms individually.

A.1. 1-halo term

To calculate the 1-halo term ₁ 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

(A.9)

Evaluating all m-integrals aside from the one over m_i (H − 1 integrals in total) leads to

(A.10)

Using Eqs. (A.6), ∫ⅆ[3]x ug^a(x | m) = 1, and , we find

(A.11)

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,

(A.12)

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.13)

and hence

(A.14)

In this limit (assumed henceforth),

(A.15)

which leads to

(A.16)

(A.17)

with k₃ = −k₁ − k₂ and

(A.18)

For mixed pairs, a ≠ b, an analogous calculation gives

(A.19)

(A.20)

In the main text, we denote ₁(k₁, k₂, k₃) as ₁(k₁, k₂) for brevity.

A.2. 2-halo term

The 2-halo term ₂ 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

(A.21)

We now evaluate all m-integrals 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

(A.22)

In the limit of an infinite volume V, Eq. (A.14), the halo two-point correlation function ξ_h and its Fourier transform P_h, gives

(A.23)

(A.24)

With this expression, neglecting unconnected terms proportional to δ_D(k₁) and δ_D(k₂) and the approximation H(H − 1)/H² ≃ 1, we find

(A.25)

(A.26)

where again k₃ = −k₁ − k₂ and

(A.27)

For mixed pairs, a ≠ b, a similar calculation yields

(A.28)

In the main text we just denote ₂(k₁, k₂).

A.3. 3-halo term

The 3-halo term for unmixed pairs a = b is given by

(A.29)

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

(A.30)

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

(A.31)

(A.32)

As before, all terms not proportional to δ_D(k₁ + k₂ + k₃) are unconnected and therefore neglected. Using this expression and the approximation H(H − 1)(H − 2)/H³ ≃ 1, we obtain

(A.33)

(A.34)

For mixed pairs a ≠ b a similar calculation yields

(A.35)

Again, we just denote ₃(k₁, k₂) 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 non-Poisson 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 halo-mass range, such as between 2 × 10¹¹ M_⊙ ≲ m ≲ 10¹² 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 ⟨𝒩^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

(B.1)

where σ²(|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 1-halo term, containing ⟨( − 1)|m⟩, dominates here. Smaller scales are dominated by the halo term containing ⟨|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, that is, r^ab = 1 for perfectly correlated and , 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. 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 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 χ² of the fits. All fits have 76 d.o.f. (the same as the fit to the colour-selected samples in Sect. 6.1), so a χ² < 97.4 signifies an 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 χ² and p-values of the fits are given in Table C.2. Figure C.4 shows the model fits to the G3L aperture statistics from stellar mass-selected lenses in KV450× GAMA and the χ² of the fits. The parameters of the best-fitting models are given in Tables C.3 and C.4.

All Tables

All Figures