KiDS-1000: Combined halo-model cosmology constraints from galaxy abundance, galaxy clustering and galaxy-galaxy lensing

We present constraints on the flat $\Lambda$CDM cosmological model through a joint analysis of galaxy abundance, galaxy clustering and galaxy-galaxy lensing observables with the Kilo-Degree Survey. Our theoretical model combines a flexible conditional stellar mass function, to describe the galaxy-halo connection, with a cosmological N-body simulation-calibrated halo model to describe the non-linear matter field. Our magnitude-limited bright galaxy sample combines 9-band optical-to-near-infrared photometry with an extensive and complete spectroscopic training sample to provide accurate redshift and stellar mass estimates. Our faint galaxy sample provides a background of accurately calibrated lensing measurements. We constrain the structure growth parameter $S_8=\sigma_8\sqrt{\Omega_{\mathrm{m}}/0.3}=0.773^{+0.028}_{-0.030}$, and the matter density parameter $\Omega_{\mathrm{m}}=0.290^{+0.021}_{-0.017}$. The galaxy-halo connection model adopted in the work is shown to be in agreement with previous studies. Our constraints on cosmological parameters are comparable to, and consistent with, joint $3\times2{\mathrm{pt}}$ clustering-lensing analyses that additionally include a cosmic shear observable. This analysis therefore brings attention to the significant constraining power in the often-excluded non-linear scales for galaxy clustering and galaxy-galaxy lensing observables. By adopting a theoretical model that accounts for non-linear halo bias, halo exclusion, scale-dependent galaxy bias and the impact of baryon feedback, this work demonstrates the potential and a way forward to include non-linear scales in cosmological analyses. Varying the width of the satellite galaxy distribution with an additional parameter yields a strong preference for sub-Poissonian variance, improving the goodness of fit by 0.18 in reduced $\chi^{2}$ value compared to a fixed Poisson distribution.


Introduction
In the last quarter of the century we have seen the rise and establishment of the concordance cosmological model which describes the formation and subsequent evolution of the cosmic structure.In this concordance model the Universe at present time is modelled as a flat, cold dark matter (CDM) and cosmological constant (Λ) dominated Universe with a negligible contribution from neutrinos, and gravity described by the law of General Relativity.Furthermore, the initial power spectrum of density fluctuations is assumed to be a single power-law.Such ΛCDM models are described using only six free parameters, which govern the energy densities of baryons, Ω b , and cold dark matter, Ω dm , the spectral index, n s , and normalisation, σ 8 , of the density perturbations power spectrum, the Hubble parameter, H 0 , and the optical depth of reionisation.The flat geometry, implying that Ω Λ = 1 − Ω b − Ω dm , is strongly supported by high precision early-Universe measurements of the cosmic microwave background (CMB) temperature fluctuations combined with supernova and/or baryon acoustic oscillation distance measurements (Planck Collaboration et al. 2020;Scolnic et al. 2018;Alam et al. 2021).Formally the models also include the total neutrino mass, but the value of the parameter is too small for the current precision of observations (Gerbino & Lattanzi 2018).
A powerful probe of late-time cosmology is the large-scale distribution of galaxies.Even though the stars contribute negligibly to the overall energy density of the Universe, the light from stars in galaxies can be used to trace the evolution of the underlying distribution of dark matter in two complementary ways.Firstly, the light-path from distant galaxies is impacted by the distribution of foreground mass.This 'gravitational lensing' effect leads to a correlation between the observed shapes of galaxies, commonly referred to as cosmic shear.This observable can be used to probe the statistical properties of the total matter distribution in the Universe, typically quantified through the shape and amplitude of the matter power spectrum (Heymans et al. 2013;Hikage et al. 2019;Hamana et al. 2019;Asgari et al. 2021b;Secco et al. 2022;Amon et al. 2022a;van den Busch et al. 2022).Secondly, galaxies are expected to reside within dark matter haloes which form from the highest density peaks in the initial Gaussian random density field (e.g.Mo et al. 2010, and the references therein).Galaxies are therefore tracers of the underlying dark matter distribution, and with an accurate understanding of how biased these tracers are, the measurement of galaxy clustering as a function of redshift and scale provides strong constraints on the properties of the ΛCDM model (see for example Alam et al. 2021).It is becoming increasingly common to combine these two different 'two-point' (2pt) statistics, along with a third measurement of the gravitational lensing of background galaxies by foreground galaxies, otherwise known as galaxy-galaxy lensing.These joint '3×2pt' large-scale structure cosmological studies already have the precision to directly constrain some cosmological parameters independently of CMB measurements (Heymans et al. 2021;DES Collaboration et al. 2022).
In this analysis we focus on exploiting the significant precision recovered with small-scale measurements of galaxy clustering and galaxy-galaxy lensing.These non-linear scales are typically excluded from cosmological analyses owing to insufficient or uncertain modelling of the complex relationship between galaxies and the underlying matter distribution on these scales (Davis et al. 1985;Dekel & Rees 1987).Galaxy bias is scaledependent, stochastic and changes as a function of galaxy luminosity, colour and morphological type (Dekel & Lahav 1999;Zehavi et al. 2011;Cacciato et al. 2012;Dvornik et al. 2018).Based on these facts, it is not surprising that galaxy bias is generally considered a nuisance to be marginalised over in the recovery of cosmological constraints.Many studies limit their analyses to scales where the galaxy bias is considered to be linear and scale independent (see for example van Uitert et al. 2018;Yoon et al. 2019;DES Collaboration et al. 2022).An alternative approach uses perturbation theory (Desjacques et al. 2018), to expand galaxy bias modelling into the mildly non-linear regime (e.g.Mandelbaum et al. 2013;Sánchez et al. 2017;Heymans et al. 2021;Pandey et al. 2022).However, as a measurement of small-scale galaxy bias also contains a wealth of information regarding galaxy formation, we argue that it is preferential to utilise all the data, along with an appropriate galaxy bias model, to facilitate joint constraints on both cosmology and galaxy bias.
Given our previous attempts to shine the light on the galaxy bias and its properties (Dvornik et al. 2018), in this analysis we adopt a realistic and physically motivated halo model for galaxy bias.Under the assumption that all galaxies reside in dark matter haloes, we adopt a halo occupation distribution (HOD) model, a statistical description for how galaxies are distributed between and within the dark matter haloes (Peacock & Smith 2000;Scoccimarro et al. 2001;Mo et al. 2010;Yang et al. 2009;Cacciato et al. 2013;van den Bosch et al. 2013;Cacciato et al. 2014).When combined with the halo model, which describes the non-linear matter distribution as a sum of spherical dark matter haloes (Seljak 2000;Cooray & Sheth 2002), these models provide a fairly complete, broadly accurate1 , and easy to understand description of galaxy bias, halo masses, and galaxy clustering (Cacciato et al. 2013).
Our approach builds on the cosmological analysis presented in Cacciato et al. (2013); More et al. (2015), where the halo model is used to coherently analyse the clustering of galaxies, the galaxy-galaxy lensing signal (Guzik & Seljak 2002;Yoo et al. 2006;Cacciato et al. 2009) and the galaxy abundances as a function of luminosity or stellar mass (van den Bosch et al. 2013;Cacciato et al. 2013).Furthermore, the same approach was used to study the galaxy-halo connection exclusively, with a fixed cosmology, by Leauthaud et al. (2011); Coupon et al. (2015).This combination of probes is hereafter referred to as '2×2pt+SMF', representing the combination of the two twopoint statistics galaxy-galaxy lensing and galaxy clustering with the one-point stellar mass function.Analysis showing how much more information is in the 2×2pt+SMF compared to the 2×2pt is nicely presented in the paper from More et al. (2013).The degeneracy breaking of model parameters arising from the inclusion of a stellar mass function is shown in Mahony et al. (2022).
Since early applications, there has been significant interest in using the halo model to interpret large-scale structure probes Seljak et al. (2005); Cacciato et al. (2009); Li et al. (2009).The analysis of 2×2pt statistics, down to non-linear scales, has been shown to lead to tight constraints for both Ω m and σ 8 (Cacciato et al. 2013;Mandelbaum et al. 2013;More et al. 2015;Wibking et al. 2019).The halo model can also constrain extensions to the standard ΛCDM cosmologies, such as the equation of state of dark energy and neutrino mass (More et al. 2013;Krause & Eifler 2017).The choice of observables is motivated by the focus on a feasibility study on smaller scales which achieves similar precision, thus allowing for a direct and/or independent comparison to cosmic shear studies.In the era of high-precision cosmology, however, Miyatake et al. (2022) showed that the use of only a 'broadly accurate' standard halo model leads to significant offsets in the recovered cosmological parameters from a 2×2pt analysis of HOD-populated numerical simulations.Consistency studies between the observed small-scale clustering and galaxygalaxy lensing signals cast similar doubt on the accuracy of any standard halo model analysis (Leauthaud et al. 2017;Lange et al. 2021;Amon et al. 2022b).
Arguably the two most flawed approximations in the standard halo model formalism are that (i) haloes, and therefore galaxies, can overlap, and (ii) that haloes trace the matter distribution with a linear and scale independent halo bias.Previous attempts to improve these approximations have used halo exclusion formulations (Giocoli et al. 2010) to solve the first problem, combined with radial bias functions which add scale dependence to the halo bias model (Tinker et al. 2005;van den Bosch et al. 2013;Cacciato et al. 2013).
In this analysis we instead follow the method proposed by Mead & Verde (2021), accounting for both scale-dependent nonlinear halo bias and halo exclusion by incorporating the halo bias measured directly from the DarkEmulator suite of cosmological simulations (Nishimichi et al. 2019).As shown by Mahony et al. (2022), this necessary upgrade to the standard halo model leads to sufficient accuracy in the recovered cosmological parameters from a 2×2pt+SMF analysis, for the statistical power of current imaging surveys.
Other approximations in a halo model analysis include that the halo mass is the sole variable that determines the properties of the haloes and their occupying galaxies.Galaxy properties and the clustering of haloes are, however, expected to have a secondary dependence, on their local environment and assembly history (see Wechsler & Tinker 2018, and references therein).Furthermore the adopted halo density profile is modelled from dark matter-only numerical simulations, even though hydrodynamical simulations show that these profiles are modified by the presence of active galactic nuclei (Schaller et al. 2015;Wang et al. 2020).In Debackere et al. (2021), it is shown that to account for baryon physics, it is sufficient to leave the concentration of dark matter haloes free.Amon et al. (2022b) review the literature studies on the impact of these two approximations on 2×2pt halo model studies with luminous red galaxies.Motivated by their conclusions, we choose to adopt nuisance parameters in our halo model to encapsulate the uncertainty of the impact of assembly bias and baryon feedback within our error budget.
In this paper, we analyse the most recent data release from the Kilo-Degree Survey (KiDS-1000, Kuijken et al. 2019;Giblin et al. 2021;Hildebrandt et al. 2021), spanning over 1000 square degrees of imaging in nine bands from the optical through to the near-infrared.Our main 'KiDS-Bright' galaxy sample (Bilicki et al. 2021) benefits from the 180 square degree overlap between KiDS and the spectroscopic Galaxy And Mass Assembly (GAMA) survey (Driver et al. 2011).As an essentially complete spectroscopic survey to r < 19.8, GAMA serves as an extensive training set for machine learning and the calibration of different sample selections.The resulting GAMA-trained photometric redshifts and stellar mass estimates for the KiDS-Bright sample have an enhanced accuracy and precision that benefits this galaxy-galaxy lensing and galaxy clustering study.In order to simultaneously constrain cosmology and galaxy bias, we use the 2×2pt+SMF combination of galaxy clustering and galaxygalaxy lensing, as well as constraints on galaxy abundances in the form of the stellar mass function.
We improve previous related 2×2pt+SMF studies by (i) using a more accurate analytical model with the addition of nonlinear halo bias (Mead & Verde 2021), taking into account the halo exclusion and scale dependence, (ii) taking the latest lensing and clustering data from a single survey, and (iii) using the full analytical covariance matrix including the cross-variance between all observables.Our analysis is highly complementary to the emulator based 2×2pt halo model analysis of the Hyper Suprime Camera Survey (HSC, Miyatake et al. 2021).
Throughout this paper, all radii and densities are given in comoving units, 'log' is used to refer to the 10-based logarithm, and 'ln' for natural logarithm.All the quantities that depend on the Hubble parameter adopt units of h, where h = H 0 /100 km s −1 Mpc −1 .We also use ρ m as the present day mean matter density of the Universe, ρ m = Ω m,0 ρ crit , where ρ crit = 3H 2 0 /(8πG) and the halo masses are defined as M = 4πr 3 ∆ ∆ ρ m /3 enclosed by the radius r ∆ within which the mean density of the halo is ∆ times ρ m , with ∆ = 200.
This paper is organised as follows.In Sec. 2, we review our analytical model to compute the galaxy stellar mass function, the galaxy-galaxy correlation function, and the galaxy-galaxy lensing signal using the halo model combined with a model that describes halo occupation statistics as a function of galaxy stellar mass.In Sec. 3, we introduce the 2×2pt+SMF KiDS measure-ments, specifics of the covariance calculation and our Bayesian analysis methodology.Our main results are presented in Sec. 4, concluding our findings in Sec. 5.

The Halo model
The halo model is an analytic framework that can be used to describe the clustering of matter and its evolution in the Universe (Seljak 2000;Peacock & Smith 2000;Cooray & Sheth 2002;van den Bosch et al. 2013;Mead et al. 2015).It is built upon the statistical description of the properties of dark matter haloes (namely the average density profile, large scale bias and abundance) as well as on the statistical description of the galaxies residing in them, using halo occupation distributions (HOD).The model is sufficiently flexible to consistently describe the statistical weak lensing signal around a selection of galaxies, their clustering, abundances and cosmic shear signal.

Halo model ingredients
We assume that dark matter haloes are spherically symmetric on average, and have density profiles, ρ(r|M) = M u h (r|M), that depend only on their mass M, and u h (r|M) is the normalised density profile of a dark matter halo.Similarly, we assume that satellite galaxies in haloes of mass M follow a spherical number density distribution n s (r|M) = N s u s (r|M), where u s (r|M) is the normalised density profile of satellite galaxies.All central galaxies are positioned at the centre of their halo: r = 0. We assume that the density profile of dark matter haloes follows a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997).Since centrals and satellites are distributed differently, we write the galaxygalaxy power spectrum, P gg (k), as a combination of the central 'c', satellite 's' and cross power spectrum, with and the galaxy-matter power spectrum, P gm (k), Here A c = n c /n g and A s = n s /n g = 1 − A c are the central and satellite fractions, respectively, and the average number densities n g , n c and n s follow from: where 'x' stands for 'g' (for galaxies), 'c' (for centrals) or 's' (for satellites), ⟨N x |M⟩ is the average number of galaxies given halo mass M, and n(M) is the halo mass function in the following form: with ν = δ c /σ(M) being the peak height.Here δ c is the critical overdensity required for spherical collapse at redshift z, and σ(M) is the mass variance.For f (ν) we use the fitting function to the numerical simulations presented in Tinker et al. (2010).In addition, it is common practice to split two-point statistics into a 1-halo term (both points are located in the same halo) and a 2-halo term (the two points are located in different haloes).The 1-halo terms are: and all other terms are given by: Here 'x' and 'y' are either 'c' (central), 's' (satellite), or 'm' (matter), P is a Poisson parameter that captures the scatter in the number of satellite galaxies at fixed halo mass (in this case a free parameter -we define the P in detail using equations 24 and 25) and we have defined the mass, central and satellite profiles as and with ũh (k|M) and ũs (k|M) the Fourier transforms of the halo density profile and the satellite number density profile, respectively, both normalised to unity [ũ(k=0|M)=1].The various 2-halo terms are given by: where P lin (k) is the linear power spectrum, obtained using the Eisenstein & Hu (1998) β NL is measured using the DarkQuest emulator (Nishimichi et al. 2019;Miyatake et al. 2022;Mahony et al. 2022), by measuring the non-linear halo-halo power spectrum and then dividing it by the linear matter power spectrum multiplied with the product of linear bias factors (Mead & Verde 2021, Eq. 23).Due to the definition of β NL , this measurement also holds true for galaxy-galaxy and galaxy-matter correlations.As shown in Mahony et al. (2022), this function is cosmology dependent, but does not account for assembly bias effects.In this paper, owing to the volume-limited mix of all types of galaxies used in our analysis, we consider any assembly bias to be a subdominant effect as the secondary properties are unlikely to manifest for a non-specific galaxy type selection (Wechsler & Tinker 2018).
Numerically, the integrals in the halo model are not integrated from zero to infinity, but rather between a wide range of halo masses.Special care has to be taken to account for the masses outside of the integration limits, for which an appropriate correction is applied (as derived in Cacciato et al. (2009), Eqs. 24 and25, andin Mead et al. (2020); Mead & Verde (2021), Appendices A in both papers).The two-point correlation functions corresponding to these power-spectra are obtained by Fourier transformation: For the halo bias function, b h , we use the fitting function from Tinker et al. (2010), as it was obtained using the same numerical simulation from which the halo mass function was obtained.We have adopted the parametrisation of the concentration-mass relation, given by Duffy et al. (2008): We allow for an additional normalisation f h,s , such that where f h is the normalisation of the concentration-mass relation for dark matter haloes ũh (k|M), and f s is the normalisation of the concentration-mass relation for the distribution of satellite galaxies ũs (k|M).The profiles ũh (k|M) and ũs (k|M) are both assumed to be non-truncated NFW profiles, with the same virial mass.
Our adoption here of separate concentration-mass relations for dark matter haloes and satellite galaxies provides enough flexibility in the model to capture the uncertain impact of baryon feedback (for the scales adopted, Debackere et al. 2020Debackere et al. , 2021;;Amon et al. 2022b), and it has been used in the literature (Cacciato et al. 2013;Viola et al. 2015;van Uitert et al. 2016;Dvornik et al. 2018) to account for such effects.This additional flexibility is motivated by the fact that in the simulations, the AGN feedback pushes the baryons and dark matter from halo centres towards outskirts, and by that effectively changing the concentration of the matter distribution (Debackere et al. 2020;Mead et al. 2020).This is also supported by observations (Viola et al. 2015), which showed that the preferred value for the concentration normalisation is lower than 1.Using the halo model with these extra parameters is a benefit over the emulators that are based on darkmatter only simulations (as for instance the DarkQuest emulator Nishimichi et al. 2019;Miyatake et al. 2022;Mahony et al. 2022), since they do not offer a simple way to accommodate for such flexibility, nor require simulations (Schneider & Teyssier 2015).
In the halo model we do not consider the mis-centred central term, as for a selection of galaxies the signature is accounted for through the terms for satellite galaxies, which do not reside in the centres of haloes by definition.What is more, the satellite galaxies populate haloes regardless of the existence of a central galaxy, which further removes the need for a mis-centred term (no central condition is enforced).

Conditional stellar mass function
We model the galaxy stellar mass function and halo occupation statistics using the Conditional Stellar Mass Function (CSMF, motivated by Yang et al. 2008;Cacciato et al. 2009Cacciato et al. , 2013;;Wang et al. 2013;van Uitert et al. 2016).The CSMF, Φ(M ⋆ |M), specifies the average number of galaxies of stellar mass M ⋆ that reside in a halo of mass M. In this formalism, the halo occupation statistics of central galaxies are defined via the function: In particular, the CSMF of central galaxies is modelled as a lognormal, and the satellite term as a modified Schechter function, where σ c is the scatter between stellar mass and halo mass and α s governs the power law behaviour of satellite galaxies.Note that M * c , σ c , ϕ * s , α s and M * s are, in principle, all functions of halo mass M, but here we assume that σ c and α s are independent of the halo mass M. Inspired by Yang et al. (2008), who studied the halo occupation properties of galaxies in the Sloan Digital Sky Survey, we parametrise M * c , M * s and ϕ * s as: and where m 13 = M/(10 13 M ⊙ h −1 ).In their analysis of the stellar-tohalo mass relation of GAMA galaxies van Uitert et al. (2016) find that varying the pre-factor of 0.56 in Equation 20does not significantly affect the results, therefore we retain this normalisation in our analysis.We can see that the stellar to halo mass relation for M ≪ M 1 behaves as M * c ∝ M γ 1 and for M ≫ M 1 , M * c ∝ M γ 2 , where M 1 is a characteristic mass scale and M 0 is a normalisation.Here γ 1 , γ 2 , b 0 and b 1 are all free parameters that govern the two slopes of the stellar-to-halo mass relation and the normalisation of the Schechter function.The choice of functional form of the CSMF is motivated by the good performance as seen in previous lensing and combined lensing and clustering studies.In Eq. 19, we adopt an effective stellar-to-halo mass relation for our mixed-population of red and blue galaxies.Bilicki et al. (2021) demonstrate a strong colour-dependence to this relationship, and future studies will investigate including a red/blue galaxy split in our analysis, which can also help to improve the modelling of intrinsic galaxy alignments (e.g.Li et al. 2021).
From the CSMF it is straightforward to compute the galaxy stellar mass function (SMF) and the halo occupation numbers.The galaxy stellar mass function is in this case given by and the average number of galaxies with stellar masses in the range M ⋆,1 ≤ M ⋆ ≤ M ⋆,2 is given by: where 'x' is either 'c' (central), 's' (satellite), or the total contribution from all galaxies.In order to predict the satellitesatellite term for the galaxy clustering power spectra (Eq.6), we use that where P(M) is the mass-dependent Poisson parameter defined as: which is unity if ⟨N s |M⟩ is given by a Poisson distribution, larger than unity if the distribution is wider than a Poisson distribution (also called super-Poissonian distribution) or smaller than unity if the distribution is narrower than a Poisson distribution (also called sub-Poissonian distribution).In our fiducial analysis we limit ourselves to cases in which P(M) is independent of halo mass, i.e., P(M) = P, and we treat P as a free parameter.In Appendix 4.4 we present an extension to our fiducial analysis, allowing for mass-dependence in the Poisson parameter, based on the observed distribution of satellite galaxies in the GAMA group catalogue (Robotham et al. 2011).Our findings are sensitive to the selection criteria chosen for the GAMA group catalogue.We are nevertheless able to conclude that assuming the Poisson parameter is independent of halo mass impacts our primary cosmological parameter constraints (mostly S 8 ) at an acceptable ∼ 1σ level.
Overall, all the free parameters used to describe the halo occupation distributions and the connection with the dark matter are: Priors on these parameters are broad, assuming wide uniform distributions, similar to the priors used in two studies of the galaxy-halo connection that both used GAMA and KiDS data (van Uitert et al. 2016;Dvornik et al. 2018).The halo occupation distribution parameters could in principle also depend on redshift and halo mass, but furthering the complexity of the model, by increasing the number of parameters, would not be justified by the data.Our parameters describe an effective model over the redshift range in the analysis.Priors and their ranges can be found in Table 2.

Projected lensing and clustering functions
Once P gg (k) and P gm (k) have been determined, it is fairly straightforward to compute the projected galaxy-galaxy correlation function, w p (r p ), and the excess surface density (ESD) profile, ∆Σ(r p ).The projected galaxy-galaxy correlation function, w p (r p ), is related to the real-space galaxy-galaxy correlation function, ξ gg (r), according to Here ξ gg (r p , r π , z) is the redshift-space galaxy-galaxy correlation function, r π is the redshift-space separation perpendicular to the line-of-sight and r π,max is the maximum integration range used for the data (here we use is the separation between the galaxies, L l (x) is the l th Legendre polynomial, and ξ 0 , ξ 2 , and ξ 4 are given by where and with a = 1/(1 + z) the scale factor, D(z) the linear growth factor, and the mean bias of the galaxies in consideration.Eq. 27 accounts for the large-scale redshift-space distortions due to infall (the 'Kaiser'-effect), which is necessary because the measurements for w p (r p ) are obtained for a finite r max .We note that whilst this Kaiser (1987) formalism is only strictly valid in the linear regime, we adopt the non-linear galaxy-galaxy correlation function, ξ gg (r), in Eqs.28 -30, with the non-linearities captured through the halo model power spectra in Eq. 13. van den Bosch et al. (2013) show that this modification provides a more accurate correction for the residual redshift space distortions, and that ignoring the presence of residual redshift space distortions leads to systematic errors that can easily exceed 20 percent on scales with r p > 10h −1 Mpc (Cacciato et al. 2013).
The excess surface density profile (ESD), ∆Σ(r p ), is defined as Here Σ(r p ) is the projected surface mass density, which is related to the galaxy-dark matter cross correlation, ξ gm (r), according to The final model predictions and the covariance matrix are binaveraged to the bin widths of the data vectors.

Cosmological parameters
The cosmological parameters in our model are described by the vector: (36) As mentioned in Sect. 1, the goal of this paper is to use the ESD, w p and SMF data to constrain σ 8 and Ω m .Because of that, we set the priors for those two parameters to be uninformative and set their ranges following the latest KiDS cosmic shear analysis (Asgari et al. 2021b).The last 3 cosmological parameters are shown to be poorly constrained using the ESD, w p , and SMF data (Cacciato et al. 2013;Mandelbaum et al. 2013), thus they form a set of secondary cosmological parameters with informative priors.Priors and their ranges can be found in Table 2 2 .In Appendix B we verify that our choice of priors do not inform the main cosmological parameters.
2 Our prior range is larger than the range of available nodes in the DarkQuest emulator.Owing due to the iterative updates to the β NL es-

Data and sample selection
In this analysis we combined three observables from the Kilo-Degree Survey (KiDS): galaxy abundances in the form of the galaxy stellar mass function, galaxy clustering in the form of the projected galaxy correlation function, and galaxy-galaxy lensing in the form of excess surface density profiles.Our KiDS observations were taken with OmegaCAM (Kuijken 2011), a 268million pixel CCD mosaic camera mounted on the VLT Survey Telescope.These instruments were designed to perform weak lensing measurements, with the camera and telescope combination providing a fairly uniform point spread function across the field-of-view (de Jong et al. 2013).
We analysed the latest data release of the KiDS survey (KiDS-1000, Kuijken et al. 2019), containing observations from 1006 square-degree survey tiles.Specifics of the survey, the calibration of the source shapes and photometric redshifts are described in Kuijken et al. (2019); Giblin et al. (2021); Hildebrandt et al. (2021), respectively.The companion VISTA-VIKING (Edge et al. 2013) survey has provided complementary imaging in near-infrared bands (ZYJHK s ), resulting in a unique deep, wide, nine-band imaging dataset (Wright et al. 2019).The default photo-z estimates provided as part of the KiDS survey were derived with the Bayesian Photometric Redshift approach (BPZ, Benitez 2000).
We used shape measurements based on the r-band images, which have an average seeing of 0.66 arcsec.The galaxy shapes were measured using lensfit (Miller et al. 2013), which has been calibrated using image simulations described in Kannawadi et al. (2019).This provides galaxy ellipticities (ϵ 1 , ϵ 2 ) with respect to an equatorial coordinate system, and an optimal weight.
The galaxies used for our lens and clustering sample were taken from the 'KiDS-Bright' sample (Bilicki et al. 2021).This sample mimics the selection of GAMA galaxies (Driver et al. 2011), by applying the condition m r < 20.0.For these galaxies a different method of determining the photometric redshifts was employed using the ANNz2 (Artificial Neural Network) machine learning method (Sadeh et al. 2016), with the spectroscopic GAMA survey, which is 98.5% complete to r < 19.8, as a training set (Bilicki et al. 2018(Bilicki et al. , 2021)).Comparing the obtained redshifts with the spectroscopic redshifts from the matched galaxies between KiDS-Bright and GAMA, Bilicki et al. (2021) concluded that the ANNz2 photo-z are highly accurate with a mean offset of δ z = 5 × 10 −4 , and a scaled mean absolute deviation scatter of σ z = 0.018(1 + z).
Stellar mass estimates for the KiDS-Bright sample are obtained using the LePhare template fitting code (Arnouts et al. 1999;Ilbert et al. 2006).In these fits, ANNz2 photo-z estimates are used as input redshifts for each source, treating them as if they were exact, neglecting the percent error associated with the ANNz2 redshift.In practice, this error has little impact on the fidelity of the stellar mass estimates (Taylor et al. 2011).The estimates assume a Chabrier (2003) initial mass function, the Calzetti et al. (1994) dust-extinction law, Bruzual & Charlot (2003) stellar population synthesis models, and exponentially declining star formation histories.The input photometry to LePhare is extinction corrected using the Schlegel et al. (1998) maps with the Schlafly & Finkbeiner (2011) coefficients, as described in Kuijken et al. (2019).Bilicki et al. (2021) found that the KiDS-Bright stellar mass estimates are in excellent agreement with independent stellar mass estimates from Wright et al. (2016) that combine GAMA timation and the quick convergence we find towards parameters within the emulator's range, this does not pose an issue.2021) estimated the overall systematic uncertainty on the stellar mass estimates of the KiDS-Bright sample, combining the uncertainty arising from the LePhare model fit, the photometric redshift scatter, and the difference found when exchanging elliptical aperture magnitudes for Sérsic model magnitudes.They estimated an overall uncertainty of σ M * = 0.12 dex for the KiDS-Bright sample.This systematic uncertainty also includes the estimated Eddington systematic bias of ∼ 0.027 dex (Brouwer et al. 2021), which is estimated from the population of red and blue galaxies and it is considered a worst-case scenario.We choose to account for both statistical and systematic uncertainty in the stellar mass estimates through the nuisance parameter σ c , in equation 17, which provides the freedom to model both the intrinsic and measurement noise scatter in the stellar-to-halo mass relation (Leauthaud et al. 2012;Bilicki et al. 2021).Furthermore, as the systematic and statistical uncertainties are comparable in power, the entries in the SMF and cross-covariances are inflated by a factor of 2 to account for the uncertainty arising from Eddington bias and the systematic shift in stellar masses, and not only through the σ c parameter.Due to the weak cosmology dependence of the SMF, this primarily increases only the uncertainty of our HOD parameters, as the SMF is in the first place used to break degeneracies in our HOD part of the halo model.

Stellar mass function measurements and sample selection: SMF
Our SMF measurements are performed using the maximumvolume weighting method (Schmidt 1968;Saunders et al. 1990;Cole 2011;Baldry et al. 2012;Wright et al. 2017).We weight each galaxy i by the inverse of the comoving volume over which the galaxy would be visible, given the magnitude limit of the whole sample, 1/V max,i .To estimate the number density Φ(M ⋆ ), we have to derive M ⋆,lim (z), the completeness in stellar mass as a function of redshift for our flux-limited sample.For the 1/V max technique, we need to know z max,i , the maximum redshift beyond which galaxy i with stellar mass M ⋆,i would no longer be part of the subsample (Weigel et al. 2016).This is done by determining the point at which the sample begins to become incomplete.Usually this process contains a potentially biased visual inspection.To avoid any bias, we instead adopt the automated method presented by Wright et al. (2017), using the MassFunc-FitR package.The algorithm estimates the turn-over point of the number density distribution in bins of comoving distance and stellar mass independently.In each fine bin of comoving distance, we take the mass at the peak density as the mass turn-over point.In each fine bin of stellar mass, we take the largest comoving distance at median stellar mass density as the distance turn-over point.The obtained turn-over points are then fit with a high-degree polynomial resulting in a smooth form for the stellar mass limit as a function of redshift.This limit can be compared to the M * − z ANNz2 distribution of the full KiDS-Bright galaxies in Fig. 1.
In Fig. 2 we present the stellar mass function of the volumelimited KiDS-Bright sample from the galaxies in the 6 stellar mass bins, Φ(M ⋆ ), which has a median redshift of z = 0.25.This is determined from the galaxy counts within the stellar mass limit, with errors derived analytically in Appendix A. We find good agreement between the KiDS measurement and the stellar mass function from Wright et al. (2018), evaluated at the median redshift of our sample.Wright et al. (2018)   The full sample is shown with a logarithmic hexagonal density plot.The blue line shows the stellar mass limit determined using the automated method presented by Wright et al. (2017).Red boxes show the six stellar mass bins used in the analysis, with individual galaxies plotted as black dots.The bin ranges were chosen in such a way as to achieve a good signal-to-noise ratio in all bins for our galaxygalaxy lensing and galaxy clustering measurements.analysis using spectroscopic data from GAMA, COSMOS and HST.This comparison therefore demonstrates the agreement in the SMF between spectroscopic data and our photometric KiDS-Bright sample of galaxies, demonstrating that our stellar mass estimates are robust to the uncertainty in the photometric redshifts (Taylor et al. 2011;Bilicki et al. 2021;Brouwer et al. 2021).
As galaxy bias is inherently dependent on the stellar mass of the galaxy (Dvornik et al. 2018), we analyse the weak lensing and galaxy clustering of the KiDS-Bright galaxies grouped into 6 stellar mass bins.We choose to limit our analysis to galaxies within the stellar mass range of 9.1 < log(M ⋆ /h −2 M ⊙ ) ≤ 11.3, with the number of bins, and bin limits chosen in such a way to achieve a similar and significant signal-to-noise ratio in all bins.Using the redshift-dependent stellar mass limit, we define upper redshift bounds to ensure each stellar mass bin is volumelimited, as indicated with red boxes in Fig. 1.The lower redshift bound is set to contain 95 percent of the volume-limited sample.The number of galaxies, median stellar mass and redshift of each bin is reported in Table 1.

Galaxy-galaxy lensing measurement: ESD
As shown in Bilicki et al. (2021), the excellent ANNz2 photometric redshift estimates for the galaxies in the KiDS-Bright sample allow for robust estimates of their physical characteristics, in particular the stellar mass.In this section we combine this information with accurate shape measurements for more distant KiDS sources from Giblin et al. (2021) to measure the galaxy-  2018), evaluated at the median redshift of our sample (black dashed).The blue line and shaded region indicates the best fit model and 68% confidence levels of our best fit halo model (Eq. 22).We caution that the quality of the fit cannot be judged by eye, because of the covariance in the data, between the data points and between the other observables.The reduced χ 2 value for this observable is 1.05 (DoF = 14.58, p-value = 0.39), estimated using the method presented in Appendix C. Lower panel: The fractional errors on the data and the model, ∆Φ/δΦ.galaxy lensing signal.To quantify the weak gravitational lensing signal we use source galaxies from KiDS DR4 with a BPZ photo-z in the range 0.1 The lensing signal of an individual lens is too small to be detected, and hence we compute a weighted average of the tangential ellipticity ϵ t as a function of projected distance r p using a large number of lens-source pairs.In the weak lensing regime this provides an unbiased estimate of the tangential shear, γ t , which in turn can be related to the excess surface density (ESD) ∆Σ(r p ), defined as the difference between the mean projected surface mass density inside a projected radius r p and the mean surface density at r p (as in Eq. 34, for more details see appendix C from Dvornik et al. 2018).
We compute a weighted average to account for the variation in the precision of the shear estimate, captured by the lensfit weight w s (see Fenech Conti et al. 2017;Kannawadi et al. 2019, for details), and the fact that the amplitude of the lensing signal depends on the source redshift.The weight assigned to each Table 1: KiDS-Bright stellar mass samples: overview of the number of galaxies/lenses, median stellar masses M ⋆,med and median redshifts z med .Stellar masses are given in units of log(M ⋆ /h −2 M ⊙ ).The final sample of galaxies used is a small subsample of all KiDS-Bright galaxies (∼ 1 million).
the product of the lensfit weight w s and the square of Σ −1 cr,ls -the effective inverse critical surface mass density, which is a geometric term that downweights lens-source pairs that are close in redshift (e.g.Bartelmann & Schneider 2001).
We compute the effective inverse critical surface mass density for each lens using the photo-z of the lens z l and the full normalised redshift probability density of the sources, n(z s ).The latter is calculated employing the self-organising map calibration method (Wright et al. 2020) as applied to KiDS DR4 in Hildebrandt et al. (2021).The resulting effective inverse critical surface density can be written as: where D(z l ), D(z s ), D(z l , z s ) are the angular diameter distances to the lens, source, and between the lens and the source, respectively.For the lens redshifts z l we use the ANNz2 photo-z of the KiDS-Bright foreground galaxy sample.We implement the contribution of z l by integrating over the redshift probability distributions p(z l ) of each lens.The lensing kernel is wide and therefore the resulting ESD signals are not sensitive to the small wings of the lens redshift probability distributions.We can thus safely approximate p(z l ) as a normal distribution centred at the lenses photo-z, with a standard deviation σ z /(1 + z l ) = 0.018 (Bilicki et al. 2021).From previous KiDS GGL studies we know that the error on the mean and width of source n(z s ) are not biasing the galaxy-galaxy lensing signal (as shown in Dvornik et al. 2017).
For the source redshifts z s we follow the method used in Dvornik et al. (2018), by integrating over the part of the redshift probability distribution n(z s ) where z s > z l .The galaxy source sample is specific to each lens redshift z l , with a minimum photometric redshift z s = z l + δ z , with δ z = 0.2 that is used to remove sources that are physically associated with the lenses.Thus, the ESD can be directly computed in bins of projected distance r p to the lenses as: where Σ ′ cr,ls ≡ 1/ Σ −1 cr,ls , the sum is over all source-lens pairs in the distance bin, and is an average correction to the ESD profile that has to be applied to account for the multiplicative bias m in the lensfit shear estimates.The sum goes over thin redshift slices for which m i is obtained using image simulations (Kannawadi et al. 2019), weighted by w ′ = w s D(z l , z s )/D(z s ) for a given lens-source sample.The value of m is −0.003 for the 6 stellar mass bins, independent of the scale at which it is computed.The uncertainty in m is not marginalised over, as the contribution of the central m value is at most a percent of the total error budget of the galaxy-galaxy lensing signal.We note that the measurements presented here are not corrected for the contamination of the source sample by galaxies that are physically associated with the lenses (the so-called 'boost correction').The impact on ∆Σ is minimal, because of the weighting with the inverse square of the critical surface density in Eq. ( 38), (see for instance the bottom panel of fig.A4 in Dvornik et al. 2017) and the removal of the sources physically associated with the lens from our signal measurements.The effect of using photometric lenses in the ESD measurements is directly accounted for in our estimator and the covariance matrix.We subtract the signal around random points, which suppresses any large-scale systematics and sample variance (Singh et al. 2017).This empirical 'random' correction for large-scale sample variance has been shown to improve robustness on the measurement scales which are particularly relevant to constrain linear bias (Dvornik et al. 2018).We find the random correction for the KiDS-Bright sample becomes significant at scales R ≳ 3h −1 Mpc, rising to more than 100% of the ESD signal in the three lowest stellar mass bins, and it thus dictates the range of measurement scales we use in the analysis.On these large scales the random correction is more than four times larger than the statistical uncertainty (see Appendix D for details).The resulting random-corrected galaxy-galaxy lensing ESD measurements for the six stellar mass bins are shown in Fig. 3 Here we count the number of galaxy-galaxy (DD), randomrandom (RR) and galaxy-random (DR) pairs, as a function of the pair's transverse r p and radial r π comoving separation.The accuracy of galaxy clustering measurements with this estimator depends critically on the quality of the random, R, catalogues.We use the Johnston et al. (2021b) organised random methodology that has been shown to recover unbiased clustering measurements in a series of mock galaxy catalogue analyses for the KiDS-Bright sample.Using machine learning, we infer the high-dimensional mapping between the observed onsky galaxy number density and three systematic-tracer variables; atmospheric seeing, point spread function ellipticity and limiting magnitude.Systematically-induced density variations across the survey footprint can then be defined.We randomly distribute clones of the real galaxies across the survey footprint, preserving the on-sky systematic density patterns, and matching the on-sky systematic-tracer properties to that of the clone's parent galaxy.By retaining the photometric properties of the parent for each clone, selection effects are accurately mirrored in the organised randoms for any galaxy sub-sample, for example the 6 different Fig. 3: Galaxy-galaxy lensing: The stacked ESD profiles of the six stellar mass bins in the KiDS-Bright galaxy sample defined in Table 1.The solid lines represent the best-fitting fiducial ESD halo model (Sec.2.3, Eq. 34) as obtained using an MCMC fit, with the 68 percent confidence interval indicated with a shaded region.We caution that the quality of the fit cannot be judged by eye, because of the covariance in the data between the observed bins and also between the observables.The reduced χ 2 value for this observable is 1.28 (DoF = 73.18,p-value = 0.05), estimated using the method presented in Appendix C. stellar mass bins in our analysis.We use 20 times more randoms than data points as presented by Johnston et al. (2021b).
The projected clustering correlation function is estimated through an integral over the line-of-sight separation, limited by a maximum defined distance r π,max , ŵp (r p ) = r π,max −r π,max ξgg (r p , r π ) dr π . (42) When analysing spectroscopic data, this continuous integral is estimated using a discrete sum, typically adopting uniform bins in r π , with r π,max ranging from 40h −1 Mpc to 100h −1 Mpc (as in for instance Mandelbaum et al. 2010;Farrow et al. 2015).Here the r π,max limits are chosen to maximise the number of correlated galaxy pairs along the line-of-sight in the presence of redshift space distortions, whilst minimising the noise arising from the inclusion of uncorrelated objects.With our KiDS-Bright photometric sample we have an additional uncertainty in the true redshift, σ z = 0.018(1 + z), which translates into an uncertainty on the radial distance of the order ∼ 100h −1 Mpc, This renders the approach taken for spectroscopic samples sub-optimal in terms of signal-to-noise.We therefore choose to follow the approach of Johnston et al. (2021a) who optimised the projected galaxy clustering analysis of the photometric Physics of the Accelerating Universe Survey (PAUS), using dynamic binning in r π out to a maximum r π = 233h −1 Mpc.This is motivated by the fact that PAUS photometric redshifts show a similar uncertainty as the KiDS-Bright sample.Using a mock galaxy catalogue, Johnston et al. (2021a) demonstrated that by allowing for an increase in the bin size from small to large values of r π , their approach maximises the count of physically associated objects, whilst minimising noise at large-r π with the broader bin size.Given the similar photometric redshift properties of KiDS-Bright and PAUS, we adopt their 12-r π -bin adapted Fibonacci sequence in our estimator.Johnston et al. (2021a) analysed mock GAMA galaxy catalogues with PAU-like photometric redshifts to compare the projected clustering correlation function estimator ŵp (r p ) with the measurements using spectroscopic redshifts.Adopting dynamic binning and random galaxy catalogues that mimic both the position and photometric redshift uncertainty of the real galaxy sample, they found a roughly scale-independent bias with ŵp /w p spec ≃ 0.8.As such the dynamic binning and organised randoms only partially correct the correlation functions for the dilution introduced by photometric redshift uncertainty.Future work will focus on accounting for this dilution effect accurately in the theoretical prediction.For the purposes of this analysis, however, we choose to include a free dilution parameter D, which is used to correct the galaxy clustering measurements in the following way: ŵp,corr (r p ) = [1 + D] ŵp (r p ) . (43) We adopt a uniform prior for D with the range between 0 and 0.3 and use a single parameter to scale all six stellar mass bins.This prior was motivated by a series of mock KiDS-Bright galaxy clustering analysis using MICE2 (Fosalba et al. 2015a,b;Crocce et al. 2015;Carretero et al. 2015;Hoffmann et al. 2015), where we confirmed the findings of Johnston et al. (2021a) and found no strong dependence of the dilution effect on stellar mass.We note that a similar correction was applied to the Dark Energy Survey (DES) photometric clustering measurements (Pandey et al. 2022;DES Collaboration et al. 2022, referred therein as X lens ).The prior and motivation behind the introduction of their systematic nuisance parameter differs, however.The resulting projected clustering measurements for the six stellar mass bins are shown in Fig. 4.

Accounting for the cosmology dependence of distance measures
To obtain estimates of the stellar mass function (SMF, Sect. 3.1, Fig. 2), the galaxy-galaxy lensing (ESD , Sect.3.2, Fig. 3), and the projected galaxy clustering (w p , Sect.3.3, Fig. 4), we have adopted a fiducial flat ΛCDM cosmology with Ω m = 0.3 to compute distances.As such our 2×2pt+SMF data vector is cosmology-dependent, with changes in the fiducial cosmology changing the distance-redshift relation, which in turn shifts Fig. 4: Galaxy clustering: The projected galaxy clustering signal of the six stellar mass bins in the KiDS-Bright galaxy sample defined in Table 1.The solid lines represent the best-fitting fiducial halo model (Sec.2.3, Eq. 27) as obtained using an MCMC fit, with the 68 percent confidence interval indicated with a shaded region.We caution that the quality of the fit cannot be judged by eye, because of the covariance in the data, between the observed bins and between the observables.The reduced χ 2 value for this observable is 1.42 (DoF = 71.62,p-value = 0.01), estimated using the method presented in Appendix C. galaxies between the stellar mass bins and lens-source pairs between the radial separation bins.
At the mean redshift of the KiDS-Bright sample, the effect of changing Ω m within our prior limits introduces changes in distance estimates at the level of few percent.The approximation that the measurements are effectively independent of cosmological parameters within their observational uncertainties (Mandelbaum et al. 2013;Cacciato et al. 2013) no longer holds for surveys with a statistical power that is similar or better than KiDS.
In this analysis we account for the cosmology dependence of our data vector following the correction procedure presented in More (2013) and More et al. (2015), which modifies the model prediction for each cosmology targeted by the likelihood sampler.
First we define a cosmology-dependent comoving separation r model p for our target model, relative to the comoving separation r p that has been used to calculate our data vectors at a fixed fiducial cosmological model, Here χ is the comoving distance to the median lens redshift z med in our target cosmological model C model , or in our fiducial cosmological model C fid .The galaxy clustering prediction for our target model is then given by where E(z) is the Hubble parameter.The galaxy-galaxy lensing prediction for our target model is given by where Σ cr is the critical surface density calculated for the median redshift of the lenses z med and a fixed source redshift z s = 0.6.Note that calculating the more precise estimate for Σ cr using Eq.38 is not necessary in this instance, as Σ cr only has a weak cosmology dependence.Finally the predictions of abundances of galaxies in the target cosmology is given by which is implicitly correcting the surveyed volume in the stellar mass function calculation.Here the z l and z u are the lower and upper redshift limits in our samples.

Covariance matrix
The covariance matrix used in this analysis is based on the analytical approach detailed in Dvornik et al. (2018) and Joachimi et al. (2021), with the addition of the analytical covariance matrix for the SMF and the cross terms between the SMF and 2-point correlation functions.The new terms for the SMF covariance and the cross covariance between the SMF and 2-point functions are presented in Appendix A. Our implementation of the analytical covariance derivation was validated against theory (Pielorz et (2020) and Joachimi et al. (2021).Survey area effects on the variance were calculated using the accurate, survey dependent and data based Healpix method presented in Joachimi et al. (2021), equation E.10.

Likelihood and iterative updates
We use Bayesian inference to determine the posterior probability distribution P(θ | d) of the model parameters θ, given the data d.
According to Bayes' theorem, P(θ | d) is: where P(d | θ) is the likelihood of the data given the model parameters, P(θ) is the prior probability of these parameters, and is the evidence for the model.Since, we do not perform model selection in this analysis, the evidence just acts as a normalisation constant which we do not need to calculate.Given this, the likelihood distribution P(d | θ) is assumed to be Gaussian: where C is the full covariance matrix for all the observables, containing their auto-and cross-correlations, |C| its determinant, m(θ) the model given the parameters θ, and n the number of observable bins.Priors can be found in Table 2.For the Bayesian inference we use the MCMC sampler emcee (Foreman-Mackey et al. 2013).
The posterior distribution in such highly multi-dimensional parameter spaces has numerous degeneracies and can be very difficult to sample from.Thus the choice of proposal distributions is very important in order to achieve fast convergence and reasonable acceptance fractions for the proposed walker positions.To do so, we combine the default stretch move in the emcee with the proposal function based on the kernel density estimator of the complementary ensemble of walkers (Foreman-Mackey et al. 2013) 3 in such way that at every step of the sampler run, there is a 50% chance to use one of the proposal methods.This setup has one downside, and that is that it uses many walkers, and thus computing power.On the other hand, the convergence is faster and the resulting auto-correlation times are shorter, giving us shorter MCMC chains overall.
During the MCMC runs we iteratively update the β NL measurement (as it is cosmology dependent), as running the emulator at each step of the chain is computationally not feasible.Thus the β NL measurement is evaluated using the median of the current position of the walkers in the parameter space.This returns an effective value for the non-linear halo bias correction that is, over the run of the MCMC, representative of the median of corrections that would be applied to every single model iteration in the chain.In our pipeline, the number of steps between iterations can be set by the user and we find that updating the β NL values every 20 steps allows for a reasonable run time while providing enough updates to the β NL correction.On the other hand, the covariance matrix is only re-evaluated with the new parameters at the end of the MCMC run and checked.We find that the updated covariance matrix and halo model parameters do not affect the results of our fit as our initial cosmological and HOD parameters were set to the one of Heymans et al. (2021) and van Uitert et al. (2016) which are close to our final results.The covariance matrix is dominated by the shot and/or shape noise on the majority of scales.
To report our result we use two methods to estimate our constraints and parameter values.One method uses the maximum statistics of the marginal posterior distributions for each parameter (MMAX).Here the asymmetric errors are estimated around the maximum point in iso-distribution levels to cover 68% of the marginal distribution.For the second method we use the full posterior distribution to find the best-fitting parameters (maximum a posteriori point, MAP) and use the methodology presented in Joachimi et al. (2021) to associate an error to this measurement with the projected highest posterior density (PJ-HPD) approach.While the former method produces more stable parameter errors, especially when the likelihood surface is sparsely sampled, its point estimates, in general, do not correspond to the best fitting parameter values.In contrast the latter method will in general produce noisier error estimates with unbiased parameter values.

Results
Now we turn the focus to our results, presenting cosmological parameter constraints in Sect.4.1, large-scale analysis in Sect.4.2, constraints on the galaxy-halo connection in Sect.4.3, and effect of modelling of satellite galaxies in Sect 4.4.Further details are presented in Appendix E. To recap, our theoretical 2×2pt+SMF model consists of 17 free parameters, two of which are our main cosmology parameters, with 3 more secondary cosmology parameters that are harder to constrain given the combination of observables, and 11 parameters describing the galaxyhalo connection in the form of the CSMF.With six stellar mass bins, and three observables, our combined data vector consists of 156 data points.In Appendix C we use mock data realisations to estimate the effective number of degrees of freedom for our analysis finding ν eff = 147.55.Following the likelihood analysis described in Sect.3.6 we are able to constrain 12 parameters, listed in Table 2 along with their prior ranges.We find the MAP to provide a good fit4 to the data, with a reduced χ 2 value of 1.07 and p(χ 2 |ν eff ) = 0.27.
We compare the prediction from our fiducial model, and its 68% confidence regions, to the measured galaxy abundance SMF in Fig. 2, the galaxy-galaxy lensing ESD in Fig. 3, and the galaxy clustering w p in Fig. 4, for all of the six stellar mass bins.We find that the model reproduces the overall trends in the data, such as the presence of the bump at ∼ 1h −1 Mpc in ESD due to satellite galaxies, and the fact that the stronger signal is present where galaxies have higher stellar mass, showing that massive galaxies reside in more massive haloes.We note that some caution is needed when interpreting the results, as the quality of the fit cannot be judged by eye due to highly correlated data points.We find acceptable fits to each component of our 2×2pt+SMF data vector (for details see Appendix C).We note that the poorest fit is found for the w p section of our joint data vector with p[χ 2 (w p )|ν w p eff ] = 0.01.Whilst a formally acceptable fit, this may indicate that our model is lacking the ability to correctly describe the photometric redshift dilution effect discussed in Sect.3.3.

Cosmology constraints
We find the following cosmological parameter constraints from our simultaneous 2×2pt+SMF analysis of the ESD, w p and SMF  Notes: This table lists all the free parameters in our model: the energy density of cold matter Ω m , the normalisation of power spectrum σ 8 , the dimensionless Hubble parameter h, the spectral index n s , the energy density of baryonic matter Ω b , the derived parameter S 8 , the normalisation of the concentration-mass relation for dark matter haloes f h , the normalisation of stellar-tohalo mass relation M 0 , the characteristic scale of the stellar-tohalo mass relation M 1 , the slope parameters of the stellar-to-halo mass relation γ 1 and γ 2 , the scatter between stellar mass and halo mass σ c , the normalisation of the concentration-mass relation for distribution of satellite galaxies f s , the power law behaviour of satellites α s , the normalisation constants of the Schechter function b 0 and b 1 , and the Poisson parameter P. Parameters are deemed unconstrained when the marginal probability at 2σ level exceeds 13% of the peak probability (see appendix A of Asgari et al. 2021b).In cases where one side is constrained we report the 1σ lower/upper limit.The MMAX estimate is the marginal maximum statistic, reporting the point of maximum marginal posterior distribution to the iso-posterior levels above and below the maximal point.The MAP+PJ-HPD (maximum posterior with projected joint highest posterior density) estimates are calculated following Joachimi et al. (2021).where S 8 = σ 8 √ Ω m /0.3, and we quote the maximum statistics of the marginal posterior distributions (MMAX).The remaining cosmological parameters are unconstrained by our analysis, and informed by our choice of prior (see Fig . E.2).In Fig. 5 we present the 68% and 95% confidence levels of the joint twodimensional, marginalised posterior distribution in the S 8 − Ω m plane.The 2×2pt+SMF constraints are shown to be in good agreement with constraints from KiDS cosmic shear and KiDS with BOSS 3×2pt constraints (Asgari et al. 2021b;Heymans et al. 2021).They are formally consistent, but in some mild tension with the Planck Collaboration et al. (2020) TT,TE,EE+lowE CMB results.Using the Hellinger distance as a tension measure (see Heymans et al. 2021), the mild tension between our fiducial results and Planck is 1.9σ in S 8 .
In Fig. 6 we compare our constraints to a broader range of joint-probe large-scale structure analyses, finding consistency with all.Cacciato et al. (2013) and More et al. (2015) adopt a similar methodology to this study, using a halo model formalism to jointly analyse galaxy-galaxy lensing, galaxy clustering and galaxy abundance observables.Cacciato et al. (2013) analysed the Sloan Digital Sky Survey (SDSS-DR7).More et al. (2015) combined data from the Baryon Oscillation Spectroscopic Survey (BOSS) and the Canada France Hawaii Telescope Lensing Survey (CFHTLenS).The improved constraining power in this KiDS analysis reflects the significant increase in depth for the lensing sample relative to SDSS-DR7, and the ten-fold increase in area relative to CFHTLenS.We can also compare to 2×2pt analyses that combine galaxy-galaxy lensing and galaxy clustering that introduce conservative scale cuts to reflect known limitations their adopted galaxy bias models.These include the fiducial analyses from DES, which adopt a linear galaxy bias model5 (Porredon et al. 2022;Pandey et al. 2022), and from the Hyper Suprime Camera survey (HSC, Miyatake et al. 2021).The HSC analysis is highly complementary to our study as both analyses used a form of halo model with a halo occupation distribution.Miyatake et al. (2021) use the Zheng et al. (2005) HOD built into the dark-matter N-body DarkQuest emulator to predict the 2×2pt observables (Nishimichi et al. 2019;Miyatake et al. 2022).We use the same emulator to calibrate the non-linear halo bias in our model (see Sect. 2.1).Compared to this analysis, HSC adopts more conservative scale cuts arising from concern over unmodeled baryon feedback on small-scales, which our methodology can account for through the free normalisation  2021) also choose scale cuts to mitigate small-scale assembly bias for the relatively rare luminous red galaxies in their sample, which cannot be modelled using an HOD.Owing to the volumelimited mix of all galaxies used in our KiDS-Bright analysis, we consider any assembly bias to be a subdominant effect in our theoretical model.The constraining power of the KiDS 2×2pt+SMF analysis in the S 8 − Ω m plane is the same as that of the 3×2pt studies from DES Collaboration et al. ( 2022), with σ S 8 σ Ω m [DES 3×2pt] /σ S 8 σ Ω m [KiDS 2×2pt+SMF] = 0.97.This may be surprising given the five-fold increase in area for DES relative to KiDS, and the addition of the cosmic shear probe in the 3×2pt analysis.This comparison therefore highlights the significant constraining power from the inclusion of non-linear scales in the 2×2pt+SMF analysis that are excluded from the DES 3×2pt analysis.Comparing KiDS 2×2pt+SMF constraints to the KiDS with BOSS 3×2pt analysis (Heymans et al. 2021), we first review the BOSS spectroscopic clustering constraints where 03. Finding the same constraining power between these analyses may again be surprising, given the nine-fold increase in area for BOSS relative to KiDS.As such it demonstrates the significant constraining power from non-linear scales when the galaxy bias can be constrained using galaxy-galaxy lensing and galaxy abundance.Comparing to the full 3×2pt analysis we find 39, where the extra constraining power in the 3×2pt analysis is driven by the cosmic shear.Future studies with KiDS will combine 2×2pt+SMF with cosmic shear data, including further development and validation of our adopted halo model methodology.
In Appendix E we explore a number of extensions to our fiducial analysis.In Appendix E.1 we quantify the expected contamination to our observables from intrinsic galaxy alignments and magnification.The contamination levels are found to be negligible relative to our statistical errors, justifying our choice to not account for these astrophysical effects in our model.In Appendix E.2 we demonstrate that our S 8 and Ω m constraints are insensitive to our choice of prior on n s .In Appendix E.3 we quantify the bias in Ω m without the inclusion of our nuisance parameter D to model photometric redshift dilution in our galaxy clustering measurement.In Appendix 4.4 we quantify the impact of assumptions governing the behaviour of satellite galaxies.

Large-Scale Analysis
Amon et al. (2022b) present a detailed analysis of the uncertainty on the amplitude of the small-scale galaxy-galaxy lensing and galaxy clustering signal for BOSS galaxies that arises from our imperfect knowledge of baryon feedback and assembly bias.They conclude that the introduction of scale cuts with r p > 5h −1 Mpc fully isolates these effects.Following the methodology in Appendix C we determine ν r p >5 eff = 60.31 for a largescale only 2×2pt+SMF data vector analysis, calculating the large-scale goodness of fit of our fiducial best-fit model (see Table 2) to be p[χ 2 (r p > 5)|ν r p >5 eff ] = 0.19.As such we find no sign of tension between our fiducial all-scale analysis and a restricted large-scale analysis.
The majority of the information in our (r p > 5) data vector comes from the SMF as there are only 4 highly correlated data points remaining in the ESD and w p measurements, per stellar mass bin.We found that a full likelihood analysis of the (r p > 5) data vector was unable to converge in our highly-flexible 17parameter model space.Where larger scales are used exclusively, either more precise data or the use of a less flexible model is necessary (such as in More 2013; Miyatake et al. 2021;Amon et al. 2022b, amongst others).The extra flexibility afforded in  Leauthaud et al. 2012Moster et al. 2013Coupon et al. 2015van Uitert et al. 2016 Fig. 7: The stellar-to-halo mass relation, as defined by Eq. 19 using the best-fit HOD parameters from our 2×2pt+SMF analysis (blue).The result can be compared to results from Leauthaud et al. (2012, green), Moster et al. (2013, black), Coupon et al. (2015, orange), andvan Uitert et al. (2016, purple).The grey area shows the range in stellar masses where the obtained stellar-to-halo mass relation is extrapolated beyond the range of median stellar masses used in this analysis.our model is, however, essential when analysing small scales in order to capture baryonic effects (Debackere et al. 2020(Debackere et al. , 2021) )

Galaxy-halo connection
The powerful aspect of the method used in this work is that it is able to simultaneously constrain both cosmological parameters as well as the halo occupation statistics.The full results are listed in Table 2 and the marginalised posterior distributions of all the parameters are shown in Fig. E.2.The HOD parameters are tightly constrained, with some strong degeneracies between the parameters M 0 , M 1 and γ 1 that govern the characteristic mass of the stellar-to-halo mass relation knee and the high mass slope of centrals, and between the b 0 and b 1 parameters, which govern the normalisation of the satellite CSMF.The first degeneracy is somewhat expected given the data, as the stellar mass function at the high mass end is highly uncertain and dominated by the Eddington bias.The second degeneracy also arises from the fact that both parameters compete for the overall normalisation of the satellite CSMF.
We find the parameters of our HOD model to be in good agreement with the previous studies using GAMA-like galaxies (van Uitert et al. 2016;Dvornik et al. 2018;Bilicki et al. 2021).In order to show the agreement in a more intuitive way, we take the parameters of the stellar-to-halo mass relation and combine them using the same functional form (Eq. 19), which results in the relation shown in Fig. 7.In the same figure we show good agreement with the results from van Uitert et al. (2016), where the HOD parameters were constrained using galaxy-galaxy lensing combined with a SMF for KiDS and GAMA data, adopting a fixed Planck cosmology.We also find qualitative agreement with constraints from abundance matching to numerical simulations (Moster et al. 2013) and constraints from a 2×2pt+SMF analysis of COSMOS with a fixed WMAP5 cosmology (Leauthaud et al. 2012, note here we compare to constraints from the most similar 0.22 < z < 0.48 COSMOS sample) as well as the CFHTLenS/VIPERS analysis by Coupon et al. (2015).We find our constraints on the scatter of the central CSMF, σ c , and the low mass end slope of the satellites, α s , to be in agreement with Yang et al. (2009);Cacciato et al. (2013);van Uitert et al. (2016); Dvornik et al. (2018); Bilicki et al. (2021).The Eddington bias at the high mass end is captured by the σ c parameter, leaving the other parameters mostly unaffected.
We account for the impact of baryon feedback in our model by allowing for freedom in the normalisation of the massconcentration relation for both the haloes, f h , and the satellite galaxy distribution, f s (Eq.15).With these independent free parameters we can capture the expected small-scale baryon feedback power suppression shown in hydrodynamical simulations (Debackere et al. 2021;Amon et al. 2022b).We find f h and f s to be consistent with 1, with a preference for lower values, with 1σ lower limits f h > 0.65 and f s > 0.38.Our results are consistent with Viola et al. (2015), indicating that the concentrations of real haloes and satellite distributions are smaller than the haloes in dark matter only simulations (see also Debackere et al. 2020Debackere et al. , 2021)).Future work will compare direct measurements from hydrodynamical simulations (e.g.McCarthy et al. 2017) with our halo model approach to account for the mass dependence of baryonic effects on the radial profiles of dark matter haloes.

Modelling satellite galaxies
In our fiducial model we used the findings of Dvornik et al. (2018) that showed that the occupation distribution of satellite galaxies does not follow a Poisson distribution, and that generally the parameter P (Eq.25) is not unity, with our fiducial run preferring a sub-Poissonian behaviour.Following Cacciato et al. (2013), we quantify the impact of removing this flexibility in the model, by fixing the parameter P to unity, thus assuming that the satellite galaxies obey the Poisson distribution.We run another set of MCMC chains using the same setup as in the fiducial case, but with one less parameter to constrain.The resulting constraints are shown in Fig. E.2, with marginalised constraints quoted in Table E.1.We find significant shifts for the two main cosmological parameters with Ω m = 0.330±0.019and S 8 = 0.951 +0.037 −0.036 and a formally acceptable fit with p(χ 2 |ν eff ) = 0.02 for the whole data vector.In this case we find that the fixed Poisson parameter non-trivially affects the other parameters governing the satellites in the halo model.Specifically the normalisation of the satellite conditional stellar mass function, b 0 and b 1 , shifts, and these parameters are in turn non-trivially correlated with the main cosmological parameters.Cacciato et al. (2013) argues that flexibility in the form of the satellite galaxy model is critically important in order to both constrain the galaxy bias (Cacciato et al. 2012;Dvornik et al. 2018;Asgari et al. 2021a), and to obtain unbiased cosmological parameters.
There are several reasons to reject the results of our Poissonian satellite distribution model analysis.First, whilst we find an acceptable fit of the P = 1 model to our full 2×2pt+SMF data vector, there is an unacceptable fit to the w p part of the data vector, with p[χ 2 (w p )|ν w p eff ] ∼ 10 −4 .Secondly, there is observational evidence from the GAMA group catalogue (Robotham et al. 2011) that for haloes with masses below 10 13 h −1 M ⊙ , the number of satellites exhibit sub-Poisson behaviour, where in Fig. 8 we measure P(M) as a function of the dynamical mass M dyn (Driver et al. 2022).We relate the parameter P with the observed mean and variance of the number of GAMA group members in narrow bins of dynamical mass M dyn as: Here the mean ⟨N s |M dyn ⟩ and variance Var[N s |M dyn ] are directly obtained from the GAMA groups catalogue, where we select groups with the number of members that is equal to or larger than 3 (Robotham et al. 2011).We find the satellite distribution to be sub-Poissonian for M dyn < 10 13 h −1 M ⊙ , ranging from P(M dyn ) = 0.15 at M dyn = 10 11 h −1 M ⊙ , to P(M dyn ) = 0.76 at M dyn = 10 13 h −1 M ⊙ , (shown blue in Fig. 8).Assuming the GAMA dynamical mass is a reasonable estimate of the halo mass M h , using Fig. 7 to define our stellar mass range, we expect the satellite distribution to be sub-Poissonian across the full stellar mass range of our 2×2pt+SMF analysis.This sub-Poissonian behaviour is recovered in our fiducial analysis where P = 0.403±0.029.In Dvornik et al. (2018) analysis the behaviour of satellite galaxies is recovered to be super-Poissonian, which is consistent with the trend seen in Fig. 8, as the halo masses of GAMA galaxies in that sample were above 10 13 h −1 M ⊙ .Finally, hydrodynamical simulations show that for haloes with masses above ∼ 5 × 10 12 h −1 M ⊙ , the number of satellites exhibit super-Poisson behaviour (Hadzhiyska et al. 2022) and for haloes with masses below ∼ 5 × 10 12 h −1 M ⊙ , the number of satellites exhibit sub-Poisson behaviour (Kravtsov et al. 2004;Jiang & van den Bosch 2017;Bhowmick et al. 2018;Beltz-Mohrmann et al. 2020).The studies from Beltz-Mohrmann et al. (2020) and Hadzhiyska et al. (2022) show a large uncertainty on the Poisson number and can be to an extent also well described by a model that assumes a Poisson distribution, while the results from Kravtsov et al. (2004), Jiang & van den Bosch (2017) and Bhowmick et al. (2018) show a clear sub-Poisson behaviour for low mass haloes, with the analysis of Jiang & van den Bosch (2017) also showing super-Poisson behaviour for high mass haloes with a transition period where satellite galaxies behave completely Poissonian.What is more, a recent analysis (Linke et al. 2022) of a semi-analytic simulation (Henriques et al. 2015) using galaxy-matter bispectrum shows that parameter P indeed varies from sub-Poisson behaviour to super-Poisson behaviour as halo mass increases, and it seems to be furthermore dependant on the stellar mass as well (Appendix B therein).
Given the significant impact of the form of the satellite galaxy model on our cosmology constraints, we explore the satellite distribution further, noting that Fig. 8 reveals a strong mass dependence that is missing from our fiducial model.We determine the impact of neglecting this mass-dependence in our analysis, by including an explicit halo mass dependence on the Poisson parameter P. We fit a power law function to the GAMA data as: finding A = 0.43, B = 0.39 and M piv = 12.54 when M = M dyn (shown in orange in Fig. 8).In our extended 2×2pt+SMF analy- Fig. 8: The Poissonian parameter P (Eq.52), as a function of dynamical group mass, M dyn , for GAMA galaxy groups with 3 or more members (blue), and groups with 5 or more members (white).The orange/green line shows the best fit power law (Eq.53) to the GAMA/2×2pt+SMF data, with the green region showing the uncertainty on the 2×2pt+SMF fit.The fiducial, single value, model for P is shown with the grey horizontal band.With the wiggly greyed our areas we show the approximate range of dynamical group masses that are outside of halo mass range of our analysis.sis, we model P(M) using Eq.53, fixing the normalisation A and slope B of the power law to the GAMA best-fit values.We treat M piv as a free parameter, however, to account for the uncertainty in the relationship between the GAMA dynamical mass, M dyn , and the true lensing halo mass used in the halo model.The best fit M piv from the 2×2pt+SMF is nevertheless found to be in good agreement with the GAMA fit (shown in green in Fig. 8).
Using the P(M) model we obtain the parameter constraints listed in Table 2 and shown in Fig E .2.We find a good fit of the model to the data with p(χ 2 |ν eff ) = 0.19.All the parameters are consistent with the constraints from our fiducial analysis that assumes no mass dependence.We find a preference for lower values for the two primary cosmological parameters, with S 8 = 0.718 +0.045 −0.040 and Ω m = 0.276 +0.024 −0.021 corresponding to a 1.3σ and 0.7σ difference from the fiducial result (see also Fig. E.1).
This extension analysis shows that our results are sensitive to how we model the mass-dependence of the satellite group distribution, at an acceptable level of ∼ 1σ in our primary cosmological parameters.We chose not to use this extended P(M) model in our fiducial analysis, as in Fig. 8, we show how sensitive the GAMA-measured P(M dyn ) relationship is to the group selection criteria.We find that the behaviour changes when groups are defined with a number of members that is equal or larger than 5 (white data points, compared to the blue data points for our original selection criteria).In future higher signal-to-noise studies we will explore keeping the parameters A and B free in the P(M) model, and implement a more complex model that uses the non-Poisson behaviour directly in the definition of the HOD (for instance a negative binomial distribution as shown in Boylan-Kolchin et al. 2010).Further considerations need to be taken into account for possible stellar mass dependency of the Poisson parameter P (Linke et al. 2022).

Discussion and conclusions
In this paper we combined measurements of galaxy clustering, galaxy-galaxy lensing and galaxy abundances in the form of the stellar mass function, in order to simultaneously set constraints on cosmological parameters and galaxy bias.Using a flexible halo model, we analysed the fourth data release of the Kilo-Degree Survey (KiDS) (KiDS-1000, Kuijken et al. 2019), where the source sample, used to measure the galaxy-galaxy lensing signal, has undergone a rigorous study to assess robustness and accuracy (Asgari et al. 2021b;Giblin et al. 2021;Hildebrandt et al. 2021).For our lens sample we used the KiDS-Bright sample (Bilicki et al. 2021) whose selection was calibrated against a complete and representative spectroscopic sample from the GAMA survey (Driver et al. 2011), with the photometric redshifts calibrated using a neural network (ANNz2, Sadeh et al. 2016).The resulting accuracy of the estimated redshifts for the KiDS-bright sample is sufficient for galaxy-galaxy lensing and galaxy clustering studies.
We used the halo model to analyse our data, building upon the cosmological analyses presented in Cacciato et al. (2013) and More et al. (2015).We used a single halo occupation model to compute the clustering of galaxies and the galaxy-galaxy lensing signal (Guzik & Seljak 2002;Yoo et al. 2006;Cacciato et al. 2009), and the galaxy abundances (van den Bosch et al. 2013;Cacciato et al. 2013).This model was shown to be able to simultaneously constrain cosmology and halo occupation statistics, as well as constrain the extensions to the standard ΛCDM cosmologies, such as the equation of state of dark energy and neutrino mass (More et al. 2013;Krause & Eifler 2017).
We improved upon previous studies by (i) using a more accurate N-body simulation calibrated analytical model, taking into account halo exclusion, scale dependence and the non-linear nature of halo bias (Mead & Verde 2021;Mahony et al. 2022), (ii) using the latest lensing and clustering data from a single survey (KiDS-1000), and (iii) using the full analytical covariance matrix that accounts for cross-covariance between all observables and in particular the cross-covariance between the stellar mass function and two point statistics.
We have adopted a Bayesian approach to constrain our model parameters, using the MCMC to probe the posterior distributions.For a flat ΛCDM cosmology we find Ω m = 0.290 +0.021 −0.017 and S 8 = 0.773 +0.028 −0.030 , which is consistent and comparable to constraints from 3×2pt studies that also include a cosmic shear observable (Heymans et al. 2021;DES Collaboration et al. 2022).Our results follow the trend seen in other lensing studies for a tension in S 8 when compared to Planck Collaboration et al. (2020).Using the Hellinger distance as a tension metric, this difference is at the 1.9σ level for our 2×2pt+SMF analysis.We find that our constraints are sensitive, at the ∼ 1σ-level (in S 8 ), to how we choose to model the mass-dependence of the satellite distribution within the halo model.This aspect of our analysis will require further development in future higher signal-to-noise studies.
Combining galaxy clustering and galaxy-galaxy lensing with cosmic shear measurements has been a standard approach for large-scale structure analyses in recent years (van Uitert et al. 2018;Joudaki et al. 2017;Abbott et al. 2018;DES Collabora-tion et al. 2022;Heymans et al. 2021).We anticipate that combing our halo model approach with cosmic shear data will allow for additional constraints on astrophysical systematics arising from the intrinsic alignment of galaxies and baryon feedback.So far, intrinsic alignments in cosmic shear analysis have either been included using a non-linear modification of the linear alignment model (Bridle & King 2007, NLA) or a perturbation theory approach (TATT Blazek et al. 2019).For a consistent halo model approach this effect could also be modelled within the framework adopted in this analysis (e.g.Fortuna et al. 2021).
In this analysis we have varied the halo concentration parameter to account for baryon feedback.With additional data, a more complex halo model could be adopted allowing for the inclusion of gas observables to constrain baryon feedback through the Sunyaev-Zeldovich effect (Mead et al. 2020;Tröster et al. 2022).The flexibility of the halo model also allows for extensions to the underlying cosmological model without having to employ a myriad of costly simulations to cover a large range of parameters (Cataneo et al. 2019).We therefore see a significant role for our adopted methodology in future cosmological analyses of upcoming large-scale structure surveys..1, but for the 3 different observables in our analysis.The fits to the distributions are used to determine the number of degrees of freedom for each observable which are in turn are used to determine the reduced χ 2 values for each of them respectively.Note that the χ 2 distributions are not a result of a maximising the posterior for that section of the mocks.The vertical lines show the χ 2 values from the real data, matching in colours with the histograms.
sufficiently small that we do not include it in our overall error budget.We find that the random correction applied to the data exceeds 100% in the first three stellar mass bins on large scales, and can be up to 4 times larger than the statistical uncertainty.2).The contours indicate 1σ and 2σ confidence regions.
We present the contribution for the largest stellar mass bin only, as the contribution to magnification and intrinsic alignments is of the same order for all of them.On the same figure we also show the changes to the galaxy-galaxy lensing and galaxy clustering signals if we change the two main cosmological parameters.The contributions of both intrinsic alignments and magnification are well below 1%, as also found in Unruh et al. (2020) and are well within the error budget of the data.Moreover, they are also subdominant to the changes due to the cosmological parameters of interest.Such contributions have negligible effects on the overall signal, and are unlikely to be a significant source of bias.To some extent the effects cancel each other out, as for the redshifts of our lens galaxies the magnification effect dilutes the signal, while the intrinsic alignments add a similar contribution. of parameter space.Given that there are no strong degeneracies between this set and the rest of the parameters, and that the set is already well constrained from other studies (see the discussion in appendix B of Heymans et al. 2021) we find no motivation to investigate the impact of widening the priors.Instead we investigated reducing the prior range by adopting a Gaussian prior with the mean and uncertainty fixed to the Planck Collaboration et al. ( 2020) constraint for n s .This parameter is the most interesting to investigate of the three, as any tension between small and large scales may be expected to manifest in a biased spectral index constraint (Tröster et al. 2021).and Ω m , showing 68% and 95% credible regions.We compare our fiducial analysis (blue) with an analysis neglecting the impact of photometric redshift uncertainties which dilute the estimated projected galaxy clustering signal (grey).
The results are presented in Table E.1 and Fig. E.2.We find that the marginal contours for the parameters do not change significantly with the addition of a restrictive prior on n s .For example, the MMAX estimate of S 8 , shifts by 0.13σ, which is consistent with the expected MCMC run-to-run variance (Joachimi et al. 2021).In future studies we will explore the impact of using more restrictive priors on all externally-constrained parameters which we are insensitive to, noting that this may also help to reduce projection bias for posteriors with many parameters.
Appendix E.3: The effect of unmodelled photometric redshift errors on the projected galaxy clustering measurements In Sect.3.3 we introduce a free nuisance parameter, D (Eq.43), to account for our uncertainty on the amplitude of the true projected galaxy clustering signal.The expected dilution is a result of unaccounted photometric redshift errors in our theoretical model for w p (r p ).In Figure E.4 we compare our fiducial constraints in S 8 and Ω m with the constraints from an analysis where the photometric redshift dilution effect is neglected and D = 1.We find that while the omittance of the photometric redshift dilution effect does not impact the S 8 constraints, Ω m becomes biased and more constrained.This is expected as the galaxy clustering is more sensitive to Ω m compared to σ 8 (cf.
spectroscopic redshifts with multi-wavelength imaging from 21 broadband filters from the far-UV to the far-IR.The median offset is M KiDS ⋆ /M GAMA ⋆ = −0.09± 0.18 dex.Brouwer et al. ( is based on from an

Fig. 1 :
Fig.1: Galaxy stellar mass as a function of ANNz2 photometric redshift for the KiDS-Bright sample.The full sample is shown with a logarithmic hexagonal density plot.The blue line shows the stellar mass limit determined using the automated method presented byWright et al. (2017).Red boxes show the six stellar mass bins used in the analysis, with individual galaxies plotted as black dots.The bin ranges were chosen in such a way as to achieve a good signal-to-noise ratio in all bins for our galaxygalaxy lensing and galaxy clustering measurements.

Fig. 2 :
Fig. 2: The KiDS-Bright galaxy stellar mass function and the fractional errors.Upper panel: The KiDS-Bright galaxy stellar mass function (crosses) compared to the model from Wright et al. (2018), evaluated at the median redshift of our sample (black dashed).The blue line and shaded region indicates the best fit model and 68% confidence levels of our best fit halo model (Eq.22).We caution that the quality of the fit cannot be judged by eye, because of the covariance in the data, between the data points and between the other observables.The reduced χ 2 value for this observable is 1.05 (DoF = 14.58, p-value = 0.39), estimated using the method presented in Appendix C. Lower panel: The fractional errors on the data and the model, ∆Φ/δΦ.

3. 3 .
Projected galaxy clustering measurements: w p We measure the clustering of the KiDS-Bright galaxy sample using the Landy-Szalay (Landy & Szalay 1993) estimator for the galaxy correlation function ξgg (r p , r π ) = DD − 2DR + RR RR r p ,r π .
Fig. C.2: Same as Fig. C.1, but for the 3 different observables in our analysis.The fits to the distributions are used to determine the number of degrees of freedom for each observable which are in turn are used to determine the reduced χ 2 values for each of them respectively.Note that the χ 2 distributions are not a result of a maximising the posterior for that section of the mocks.The vertical lines show the χ 2 values from the real data, matching in colours with the histograms.
Fig. E.2: The full posterior distributions of the model parameters (where the priors are listed in Table2).The contours indicate 1σ and 2σ confidence regions.

Fig. E. 3 :
Fig. E.3:  Relative difference between the lensing signal with and without the contribution of intrinsic alignments and lens magnification.The intrinsic alignment contribution is modelled with A IA = 1, which is a reasonable 'worst case' scenario for the KiDS-Bright sample, and the magnification is modelled with α d = 0.85 and α d = 2.11, typical for GAMA-like lenses at a redshift of 0.21 and 0.36.Grey areas show the error on the data.We also show the changes to the galaxy-galaxy lensing and galaxy clustering signals if we change the two main cosmological parameters.The behaviour is only shown for the largest stellar mass bin, as the contribution to magnification and intrinsic alignments is of the same order for all of them.
Fig. E.3).This motivates future work to improve the estimator for, or theoretical modelling of, the projected galaxy clustering signal in the presence of photometric redshift errors, followingJoachimi et al. (2011);Chisari et al. (2014).

Table 2 :
Marginal constraints on all model parameters listed together with their priors.
(Porredon et al. 2022;Pandey et al. 2022)d Ω m constraints.Our fiducial results (KiDS 2×2pt +SMF) can be compared to a similar 2×2pt +SMF analysis from SDSS(Cacciato et al. 2013); a series of 2×2pt studies including the latest results from DES(Porredon et al. 2022;Pandey et al. 2022)and HSC (Miyatake et al. 2021); 3×2pt analyses from KiDS with BOSS (Heymans et al. 2021) and DES (DES Collaboration et al. 2022).The last entry shows the Planck Collaboration et al. (2020, TT,TE,EE+lowE) constraints.Our results are consistent with all studies, including Planck Collaboration et al. (2020), although we find a mild tension between our S 8 constraints and those from Planck, at the level of 1.9σ.