© ESO 2020

1. Introduction

According to the hierarchical galaxy formation model, galaxy groups and clusters form by the accretion of isolated galaxies and groups. This type of assembly process tidally strips mass from the infalling satellite galaxies and haloes. Because the dark matter is dissipationless (to a good approximation), it is more easily stripped from the subhalo than the baryons, which dissipate some of their energy and sink to the centre of their potential more efficiently than the dark matter, well before forming stars (White & Rees 1978). Because the dark matter is not that centrally concentrated, it is thus more susceptible to tidal stripping than the baryons, even after a galaxy forms its stars, which are, to the first order, dissipationless as well. This model thus predicts that the satellite galaxies will be preferentially stripped of their dark matter and the effect can be observed as higher stellar mass to halo mass ratios of satellite galaxies compared to their central counterparts of a similar stellar mass. While some stars may be lost, relatively more dark matter will be stripped and the result is a higher stellar-to-halo mass (SMHM) ratio for satellites in dense environments, compared to centrals and/or less dense environments. Previous simulation studies (see for example Bower et al. 2006) show that the SMHM relation of satellite galaxies is significantly different from the SMHM relation of central galaxies. Further evidence for different SMHM relations for central and satellite galaxies comes from abundance matching methods, which show that the satellites typically have more stellar mass than the central galaxies for a given halo mass (Reddick et al. 2013), and they have similar or in some cases even larger stellar masses than at the infall time (Reddick et al. 2013; Wechsler & Tinker 2018, and the references therein), while their halo masses do not.

While the SMHM relation of central galaxies has been successfully measured by many studies (for instance by Hoekstra et al. 2005; Mandelbaum et al. 2006; More et al. 2011; van Uitert et al. 2011; Leauthaud et al. 2012), this is not the case for satellite galaxies whose SMHM relation remains essentially unconstrained (Sifón et al. 2018). Recently, several weak gravitational lensing studies using galaxy groups and clusters have been undertaken (such as the ones by Limousin et al. 2007; Li et al. 2014, 2016; Sifón et al. 2015, 2018), all of them finding that the satellite galaxies’ haloes are heavily truncated with respect to the central and field ones.

Weak gravitational lensing, through the lensing of background sources by a sample of galaxies, which is commonly called galaxy-galaxy lensing, directly measures the total mass of lensing galaxies, without assuming their dynamical state (Bartelmann & Schneider 2001; Courteau et al. 2014), and it is currently the only method available to measure the total mass of samples of galaxies directly. Measuring the lensing signal around satellite galaxies, however, can be particularly challenging for several reasons: their small contribution to the lensing signal by the host galaxy group, source blending at small separations, and sensitivity to field galaxy contamination (Sifón et al. 2018). As pointed out by Sifón et al. (2015), the latter point is rather important as the field galaxies are not stripped and their contamination therefore complicates the interpretation of the lensing signal.

These studies were based on tangentially averaging the shear, which washes out information for satellites to some extent. In Dvornik et al. (2019), we revisit the two-dimensional galaxy-galaxy lensing method (Schneider & Rix 1997; Heymans et al. 2006). This method, which tries to fit a two-dimensional shear field directly to the galaxy ellipticity measurements, was shown on simulated data to perform significantly better than the traditional one-dimensional analysis of stacked tangential shear profiles or the closely related excess surface density (ESD). One important advantage of the two-dimensional method lies in the fact that it exploits all of the information regarding the actual image configuration (the model predicts the shear for each individual background galaxy image) using the galaxies’ exact positions, ellipticities, magnitudes, luminosities, stellar masses, group membership information, etc., rather than only using the ensemble properties of statistically equivalent samples (Schneider & Rix 1997). Moreover, the clustering of the lenses is naturally taken into account, although it is more difficult to account for the expected diversity in density profiles (Hoekstra 2014).

This method went out of fashion due to the unavailability of galaxy grouping information that would accurately classify galaxies as centrals and satellites (Hoekstra 2014), that is, the same information needed to robustly study the stellar mass to halo mass relation of satellite galaxies (Sifón et al. 2015). Treating the galaxies as centrals and satellites in a statistical way when considering the stacked signal could be naturally accounted for with the halo model (Seljak 2000; Peacock & Smith 2000; Cooray & Sheth 2002), thus overcoming the observational shortcomings. In recent years, this type of galaxy grouping information has become available thanks to the power of overlapping wide-field photometric surveys with highly complete spectroscopic surveys which allow one to treat central and satellite galaxies deterministically. The Kilo-Degree Survey (KiDS, Kuijken et al. 2015; de Jong et al. 2015) in combination with the overlapping Galaxy And Mass Assembly survey (hereafter GAMA, Driver et al. 2011; Robotham et al. 2011) provide an optimal data set for this type of analysis.

In this paper we present a two-dimensional galaxy-galaxy lensing measurement of the stellar-to-halo mass relation for central and satellite galaxies by combining a sample of spectrocopically confirmed galaxy groups from the GAMA survey and background galaxies from the fourth data release of KiDS (Kuijken et al. 2019). We use these measurements to constrain the stellar-to-halo mass relation, comparing the standard one-dimensional stacked tangential shear profiles with the two-dimensional galaxy-galaxy lensing method from Dvornik et al. (2019).

The outline of this paper is as follows. In Sect. 2 we present the lens and source samples used in this analysis. In Sect. 3 we present the specific lens model used in the paper and in Sect. 4 we describe the two-dimensional galaxy-galaxy lensing formalism. The parameter inference procedure is presented in Sect. 5. We show the results in Sect. 6, compare our results with the literature in Sect. 8, and conclude with Sect. 9. Throughout the paper we use the following cosmological parameters, which enter into the calculation of the distances and other relevant properties (Planck Collaboration XVI 2013): Ω_m = 0.307, Ω_Λ = 0.693, σ₈ = 0.8288, n_s = 0.9611, Ω_b = 0.04825, and h = 0.6777. The halo masses are defined as $M=4\pi r_{\mathrm{\Delta }}^3\mathrm{\Delta }{\overline{\rho }}_m/3$, the mass enclosed by the radius r_Δ within which the mean density of the halo is Δ times the mean density of the Universe ${\overline{\rho }}_m$, with Δ = 200. All of the measurements presented in the paper are in co-moving units.

2. Data and sample selection

The foreground galaxies used in this lensing analysis are taken from GAMA, a spectroscopic survey carried out on the Anglo-Australian Telescope with the AAOmega spectrograph. Specifically, we use the information of GAMA galaxies from three equatorial regions, G9, G12, and G15 from GAMA II (Liske et al. 2015). We do not use the G02 and G23 regions, as G02 does not overlap with KiDS and G23 uses an inconsistent target selection. These equatorial regions encompass ∼180 deg², contain 180 960 galaxies (with nQ ≥ 3, where the nQ is an indicator of redshift quality), and are highly complete down to a Petrosian r-band magnitude of r = 19.8. We make use of the GAMA galaxy group catalogue by Robotham et al. (2011), which provides information on the galaxy’s group membership which is used to separate them into central and satellite galaxies. The GAMA galaxy group catalogue was constructed using a three-dimensional Friends-of-Friends (FoF) algorithm, linking galaxies in projected and line-of-sight separation. We use version 10 of the group catalogue (G3Cv10), which contains 26 194 groups with at least two members. All of the galaxies that are not grouped in any of the 26 194 groups are considered as centrals, which is shown to be a correct assumption in Brouwer et al. (2017). The GAMA survey is 98% complete down to the observed magnitude limit. The group catalogue purity reaches 90% for high multiplicity groups (Robotham et al. 2011), out of which 70–75% of centrals are correctly identified (Robotham et al. 2011; Sifón et al. 2015). We consider all galaxies whose stellar mass is between 10⁸ M_⊙ and 10¹² M_⊙. Stellar masses are taken from version 20 of the LAMBDAR stellar mass catalogue, described in Wright et al. (2017). The final selection of galaxies can be seen in Figs. 1 and 2, and all of the relevant properties we need in our analysis are presented in Table 1. The stellar mass binning is only used for the one-dimensional galaxy-galaxy lensing case in order to obtain stacks of tangential shear signal, and it is chosen in such way that we have a similar signal-to-noise ratio in each stellar mass bin. In the two-dimensional case, we directly use the relevant individual galaxy quantities in the model.

Fig. 1. Stellar mass versus redshift of galaxies in the equatorial regions of the GAMA survey that overlap with KiDS. The full sample is shown with the hexagonal density plot and the dashed lines show the cuts for the stellar mass bins used in our analysis.

Fig. 2. Stellar mass distributions in our six bins used for one-dimensional stacked tangential shear measurements. The exact bin values are presented in Table 1.

We use imaging data from the 180 deg² of the fourth KiDS data release (Kuijken et al. 2019) that overlaps with the three equatorial patches of the GAMA survey to obtain shape measurements of background galaxies. KiDS is a four-band imaging survey conducted with the OmegaCAM CCD mosaic camera mounted at the Cassegrain focus of the VLT Survey Telescope (VST); the camera and telescope combination provide us with a fairly uniform point spread function across the field-of-view. The companion VISTA-VIKING (Edge et al. 2013) survey has provided complementary imaging in near-infrared bands (Z, Y, J, H, K_s), resulting in a deep, wide, nine-band imaging dataset (Wright et al. 2019).

We use shape measurements based on the r-band images, which have an average seeing of 0.66 arcsec. The image reduction, photometric redshift calibration, and shape measurement analysis is described in detail in Hildebrandt et al. (2020) and Kuijken et al. (2019). We measure galaxy shapes using lensfit (Miller et al. 2013), which has been calibrated using image simulations described in Kannawadi et al. (2019). This provides galaxy ellipticities (ϵ₁, ϵ₂) with respect to an equatorial coordinate system, and an optimal weight.

3. Lens model

The most widely assumed density profile for dark matter haloes is the Navarro–Frenk–White (NFW) profile (Navarro et al. 1996). Using simple scaling relations, this profile can be matched to simulated dark matter haloes over a wide range of masses and was found to be consistent with observations (Navarro et al. 1996). The NFW profile is defined as:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{NFW}}(r)=\frac{{\delta }_{\mathrm{c}}{\overline{\rho }}_m}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2},\end{array}$$(1)

where the free parameters δ_c and r_s are called the overdensity and the scale radius, respectively, r is the radius, and ${\overline{\rho }}_m$ is the mean density of the Universe, where ${\overline{\rho }}_m={\mathrm{\Omega }}_m{\rho }_{\mathrm{c}}$ and ρ_c is the critical density of the Universe, defined by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{c}}\equiv \frac{3H_0^2}{8\pi G},\end{array}$$(2)

where H₀ is the present-day Hubble parameter.

Some thought is warranted when choosing how to model stripped satellites galaxies. In numerical simulations, the satellite galaxies are heavily stripped by their host halo, but the effect of stripping on their density profile is not that severe. Even though tidal stripping removes mass from the outskirts of the halo, tidal heating causes the subhalo to expand, and the resulting density profile is similar in shape to that of a central galaxy which has not been subject to tidal stripping (Hayashi et al. 2003). Similarly, Pastor Mira et al. (2011) found that the NFW profile is a better fit than truncated profiles for subhaloes in the Millennium Simulation (Springel et al. 2005), and that the reduction in mass produced by tidal stripping is simply reflected as a change in the NFW concentration of subhaloes. Following Sifón et al. (2018), we have decided to model both centrals and satellites using the NFW profile, but while allowing the concentrations to differ.

The NFW profile in its usual parametrisation has two free parameters for each halo, halo mass M_h, and concentration c, and using those is the conventional way of modelling halo profiles. However, having two free parameters for each halo is computationally very expensive. Instead, we would like to describe these parameters through relations that depend on halo properties and then fit to a few free parameters in these global relations instead of hundreds or thousands of free, halo-specific parameters. To do so, we adopt the halo mass – concentration relation of Duffy et al. (2008), with the free concentration normalisation f_c:

$$\begin{array}{c}\hfill c(M_{\mathrm{h}},z)=f_{\mathrm{c}}10.14{\left[\frac{M_{\mathrm{h}}}{(2\times {10}^{12}{M}_{\odot }/h)}\right]}^{-0.081}{(1+z)}^{-1.01}.\end{array}$$(3)

We describe the stellar mass to halo mass relation as an exponential function:

$$\begin{array}{c}\hfill M_{\mathrm{h}}/M_{\odot }={\left(\frac{\alpha -log(M_{\star }/M_{\odot })}{e^{\gamma }}\right)}^{\beta },\end{array}$$(4)

where¹α = 12.0, and β and γ are the free parameters we will be fitting. We note that the functional form presented here stems from Matthee et al. (2017), but we have redefined some of the quantities². We use separate relations for the central and satellite galaxies, thus we have two sets of β and γ parameters since we want to constrain the SMHM relation for those populations separately. The choice of this type of parametrisation for the SMHM relation is motivated by reducing the number of free parameters required for the fit, while to first order maintaining the shape of the relation that is similar to the more widely adopted double power law parametrisations (as presented by Leauthaud et al. 2012; Moster et al. 2013; van Uitert et al. 2016).

The gravitational shear and convergence profiles are calculated using Eqs. (14)–(16) from Wright & Brainerd (2000), from which the predicted ellipticities for all of the lenses are calculated according to the weak lensing relations presented in Schneider (2003). We first calculate the reduced shear for our NFW profiles:

$$\begin{array}{c}\hfill g(x_i,z_{\mathrm{s}})=\frac{\gamma (x_i,z_{\mathrm{s}})}{1-\kappa (x_i,z_{\mathrm{s}})},\end{array}$$(5)

from which the ellipticities are calculated according to the following equation:

$$\begin{array}{c}\hfill ϵ=\{\begin{array}{cc}g\hfill & \text{if}|g|\le 1\hfill \\ 1/g^{\ast }\hfill & \text{if}|g|>1,\hfill \end{array}\end{array}$$(6)

where we have assumed that the intrinsic ellipticities of the sources average to 0 due to their random nature. Intrinsic alignments are not thought to contribute significantly to the signal at the current signal-to-noise ratio (Blazek et al. 2012).

We compute the effective critical surface mass density that we need in our lens model for each lens using the spectroscopic redshift of the lens z_l and the full normalised redshift probability density of the sources, n(z_s), calculated using the direct calibration method presented in Hildebrandt et al. (2017, 2020). The effective inverse critical surface density³ can be written as:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{cr},\mathrm{ls}}^{-1}=\frac{4\pi G}{c^2}{(1+z_{\mathrm{l}})}^{2}D(z_{\mathrm{l}}){\displaystyle {\int }_{z_{\mathrm{l}}}^{\infty }\frac{D(z_{\mathrm{l}},z_{\mathrm{s}})}{D(z_{\mathrm{s}})}n(z_{\mathrm{s}})\mathrm{d}{z}_{\mathrm{s}},}\end{array}$$(7)

where D(z_l) is the angular diameter distance to the lens, D(z_l, z_s) is the angular diameter distance between the lens and the source, and D(z_s) is the angular diameter distance to the source.

The galaxy source sample is specific to each lens redshift with a minimum photometric redshift z_s = z_l + δ_z, with δ_z = 0.2, where δ_z is an offset to mitigate the effects of contamination from the group galaxies (for details, see also the methods section and Appendix of Dvornik et al. 2017). We determine the source redshift distribution n(z_s) for each sample by applying the sample photometric redshift selection to a spectroscopic catalogue that has been weighted to reproduce the correct galaxy colour-distributions in KiDS (for details see Hildebrandt et al. 2020). The accuracy of this method, which is determined through mock data analysis, is sufficient for our study (Wright et al. 2019). We correct the measured ellipticities for the multiplicative shear bias per source galaxy per redshift bin as defined in Hildebrandt et al. (2020) with the (small) correction values estimated from image simulations (Kannawadi et al. 2019).

4. Galaxy-galaxy lensing formalism

In this study of satellite galaxy-galaxy lensing, we use the two-dimensional galaxy-galaxy lensing formalism as presented in Dvornik et al. (2019); we follow their model and adapted it to KiDS+GAMA by taking into account the survey’s specific requirements. Generally, for both the one-dimensional and two-dimensional cases, the likelihood of a model with a set of parameters θ given data d can be parametrised in the following form:

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{d}|\mathit{\theta })=\frac{1}{\sqrt{{\left(2\pi \right)}^n|\mathbf{C}|}}exp[-\frac{1}{2}{(\mathbf{m}(\mathit{\theta })-\mathbf{d})}^T{\mathbf{C}}^{-1}(\mathbf{m}(\mathit{\theta })-\mathbf{d})],\end{array}$$(8)

where m(θ) is the value of d predicted by the model with parameters θ. We assume the measured data points d = [d₁, …, d_n] are drawn from a normal distribution with a mean equal to the true values of the data, and n is the dimensionality of the data. In principle, the likelihood does not need to be Gaussian, but in practice it is a very good approximation due to ellipticity distribution being nearly Gaussian as well. The likelihood function accounts for correlated data points through the covariance matrix C of shape n × n. The covariance matrix C generally consists of two parts, the first arising from shape noise and the second from the presence of a cosmic structure between the observer and the source (Hoekstra 2003):

$$\begin{array}{c}\hfill \mathbf{C}={\mathbf{C}}^{\mathrm{shape}}+{\mathbf{C}}^{\mathrm{LSS}}.\end{array}$$(9)

In the case when one wants to fit one-dimensional tangential shear profiles, which are stacked over a sample of lenses, the likelihood function can be written as:

$$\begin{array}{cc}& \mathcal{L}(g_t^{\mathrm{obs}}|M_{\mathrm{h}},M_{\star },c)\hfill \\ \hfill & ={\displaystyle \prod _{i=1}^n\frac{1}{{\sigma }_{g_t,i}\sqrt{2\pi }}exp[-\frac{1}{2}{\left(\frac{g_{t,i}(M_{\mathrm{h}},R_i,z)-g_{t,i}^{\mathrm{obs}}}{{\sigma }_{g_t,i}}\right)}^2],}\hfill \end{array}$$(10)

where we have used m_i = g_t, i(M_h, R_i, z) (see Eq. (5)) as the model prediction given halo mass M_h, radial bin R_i, and redshift of the lens z, and $d_i=g_{t,i}^{\mathrm{obs}}$ as the tangentially averaged reduced shear of a sample of lenses measured from observations. The halo mass M_h and stellar mass M_⋆ are connected through Eq. (4), and the concentration c is defined in Eq. (3). Here we have also used the uncertainty of our measurement, given by the σ_gt, i calculated from the intrinsic shape noise of sources in each radial bin. The product runs over all radial annuli i. Moreover we only account for the diagonal terms of the covariance matrix and only include the error due to the shape noise, that is to say

$$\begin{array}{c}\hfill \sqrt{|C|}={\displaystyle \prod _{i=1}^n{\sigma }_i.}\end{array}$$(11)

This is done to reduce the computational complexity of the problem, and it is justified due to the covariance matrix being shape noise dominated. The above equation describes the one-dimensional method; the two-dimensional method differs only in the following significant way:

$$\begin{array}{cc}& \mathcal{L}({ϵ}^{\mathrm{obs}}|M_{\mathrm{h}},M_{\star },c)\hfill \\ \hfill & ={\displaystyle \prod _{i=1}^n\frac{1}{{\sigma }_{ϵ,i}\sqrt{2\pi }}exp[-\frac{1}{2}{\left(\frac{g_i(M_{\mathrm{h}},x_i,z)-{ϵ}_i^{\mathrm{obs}}}{{\sigma }_{ϵ,i}}\right)}^2],}\hfill \end{array}$$(12)

where g_i(M_h, x_i, z) are the model reduced shears evaluated at each source position x_i, ${ϵ}_i^{\mathrm{obs}}$ is the observed elipticities of real galaxies, and σ_ϵ, i is the intrinsic shape noise of our galaxy sample per component, calculated from the lensfit weights following the description by Heymans et al. (2012). The same lensfit weights are used to weight ${ϵ}_i^{\mathrm{obs}}$ as well (Heymans et al. 2012). In practice, the two-dimensional fit to the ellipticities is carried out for each Cartesian component of ellipticity ϵ₁ and ϵ₂ with respect to the equatorial coordinate system. Here the product runs over all of the individual source galaxies i.

5. Parameter inference procedure

Due to the computational complexity of the analysis, we do not use the Markov chain Monte Carlo (MCMC) method for parameter inference for any of our methods. Our model parameter inference procedure and the fit to the data is performed using three steps for both the one-dimensional and two-dimensional cases. We first sample 100 points using the latin hypercube method (McKay et al. 1979) in the six-dimensional model parameter space (see Table 2 for the ranges of all of the parameters). We picked these types of ranges for the parameters in order to sample the likelihood surface in the 5σ range that we found using an MCMC fit of the 1D model in preliminary tests. This minimises the need for having a larger number of points in the latin hypercube as well as reducing the number of interpolation points at the later step. For each of the 100 points, we calculate the likelihood value ℒ(d | θ) according to Eqs. (10) and (12). We created 10 000 realisations of the latin hypercube and select one that maximises the euclidean distance between all of the points (similarly to what it was done by Heitmann et al. 2009). We have also checked the influence of different realisations on the obtained results by evaluating the one-dimensional method on different realisations, and we always obtain the same resulting likelihoods.

The second step requires the construction of an interpolator, for which we use the Gaussian process (GP) regression method with a multi-dimensional radial basis function (RBF) kernel⁴. The interpolated value is the likelihood ℒ(d | θ) as obtained from the latin-hypercube samples. Interpolation is performed with the same ranges as used in the construction of the latin hypercube using 20⁶ (64 million) equally spaced grid points. We test the accuracy of the interpolation by choosing one point as a test point, and we use the remaining 99 likelihoods to construct the interpolator. We then compare the prediction at the test point to the actual likelihood value. This is repeated for all likelihood values of our latin hypercube. We show the fractional differences between the interpolated likelihood and true likelihood value at each point in Fig. 3. The majority of sampled points are accurate to better than 1%. The use of the latin hypercube to construct the initial grid to sample the likelihood surface does not increase the uncertainty in recovering the true values due to the properties of latin hypercube sampling compared to the usual grid search minimisation (McKay et al. 1979).

Fig. 3. Fractional differences between the interpolated likelihood and true likelihood value at each sampled point in the latin hypercube. We do not show the cases where one of the interpolated points is exactly on one of the edges of the latin hypercube since the interpolator is unable to properly perform for those edge cases. The remaining high value outliers are the points in the latin hypercube that are close to the edges of the parameter space.

In the third step, we calculate the marginalised distributions from the interpolated points. First we normalise the probability grid P_6D, such that:

$$\begin{array}{c}\hfill {\displaystyle \sum _{\mathit{\theta }}{P}_{6\mathrm{D}}(\mathbf{d}|\mathit{\theta })=\sum _{\mathit{\theta }}\mathcal{L}(\mathbf{d}|\mathit{\theta })=1,}\end{array}$$(13)

which also sets the normalisation term in Bayes’ theorem, so we can, from the resulting values, calculate the one-dimensional and two-dimensional marginalised distributions of all of the parameters.

6. Results

We fit the lens model as described in Sect. 3 to the stacked tangential shear measurements in our six stellar mass bins (our 1D result) and to the full two-dimensional shear field (our 2D result). An example of a single stacked tangential shear profile for the GAMA lenses in the 10^10.5–10^10.75 M_⊙ stellar mass bin is shown in Fig. 4, with the measurements and their respective 1σ errors⁵. The measured lens model best-fit parameters (median of the marginalised posterior estimate), together with their 68% credible intervals are presented in Table 2 for both the one-dimensional and two-dimensional analysis. The constrains from the two approaches can be compared in terms of the full posterior distributions shown in Fig. 5. Even though none of our parameters are constrained to within 5σ of the preliminary test, the prior ranges are still good, given that models outside of the prior range for γ and β parameters would be unphysical.

Fig. 4. Stacked tangential shear profile for the GAMA lenses (blue points) in the 1010.5–1010.75 M⊙ stellar mass bin, compared to the best fitting lensing model for the 1D method, with contributions from both centrals and satellites. The orange band encloses the 68% credible interval.

Fig. 5. Full posterior distributions of the model parameters γcen, βcen, fc, cen, γsat, βsat, and fc, sat, for both the one-dimensional stacked tangential shear measurements (in blue) as well as the two-dimensional galaxy-galaxy lensing method (in orange). The contours indicate the 1σ and 2σ credible regions. We note that γ and β are the free parameters of the SMHM relation (Eq. (4)) and fc is the concentration-mass relation normalisation parameter (Eq. (3)).

For both methods, we find that the SMHM relations are completely described with two parameters each (two for centrals and two for satellites): the normalisation γ and slope β, for which the obtained one-dimensional values for centrals and satellites are presented in Table 2. In the same table, we present the obtained values for centrals and satellites in the case when using the two-dimensional galaxy-galaxy lensing. These constraints are consistent with those from the one-dimensional analysis with the highest discrepancy in the γ_cen parameter, which differs by 1σ. This shows that the two methods perform equally well statistically.

The normalisations of the concentration-halo mass relation are essentially unconstrained in the adopted prior range. The other values of the SMHM relation are not correlated with these values, and the prior does not influence our results in a substantive way. The wide range is somewhat driven by the imperfect modelling for both the one-dimensional and the two-dimensional case, as the model does not properly account for the two-halo term. The prior ranges are comparable to the values found in hydrodynamical simulations (Dvornik et al. 2019). These values are also consistent with the observational findings that prefer lower normalisations than expected in simulations, such as in the studies of Viola et al. (2015), Sifón et al. (2015), and Dvornik et al. (2017). Since there are no strong covariances between f_c and the other parameters, any small systematic error in f_c probably does not propagate through to a bias in other parameters of the SMHM relation.

We show a typical halo mass for a central and satellite galaxy with a stellar mass of log(M_⋆/M_⊙) = 10.6 in Fig. 6, obtained from propagating the best fit parameters through the SMHM relation. We find that the SMHM relations are different for the central and satellite galaxies, showing that the stripping of the dark matter does indeed take place; the SMHM relation of satellite galaxies is higher than the relation for the centrals, as is also seen in Fig. 6. We find that satellite galaxies are hosted by halo masses that are systematically 0.53 ± 0.39 dex (for 1D) and 0.23 ± 0.18 dex (for 2D) smaller than those of central galaxies at this stellar mass. The uncertainty of the inferred SMHM relation is similar to the intrinsic scatter present in simulations, for instance by the EAGLE hydrodynamical simulation (Schaye et al. 2015; Matthee et al. 2017, see Figs. 8 and 9). While we see the same qualitative conclusions between the one-dimensional and two-dimensional analysis, the quantitative halo masses inferred are inconsistent at the ∼1σ level.

Fig. 6. Halo masses for a galaxy with a stellar mass of log(M⋆/M⊙) = 10.6 for the one-dimensional and two-dimensional analyses for both central and satellite galaxies.

7. Assessment of completeness bias in 2D galaxy-galaxy lensing

As is shown in Dvornik et al. (2019), the two-dimensional analysis relies on a complete sample of lenses. If lenses are missing from the model, the bias in halo mass can be as much as 20%. Although we used all of the galaxies with redshifts from the GAMA data, at a given redshift, the magnitude limit implies a limit in stellar mass. Galaxies with that stellar mass but at higher redshifts are not included in the catalogue, but they do contribute to the lensing signal. The contribution to the lensing signal is small for very high redshift lenses; however, not including lenses near the magnitude limit of GAMA may bias our measurements. Moreover, the process of labelling galaxies into centrals and satellites is not perfect. The GAMA group catalogue is only pure up to 90% for high multiplicity groups, and it is known to be contaminated by the misidentification of the central galaxy in a group to such an extent that only 70 – 75% of central galaxies are correctly identified (Robotham et al. 2011). Thus the true central galaxy would be included in the satellite sample, which can introduce a bias of about 15% in the inferred masses (Sifón et al. 2015). This effect is even more pronounced in the GAMA group catalogue for pairs of galaxies, where both components are likely to be centrals and not a central galaxy and one satellite galaxy. What is more, the satellite stellar-to-halo mass relation at a high stellar mass is possibly driven by this misidentification of satellite galaxies, which should actually be classified as centrals, given the high halo masses measured; satellite galaxies with stellar masses up to log(M_⋆/M_⊙) = 12 should not be common. This is a likely consequence of the observed problem with the FoF algorithm used to identify galaxy groups in the GAMA survey, but it does not seem to substantially affect the results. The FoF algorithm separates groups into a number of smaller groups or smaller, aggregate, unrelated groups into one large group, which would then host more than one central galaxy with them being classified as a satellite (Jakobs et al. 2018).

In order to assess the possible bias due to missing galaxies in a magnitude limited survey, such as GAMA, we use the MICE-GC N-body simulation from which the MICE collaboration constructed a lightcone spanning a full octant on the sky (Fosalba et al. 2015a,b). The MICE lightcone has a maximum redshift of 1.4. The haloes found in the simulation were populated using a hybrid halo occupation distribution (HOD) and halo abundance matching (HAM) prescription (Carretero et al. 2015; Crocce et al. 2015; Hoffmann et al. 2015. For our assessment, we use the MICECATv2.0 catalogue⁶, from which we take the positions of galaxies within a 4 deg² cutout of the lightcone with a redshift limit of z <  0.5 and an SDSS r-band magnitude of m_r <  22. We also select galaxies with a stellar mass between 10⁷ M_⊙/h and 10¹³ M_⊙/h. This selection of galaxies results in a distribution of stellar masses and redshifts similar to that of GAMA, which also has a similar number density. On this sample of galaxies, we apply an additional magnitude cut of m_r <  19.8, which is the magnitude limit of the GAMA survey (Driver et al. 2011).

We generate a noiseless mock shear field resembling a typical KiDS observation using the procedure shown in Dvornik et al. (2019). We populate the mock shear field with haloes at the locations of galaxies from MICE mocks by assigning the stellar-to-halo mass relation from Matthee et al. (2017), using the redshifts we have in the MICE mocks. We fit for the concentration and SMHM normalisations, using the two mock samples (the m_r <  22 and the m_r <  19.8 magnitude limited samples) as our input lenses for the fits. The parameter inference method is the same as described in Sect. 5. We show the comparison of the inferred parameters in Fig. 7, between the full sample of MICE galaxies (blue) and analysis for lens galaxies with m_r <  19.8 (orange). The model is able to accurately recover the input relation for the full sample of galaxies, while that is no longer the case for a magnitude limited sample. Due to a smaller number of lenses, the uncertainty increases but also the retrieved concentrations and halo masses are biased towards lower values. The effect is present at the 10% level, which is consistent with what we have found for a mock dataset of randomly placed lenses at a fixed redshift (Dvornik et al. 2019), but it is more representative of a real galaxy distribution and the observed effects due to the magnitude limit. The 10% change seems to be consistent with a statistical fluctuation, given that the recovered values lie within the 1σ contours. We do need to point out that this test is performed on noiseless data, for which we can choose the size of the contours, and the change in values reflects the true bias. Given the bias in the two-dimensional analysis and the comparable statistical performance of the two methods, our one-dimensional constraints are our preferred result for this analysis.

Fig. 7. Comparison of the inferred parameters between the full sample of MICE galaxies (blue) and for galaxies with mr <  19.8 (orange). The model is able to accurately recover the input relation for the full sample of galaxies (dashed line), while the magnitude limited sample is biased at the 10% level.

8. Comparison with previous studies

In Fig. 8 we show various published determinations of the relationship between the total and stellar mass of central galaxies (Leauthaud et al. 2012; Moster et al. 2013; Velander et al. 2014; Hudson et al. 2015; Zu & Mandelbaum 2015; van Uitert et al. 2016; Mandelbaum et al. 2016) as well as the EAGLE and Illustris TNG simulations (Engler et al. 2020). We scale all of these relations to our adopted values of H₀ and the definition of halo mass, that is to say the halo mass is defined with respect to 200 times the average density in the Universe, as is adopted in this paper. Furthermore, we also compare our results with the central and satellite properties from the hydrodynamical simulation EAGLE (Schaye et al. 2015; McAlpine et al. 2016). Specifically, we use the AGN model, from which we select galaxies with stellar masses ranging from 10^9.6 M_⊙ to 10^11.2 M_⊙ and their halo properties from which we can plot the mean SMHM relation and its scatter.

Fig. 8. Stellar-to-halo mass relation for the central galaxies compared to the various published determinations of the relationship between the total and stellar mass of central galaxies. All of the relations between the total and stellar mass of central galaxies are in broad agreement and are all using galaxy-galaxy lensing to constrain the SMHM relation, excluding the simulations where those quantities are measured directly from the particle data.

All of the relations between the total and stellar mass of central galaxies are in broad agreement and they all use galaxy-galaxy lensing to constrain the SMHM relation (in the case of Zu & Mandelbaum 2015; van Uitert et al. 2016, also in combination with galaxy clustering and the stellar mass function, respectively). For Velander et al. (2014) and Mandelbaum et al. (2016), we show their SMHM relation of red galaxies; as in the GAMA sample, central galaxies are mostly identified as red. Our measurements are also in agreement with the previous results. We also need to point out, as mentioned by Sifón et al. (2018), that the measurements from Zu & Mandelbaum (2015) and Mandelbaum et al. (2016) agree perfectly when the galaxies between the two samples are matched. The SMHM relations of Leauthaud et al. (2012), Moster et al. (2013), and Hudson et al. (2015) can be considered as better comparisons than the ones from the EAGLE or Illustris TNG simulations since they are obtained at redshifts comparable to the redshifts of the GAMA galaxies.

Similarly, for satellites, we show in Fig. 9 the comparison of the results from Rodríguez-Puebla et al. (2013), Sifón et al. (2018), and the EAGLE simulation with our two methods. However, the definitions of the halo mass of satellite galaxies in both Rodríguez-Puebla et al. (2013) and Sifón et al. (2018) are not equivalent to the one we use throughout this paper and it is also hard to correct this in order to compare the same quantities. Rodríguez-Puebla et al. (2013) use the definition of the subhalo mass that is defined as a mass of a satellite halo at the observed time (present time), and Sifón et al. (2018) use the definition of the subhalo mass as the mass within a radius for which the subhalo density matches the background density of the cluster at the distance of the subhalo in question. The closest definition to ours is the definition from the EAGLE simulation. Our results are similar to the behaviour of satellite galaxies therein. As for the comparison with Rodríguez-Puebla et al. (2013) and Sifón et al. (2018), all of the studies show lower satellite masses compared to the central galaxies at the same stellar mass. The same also holds true for the overall trend as a function of stellar mass. In Fig. 10 we also show the ratio between the satellite and central halo mass as a function of stellar mass. We observe that for the low stellar mass galaxies, the ratio is around 0.4 and drops towards 0.2 for high stellar mass galaxies, although the uncertainties on the ratio are quite large. This result directly shows us that more than ~80% of the dark matter of satellite galaxies is stripped (but with a large uncertainty), when they are accreted by a massive central galaxy.

Fig. 9. Stellar-to-halo mass relation for the satellite galaxies compared to the various published determinations of the relationship between the total and stellar mass of satellite galaxies. We show the SMHM relations of Leauthaud et al. (2012) and van Uitert et al. (2016) for comparison purposes with the centrals since they capture the majority of the other relations shown in Fig. 8.

Some studies of cosmological N-body simulations find that the ratio of satellite galaxy masses over central galaxy masses is M_h, sat/M_h, cen ≪ 1 (van den Bosch & Ogiya 2018, and the references therein), which means that the majority of dark matter haloes that merge with a central halo would be completely tidally disrupted. The tidal stripping and disruption of subhaloes in numerical simulations was mostly shown to be of numerical origin, emphasising that the cosmological simulations still suffer from overmerging (van den Bosch & Ogiya 2018). Even though our ratio of halo masses between the satellites and centrals show large uncertainties, with a refined analysis, it would be possible to use this type of measurement as an independent observational confirmation for the lack (or presence) of numerical convergence in the simulations. Given our results, as presented in Fig. 10, we cannot make a clear statement, although the results hint that less stripping is present in the observational data. On the other hand, the hydrodynamical simulations, such as EAGLE, show that the artificial stripping of haloes is not significant (Chaves-Montero et al. 2016), but we need to caution for the survivor bias in this case since completely disrupted haloes would not be present in the sample.

Fig. 10. Ratio between the satellite and central halo mass as a function of stellar mass.

All of the previous results, together with our findings show that the satellite galaxies are indeed preferentially stripped of their dark matter and the effect can be observed as higher stellar mass to halo mass ratios compared to their central counterparts of a similar stellar mass. All of the measurements show a statistically different SMHM relation for the central and satellite galaxies, which furthermore shows that two-dimensional galaxy-galaxy lensing can measure the SMHM relation for different populations of galaxies.

9. Discussion and conclusions

We have measured the stellar-to-halo mass relation of central and satellite galaxies in the GAMA survey. In this analysis, we use the more advanced two-dimensional galaxy-galaxy lensing method to constrain the SMHM relation, which has potential benefits over the traditionally used stacked tangential shear method (also referred to as one-dimensional galaxy-galaxy lensing here, Dvornik et al. 2019).

We use the three equatorial GAMA patches that overlap with the KiDS data in order to measure both the tangential shear signal around the central and satellite galaxies as well as the two-dimensional galaxy-galaxy lensing constraints on the same lenses and sources. The shear signals are then used to constrain the SMHM relation of central and satellite galaxies.

We model the lensing signal using an NFW profile together with the concentration-mass relation by Duffy et al. (2008), scaled by a normalisation factor for which we took the fit into account. We assume a functional form for the SMHM relation in the form of an exponential function, which was motivated by the observed behaviour in the simulations, and we fold it through our model, thus directly fitting for the normalisation and slope of the SMHM relation. The lens model is used to calculate the tangential shear profile, which is then fitted to the measured tangential shear profile from the GAMA and KiDS data, as well as to directly predict the two Cartesian components of the galaxies’ ellipticities used in our two-dimensional method.

We find that the SMHM relation can be successfully measured using the two-dimensional method, with a comparable statistical power to the traditional one-dimensional method using the stacked tangential shear measurements. Both methods give us similar results for the SMHM relations, showing that the two-dimensional method is indeed a robust way to measure properties of the galaxy–halo connection, without using statistically equivalent samples as in the case of the one-dimensional method, nor using more complicated halo models or relying on support from other probes. The resulting SMHM relations are broadly in agreement with the literature, and our results show that the satellite galaxies are indeed preferentially stripped of their dark matter and the effect can be observed as higher stellar mass to halo mass ratios compared to their central counterparts of a similar stellar mass.

By comparing the results of this paper with the findings of our previous paper (Dvornik et al. 2019), we are able to recognize that the comparable constraining power of the one-dimensional and two-dimensional method, shown in Fig. 5, is unexpected. In the Dvornik et al. (2019) paper, we predicted a factor of 3 improvement. As seen in the results, the statistical powers of both methods are comparable. Dvornik et al. (2019) explores an idealised mock dataset as well as noiseless and complete simulations, where the exact galaxy classification was also known. As mentioned in Sect. 6, multiple effects can and will cause differences from this idealised mock. First, the increased uncertainty of the two-dimensional method is directly dependent on how well one can identify central and satellite galaxies and how robust this identification and selection is. Even though the GAMA group catalogue (Robotham et al. 2011) is highly robust, this does not mean it is perfect and even a small number of incorrect classifications of galaxies would cause excessive scatter in the resulting SMHM relations, impacting the ability of the two-dimensional model to constrain the parameters. Secondly, the results of the two-dimensional galaxy-galaxy lensing seems to be biased due to the completeness limit of the GAMA survey. We have shown that slight incompleteness of the lens sample can cause biases in the inferred parameters, which can be as large as 10% (as is shown in Fig. 7). This can be somewhat seen in Fig. 5, where the γ_cen parameters are most noticeably different. The two-dimensional analysis is still computationally and resource expensive compared to the one-dimensional method. This further reduces the usability of the method, and our preferred galaxy-galaxy lensing analysis thus remains the standard one-dimensional approach.

Acknowledgments

We thank the anonymous referee for their very useful comments and suggestions. A.D., A.H.W. and H.Hi acknowledge support from ERC Consolidator Grant (No. 770935). HHi is further supported by a Heisenberg Grant of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1). H.Ho and A.K. acknowledge support from Vici Grant 639.043.512, financed by the Netherlands Organisation for Scientific Research (NWO). K.K. acknowledges support by the Alexander von Humboldt Foundation. M.B. is supported by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12, and by the Polish National Science Center through Grant no. 2018/30/E/ST9/00698. C.H., M.A., C.L., B.G. and T.T. acknowledge support from the European Research Council under grant number 647112, and C.H. further acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. T.T. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłlodowska-Curie grant agreement No 797794. This research is 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 and 177.A-3018, and on data products produced by Target/OmegaCEN, INAF-OACN, INAF-OAPD and the KiDS production team, on behalf of the KiDS consortium. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programs including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org. This work has made use of Python (http://www.python.org), including the packages numpy (http://www.numpy.org) and scipy (http://www.scipy.org). Plots have been produced with matplotlib (Hunter 2007). This work has made use of CosmoHub. CosmoHub has been developed by the Port d’Informació Científica (PIC), maintained through a collaboration of the Institut de Física d’Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institute of Space Sciences (CSIC & IEEC), and was partially funded by the “Plan Estatal de Investigación Científica y Técnica y de Innovación” program of the Spanish government. All authors contributed to writing and development of this paper. The authorship list reflects the lead authors (A.D., H.Ho, K.K.) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have made a significant contribution either to the data products or to the scientific analysis.

References

All Tables

All Figures

	Fig. 1. Stellar mass versus redshift of galaxies in the equatorial regions of the GAMA survey that overlap with KiDS. The full sample is shown with the hexagonal density plot and the dashed lines show the cuts for the stellar mass bins used in our analysis.
In the text

	Fig. 2. Stellar mass distributions in our six bins used for one-dimensional stacked tangential shear measurements. The exact bin values are presented in Table 1.
In the text

Fig. 3. Fractional differences between the interpolated likelihood and true likelihood value at each sampled point in the latin hypercube. We do not show the cases where one of the interpolated points is exactly on one of the edges of the latin hypercube since the interpolator is unable to properly perform for those edge cases. The remaining high value outliers are the points in the latin hypercube that are close to the edges of the parameter space.

In the text

	Fig. 4. Stacked tangential shear profile for the GAMA lenses (blue points) in the 1010.5–1010.75 M⊙ stellar mass bin, compared to the best fitting lensing model for the 1D method, with contributions from both centrals and satellites. The orange band encloses the 68% credible interval.
In the text

Fig. 5. Full posterior distributions of the model parameters γcen, βcen, fc, cen, γsat, βsat, and fc, sat, for both the one-dimensional stacked tangential shear measurements (in blue) as well as the two-dimensional galaxy-galaxy lensing method (in orange). The contours indicate the 1σ and 2σ credible regions. We note that γ and β are the free parameters of the SMHM relation (Eq. (4)) and fc is the concentration-mass relation normalisation parameter (Eq. (3)).

In the text

	Fig. 6. Halo masses for a galaxy with a stellar mass of log(M⋆/M⊙) = 10.6 for the one-dimensional and two-dimensional analyses for both central and satellite galaxies.
In the text

	Fig. 7. Comparison of the inferred parameters between the full sample of MICE galaxies (blue) and for galaxies with mr <  19.8 (orange). The model is able to accurately recover the input relation for the full sample of galaxies (dashed line), while the magnitude limited sample is biased at the 10% level.
In the text

Fig. 8. Stellar-to-halo mass relation for the central galaxies compared to the various published determinations of the relationship between the total and stellar mass of central galaxies. All of the relations between the total and stellar mass of central galaxies are in broad agreement and are all using galaxy-galaxy lensing to constrain the SMHM relation, excluding the simulations where those quantities are measured directly from the particle data.

In the text

	Fig. 9. Stellar-to-halo mass relation for the satellite galaxies compared to the various published determinations of the relationship between the total and stellar mass of satellite galaxies. We show the SMHM relations of Leauthaud et al. (2012) and van Uitert et al. (2016) for comparison purposes with the centrals since they capture the majority of the other relations shown in Fig. 8.
In the text

	Fig. 10. Ratio between the satellite and central halo mass as a function of stellar mass.
In the text