© ESO 2021

1. Introduction

Observations of the cosmic microwave background (CMB) have delivered high-precision constraints for the cosmological parameters of the flat, cold dark matter, and cosmological constant model of the Universe (ΛCDM, Planck Collaboration VI 2020). With only six free parameters, this flat ΛCDM model provides an exquisite fit to observations of the anisotropies in the CMB. The same model predicts a range of different observables in the present-day Universe, including the cosmic expansion rate (Weinberg 1972), and the distribution of, and gravitational lensing by, large-scale structures (Peebles 1980; Bartelmann & Schneider 2001; Eisenstein et al. 2005). In most cases there is agreement between the measured cosmological parameters of the flat ΛCDM model, when comparing those constrained at the CMB epoch with those constrained through a variety of lower-redshift probes (see the discussion in Planck Collaboration VI 2020, and references therein). Recent improvements in the statistical precision of the lower-redshift probes have, however, revealed some statistically significant differences. Most notably, a 4.4σ difference in the value of the Hubble constant, H₀, has been reported using distance ladder estimates in Riess et al. (2019). If this difference cannot be attributed to systematic errors in either experiment, or both, this result suggests that the flat ΛCDM model is incomplete.

Many extensions have been proposed to reconcile the observed differences between high- and low-redshift probes (see for example Riess et al. 2016; Poulin et al. 2018; Di Valentino et al. 2020). All, however, require additional components to the cosmological model that move it even further away from the standard model of particle physics, a model that already struggles to motivate the existence of cold dark matter and a cosmological constant. As the statistical power of the observations continues to improve, focus has shifted to establishing a full understanding of all systematic errors and to the development of mitigation approaches, in preparation for the high-precision ‘full-sky’ imaging and spectroscopic cosmology surveys of the 2020s (Euclid, Laureijs et al. 2011; LSST Science Collaboration 2009; DESI Collaboration 2016).

We present a multi-probe ‘same-sky’ analysis of the evolution of large-scale structures, using overlapping spectroscopic and optical-to-near-infrared imaging surveys. Our first observable is the weak gravitational lensing of background galaxies by foreground large-scale structures, known as ‘cosmic shear’. Our second observable is the anisotropic clustering of galaxies within these large-scale structures, combining measurements of both redshift-space distortions and baryon acoustic oscillations. Our third observable is the weak gravitational lensing of background galaxies by the matter surrounding foreground galaxies, known as ‘galaxy-galaxy lensing’. As these three sets of two-point statistics are analysed simultaneously, this combination of probes is usually referred to as a ‘3 × 2 pt’ analysis.

Each observable in our multi-probe analysis is subject to systematic uncertainties. For a cosmic shear analysis, the observable is a combination of the true cosmological signal with a low-level signal arising from the intrinsic alignment of galaxies, as well as potential residual correlations in the data induced by the atmosphere, telescope, and camera. The signal can also be scaled by both shear and photometric redshift measurement calibration errors (see Mandelbaum 2018, and references therein). For a galaxy clustering analysis, the observable is the true cosmological signal modulated by an uncertain galaxy bias function. This function maps how the galaxies trace the underlying total matter distribution (see Desjacques et al. 2018, and references therein). It can be non-linear and evolves with redshift. The cosmological clustering also needs to be accurately distinguished from artificial clustering in the galaxy sample, arising from potentially uncharacterised inhomogeneities in the target selection (see for example Ross et al. 2012). Finally, the galaxy-galaxy lensing analysis is subject to the systematics that impact both the cosmic shear and clustering analyses.

When analysing these observables in combination the different astrophysical and systematic dependencies allow for some degree of self-calibration (Bernstein & Jain 2004; Hu & Jain 2004; Bernstein 2009; Joachimi & Bridle 2010). Adopting ‘same-sky’ surveys, in which imaging for weak lensing observables overlaps with spectroscopy for anisotropic galaxy clustering measurements, also allows for their cross-correlation. Such a survey design therefore presents a robust cosmological tool that can calibrate and mitigate systematic and astrophysical uncertainties through a series of nuisance parameters. In addition to enhanced control over systematics, this combination of probes breaks cosmological parameter degeneracies from each individual probe. For a flat ΛCDM model, this leads to significantly tighter constraints on the matter fluctuation amplitude parameter, σ₈, and the matter density parameter, Ω_m, whilst also decreasing the uncertainty on the recovered dark energy equation of state parameter in extended cosmology scenarios (Hu & Jain 2004; Gaztañaga et al. 2012).

Three variants of a joint ‘3 × 2 pt’ analysis have been conducted to date. van Uitert et al. (2018) present a joint power-spectrum analysis of the Kilo-Degree Survey (KiDS, Kuijken et al. 2015) with the Galaxy And Mass Assembly survey (GAMA, Liske et al. 2015), incorporating projected angular clustering measurements. Joudaki et al. (2018) present a joint analysis of KiDS with the 2-degree Field Lensing Survey (2dFLenS, Blake et al. 2016) and the overlapping area in the Baryon Oscillation Spectroscopic Survey (BOSS, Alam et al. 2015), incorporating redshift-space clustering measurements. Abbott et al. (2018) present a joint real-space lensing-clustering analysis of the Dark Energy Survey (DES Y1, Drlica-Wagner et al. 2018), using a high-quality photometric redshift sample of luminous red galaxies for their projected angular clustering measurements. In all three cases a linear galaxy bias model was adopted.

In this analysis we enhance and build upon the advances of previous ‘3 × 2 pt’ studies. We analyse the most recent KiDS data release (KiDS-1000, Kuijken et al. 2019), more than doubling the survey area from previous KiDS studies. We utilise the full BOSS area and the ‘full-shape’ anisotropic clustering measurements of Sánchez et al. (2017), incorporating information from both redshift-space distortions and the baryon acoustic oscillation as our galaxy clustering probe. We adopt a non-linear evolving galaxy bias model, derived from renormalised perturbation theory (Crocce & Scoccimarro 2006; Chan et al. 2012). We maximise the signal-to-noise in our KiDS-BOSS galaxy-galaxy lensing analysis, by including additional overlapping spectroscopy of BOSS-like galaxies from 2dFLenS.

This paper is part of the KiDS-1000 series. The KiDS-1000 photometry and imaging is presented in Kuijken et al. (2019). The core weak lensing data products are presented and validated in Giblin et al. (2021, shear measurements), and Hildebrandt et al. (2021, redshift measurements). Asgari et al. (2021a) conduct the cosmic shear analysis using a range of different two-point statistics, and Joachimi et al. (2021) detail the methodology behind our ‘3 × 2 pt’ analysis, with a particular focus on pipeline validation and accurate covariance matrices. In this analysis we constrain the cosmological parameters of the flat ΛCDM model. A range of different extensions to the ΛCDM model are considered in Tröster et al. (2020a), including varying dark energy, neutrino mass, spatial curvature and various modified gravity scenarios (Bose et al. 2020).

In this paper we review the data and provide a concise summary of the findings of the KiDS-1000 series of papers in Sect. 2. We present our joint cosmological constraints in Sect. 3 and conclude in Sect. 4. Appendices tabulate the galaxy properties (A), the adopted cosmological parameter priors (B), and the cosmological parameter constraints (C). They also discuss: the choice of intrinsic galaxy alignment model (D); a series of sensitivity tests (E); the expected differences between parameter constraints for overlapping weak lensing surveys (F); a range of different ‘tension’ metrics (G); the redundancy, validation and software review for our pipeline (H); and the minor analysis additions that were included after the analysis was formally unblinded (I).

2. Data and methodology

2.1. Surveys: KiDS, BOSS, and 2dFLenS

The Kilo-Degree Survey (KiDS, de Jong et al. 2013), covers 1350 deg² split into two fields, one equatorial and one southern. Matched-depth imaging in nine bands spans the optical, ugri, through to the near-infrared, ZYJHK_s, where the near-infrared imaging was taken as part of the KiDS partner survey VIKING (the VISTA Kilo-degree INfrared Galaxy survey, Edge et al. 2013). High-quality seeing was routinely allocated to the primary KiDS r-band VST-OmegaCAM observations, resulting in a mean r-band seeing of 0.7 arcsec, with a time-allocated maximum of 0.8 arcsec. This combination of full-area spatial and wavelength resolution over a thousand square degrees provides a unique weak lensing survey that allows for enhanced control of systematic errors (Giblin et al. 2021; Hildebrandt et al. 2021). This analysis uses data from the fourth KiDS data release of 1006 deg² of imaging, (hence the name KiDS-1000), which has an effective area, after masking, of 777 deg². KiDS is a public survey from the European Southern Observatory, with data products freely accessible through the ESO archive¹.

Giblin et al. (2021) present a series of null-tests to validate the KiDS-1000 shear catalogue² in five tomographic bins spanning a photometric redshift range of 0.1 < z_B ≤ 1.2 (see Appendix A for details of the properties of each bin). Meeting their requirement that any systematic detected induces less than a 0.1σ change in the inferred cosmic shear constraints on the clustering cosmological parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, they conclude that the shear catalogue is ‘science-ready’, with no significant non-lensing B-mode distortions detected. Kuijken et al. (2015) present the catalogue-level blinding methodology that we adopted to introduce ±2σ differences in the recovered value of S₈ in order to retain team ignorance over the final cosmological results until all analysis decisions were finalised (for further details on blinding see Appendix I). Hildebrandt et al. (2021) present the KiDS-1000 photometric redshift calibration. This is determined using the self-organising map (SOM) methodology of Wright et al. (2020a), and is validated with a cross-correlation clustering analysis, following van den Busch et al. (2020). The SOM identifies and excludes any galaxies that are poorly represented in the spectroscopic calibration sample, in terms of their nine-band colours and magnitudes. The resulting ‘gold’ photometric sample, with an accurately calibrated redshift distribution, is then re-simulated in the KiDS image simulations of Kannawadi et al. (2019) in order to determine the shear calibration corrections for each tomographic bin, and an associated uncertainty (see Giblin et al. 2021; Hildebrandt et al. 2021, for full details).

The Baryon Oscillation Spectroscopic Survey (BOSS, Alam et al. 2015), spans an effective area of 9329 deg², with spectroscopic redshifts for 1.2 million luminous red galaxies (LRG) in the redshift range 0.2 < z < 0.9. A range of different statistical analyses of the clustering of BOSS galaxies have been used in combination with CMB measurements, to set tight constraints on extensions to the standard flat ΛCDM model (see Alam et al. 2017; eBOSS Collaboration 2020, and references therein). We adopt the anisotropic clustering measurements of Sánchez et al. (2017) in this multi-probe analysis. BOSS only overlaps with the equatorial stripe of the KiDS survey, with 409 deg² of the BOSS survey lying within the KiDS-1000 footprint. BOSS galaxies in this overlapping region are used as lenses in our galaxy-galaxy lensing analysis, with an effective lens number density of 0.031 arcmin⁻² (see Appendix A for details). BOSS is a public survey from the third Sloan Digital Sky Survey (York et al. 2000), and we analyse data from the twelfth data release³ (DR12, Alam et al. 2015).

The 2-degree Field Lensing Survey (2dFLenS, Blake et al. 2016), spans 731 deg², with spectroscopic redshifts for 70 000 galaxies out to z < 0.9. This galaxy redshift survey from the Anglo-Australian Telescope (AAT) was designed to target areas already mapped by weak lensing surveys to facilitate ‘same-sky’ lensing-clustering analyses (Johnson et al. 2017; Amon et al. 2018; Joudaki et al. 2018; Blake et al. 2020). We use data from the 2dFLenS LRG sample that was targeted to match the BOSS-LRG selection, but with sparser sampling. 2dFLenS thus provides an additional sample of BOSS-like galaxies in the KiDS southern stripe where there is 425 deg² of overlap within the KiDS-1000 footprint. 2dFLenS galaxies in this overlapping region are used as lenses in our galaxy-galaxy lensing analysis, with an effective lens number density of 0.012 arcmin⁻² (see Appendix A for details). 2dFLenS was an AAT Large Programme that has been made public⁴.

2.2. Cosmic shear

The observed cosmic shear angular power spectrum, C_ϵϵ(ℓ), measures a combination of the distortions arising from weak gravitational lensing by large-scale structures (labelled with a subscript ‘G’) with a low-level contaminating astrophysical signal arising from the intrinsic alignment of galaxies with the large-scale structures within which they are embedded (labelled with a subscript ‘I’). These contributions can be separated as

$$\begin{array}{c}\hfill C_{ϵϵ}^{(ij)}(\ell )=C_{\mathrm{GG}}^{(ij)}(\ell )+C_{\mathrm{GI}}^{(ij)}(\ell )+C_{\mathrm{IG}}^{(ij)}(\ell )+C_{\mathrm{II}}^{(ij)}(\ell ),\end{array}$$(1)

where the indices i and j indicate cross-correlations between the five tomographic source samples. The theoretical power spectra are given by Limber-approximated projections with

$$\begin{array}{c}\hfill C_{\mathit{ab}}^{(ij)}(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{hor}}}\mathrm{d}\chi \frac{W_a^{(i)}(\chi )W_b^{(j)}(\chi )}{f_{\mathrm{K}}^2(\chi )}{P}_{\mathrm{m},\mathrm{nl}}\left(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},z(\chi )\right),}\end{array}$$(2)

where a, b ∈ {I, G}, f_K(χ) is the comoving angular diameter distance and χ is the comoving radial distance which runs out to the horizon, χ_hor. The weight functions, W(χ), encode information about how the signal scales with the KiDS-1000 survey depth (see Eqs. (15) and (16) of Joachimi et al. 2021). In the cases of power spectra that include intrinsic ‘I’ terms, the weight function also encodes the intrinsic galaxy alignment model, which we take to be the ‘NLA’ model from Bridle & King (2007). For Stage III surveys like KiDS-1000, this model has been shown to be sufficiently flexible, capturing the likely more complex underlying intrinsic alignment model, without biasing cosmological parameters (Fortuna et al. 2021, see the discussion in Appendix D). The cosmological information for cosmic shear power spectrum is contained in both the geometric weight functions, W(χ), and in the evolution and shape of the non-linear matter power spectrum, P_m, nl(k, z), which we model using the halo formalism⁵ of Mead et al. (2015, 2016). Weak lensing is therefore a very valuable cosmological probe, as it is sensitive to changes in both the distance-redshift relation and to the growth of structures.

We estimate the cosmic shear angular power spectrum through a linear transformation of the real-space two-point shear correlation function (Schneider et al. 2002). This approach circumvents the challenge of accurately determining the survey mask for a direct power spectrum estimate. Joachimi et al. (2021) detail the apodisation advances that we have adopted for the transformation, in addition to the modelling that we use to account for the minor differences between the theoretical expectation of the true angular power spectrum in Eq. (1) and the measured ‘band powers’.

Figure 1 presents the Asgari et al. (2021a) KiDS-1000 cosmic shear power spectra for the auto- and cross-correlated tomographic bins. Here we have constructed both E-mode (upper left) and B-mode (lower right) band powers in order to isolate any non-lensing B-mode distortions (see Eqs. (17)–(21) of Joachimi et al. 2021). As expected from the analysis of Giblin et al. (2021), the measured B-modes are found to be consistent with zero⁶. The measured E-modes can be compared to the theoretical expectation from Eq. (1), given the best-fit set of cosmological parameters from our multi-probe analysis in Sect. 3.

Fig. 1. KiDS-1000 cosmic shear power spectra: tomographic band powers comparing the E-modes (upper left block) with the best-fit cosmological model from our combined multi-probe analysis. The tomographic bin combination is indicated in the upper right corner of each sub-panel. The null-test B-modes (lower right block – note the reduced ordinate scale), are consistent with zero for both the full data vector and each bin combination individually. The errors are estimated analytically (Joachimi et al. 2021). See Sect. 3 for a discussion on the goodness-of-fit.

2.3. Anisotropic galaxy clustering

Galaxy clustering observations probe the 3D non-linear galaxy-galaxy power spectrum and we follow Sánchez et al. (2017) in modelling this quantity based on a perturbation theory approach with

$$\begin{array}{c}\hfill P_{\mathrm{gg}}(k,z)={\displaystyle \sum _{\alpha ,\beta }\alpha \beta P_{\alpha \beta }(k,z)+b_1{\gamma }_3^{-}{P}_{b_1{\gamma }_3^{-}}(k,z)+P_{\mathrm{noise}}(k,z).}\end{array}$$(3)

Here α, β ∈ [b₁,b₂,γ₂], introduce the linear and quadratic bias parameters b₁ and b₂, in addition to the non-local bias parameters γ₂ and ${\gamma }_3^{-}$. Each power spectrum term on the right hand side of the equation is given by different convolutions of the linear matter power spectrum in Appendix A of Sánchez et al. (2017). In the case of an effective linear galaxy bias model (see for example van Uitert et al. 2018; Abbott et al. 2018) only the b₁ bias parameter is considered to be non-zero and Eq. (3) reduces to $P_{\mathrm{gg}}(k,z)=b_1^2{P}_{\mathrm{m},\mathrm{nl}}^{\mathrm{pert}}(k,z)$, where $P_{\mathrm{m},\mathrm{nl}}^{\mathrm{pert}}(k,z)$ is the perturbation theory estimate of the non-linear matter power spectrum which is accurate at the two percent level to k ≲ 0.3 h⁻¹ Mpc (Sánchez et al. 2017).

Sánchez et al. (2017) present the anisotropic redshift-space correlation function of galaxy clustering with the galaxy pairs separated into three ‘wedges’ equidistant⁷ in μ, where μ is the cosine of the angle between the line of sight and the line connecting the galaxy pairs. As such the 3D correlation function is measured for pairs that are either mainly transverse to the line of sight, mainly parallel to the line of sight, or placed into an intermediate sample between these two cases. The redshift-space correlation function ξ_gg(s,μ,z), where s is the co-moving galaxy-pair separation, is given by

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{gg}}\left(s,\mu ,z\right)=& {\displaystyle \sum _{l=0}^2{\mathrm{L}}_{2l}(\mu )\frac{{(-1)}^l(4l+1)}{{(2\pi )}^2}{\int }_0^{\infty }\mathrm{d}k{k}^2{j}_{2l}(ks)}\hfill \\ \hfill & \times {\displaystyle {\int }_{-1}^1\mathrm{d}{\mu }_1{L}_{2l}({\mu }_1)P_{\mathrm{gg},s}(k,{\mu }_1,z),}\hfill \end{array}$$(4)

where L_i denotes the Legendre polynomial of degree i, j_i is the spherical Bessel function of order i, and P_gg, s(k, μ, z) is the 3D redshift-space power spectrum that includes the non-linear real-space power spectrum, Eq. (3), and the galaxy-velocity and velocity-velocity power spectrum (see Sánchez et al. 2017, for details, including how the Alcock-Paczynski distortions are accounted for in the modelling). The same model and ‘wedge’ approach was adopted in the Fourier-space anisotropic galaxy clustering analysis of Grieb et al. (2017), finding consistent results.

Figure 2 presents the Sánchez et al. (2017) BOSS-DR12 anisotropic clustering correlation functions in three wedges and two redshift slices⁸, with 0.2 < z ≤ 0.5, and 0.5 < z ≤ 0.75, for the scales used in this analysis with 20 < s < 160 h⁻¹ Mpc. The measured correlation functions can be compared to the theoretical expectation⁹ given the best-fit set of cosmological parameters from our joint multi-probe analysis in Sect. 3.

Fig. 2. BOSS-DR12 anisotropic clustering from Sánchez et al. (2017): the transverse (pink), intermediate (blue) and parallel (black) clustering wedges in two redshift bins, compared with the best-fit cosmological model from our combined multi-probe analysis. The errors, estimated from mock BOSS catalogues (Kitaura et al. 2016), are highly correlated, particularly at large scales (see Fig. A.1).

2.4. Galaxy-galaxy lensing

The observed galaxy-galaxy lensing angular power spectra, C_nϵ(ℓ), measure a combination of weak lensing distortions around foreground galaxies (labelled with a subscript ‘gG’) with a low-level intrinsic alignment signal arising from the fraction of the source galaxy population that reside physically close to the lenses (labelled with a subscript ‘gI’). We also consider the low-level lensing-induced magnification bias (labelled with a subscript ‘mG’). These three contributions can be separated as

$$\begin{array}{c}\hfill C_{nϵ}^{(ij)}(\ell )=C_{\mathrm{gG}}^{(ij)}(\ell )+C_{\mathrm{gI}}^{(ij)}(\ell )+C_{\mathrm{mG}}^{(ij)}(\ell ),\end{array}$$(5)

where the index i indicates the two lens galaxy samples and the index j indicates the five tomographic source samples. The theoretical power spectra are given by Limber-approximated projections following Eq. (2), with two key differences. The first is that the lens weight function is replaced by the redshift distribution of the lenses. The second is that the non-linear matter power spectrum is replaced with the non-linear cross power spectrum between the galaxy and matter distribution P_gm(k, z) (see Eqs. (28) and (29) of Joachimi et al. 2021, for the full expressions). For the magnification bias power spectrum, C_mG(ℓ), we refer the reader to Appendix B in Joachimi et al. (2021). We find that the inclusion or exclusion of this term has a negligible impact on our cosmological constraints, but we retain it nevertheless.

We adopt the non-linear galaxy bias model from Sánchez et al. (2017), and approximate the non-linear cross power spectrum as

$$\begin{array}{cc}\hfill P_{\mathrm{gm}}(k,z)=& b_1{P}_{\mathrm{m},\mathrm{nl}}(k,z)+\{b_2{\mathcal{F}}_{b_2}(k)-{\gamma }_2{\mathcal{F}}_{{\gamma }_2}(k)\hfill \\ \hfill & -{\gamma }_3^{-}{\mathcal{F}}_{{\gamma }_3^{-}}(k)\}{P}_{\mathrm{m},\mathrm{lin}}^2(k,z).\hfill \end{array}$$(6)

Here P_m, nl(k, z) is the non-linear matter power spectrum modelled using Mead et al. (2016), in contrast to the less accurate perturbation theory estimate $P_{\mathrm{m},\mathrm{nl}}^{\mathrm{pert}}(k,z)$ used in Eq. (3). The logarithm of the functions ℱ_α(k) are second-order polynomial fits that we use to model $P_{\alpha }(k,z_{\text{ref}})/P_{\text{m},\text{lin}}^2(k,z_{\text{ref}})$, the ratio between the different bias terms in the full perturbation theory model in Eq. (3), and the square of the linear matter power spectrum. This approach permits a reasonable extrapolation of the Sánchez et al. (2017) perturbation model into the non-linear regime beyond k = 0.3 h⁻¹ Mpc. This is necessary in order to carry out the redshift-weighted projection of the 3D model to estimate the 2D galaxy-galaxy lensing observable, Eq. (2). No matter which ℓ-scales we restrict our analysis to, high-k scales will contribute to all angular scales at some level (Joachimi et al. 2021; Asgari et al. 2021b). This approach also decreases the compute time of the galaxy-galaxy lensing likelihood evaluations, by several orders of magnitude, in comparison to the direct perturbative calculation.

Figure 3 presents the KiDS-1000 galaxy-galaxy lensing power spectra, around lenses from the BOSS and 2dFLenS surveys (see Blake et al. 2020, for the real-space KiDS-1000 galaxy-galaxy lensing measurements for BOSS and 2dFLenS separately). Each panel presents the cross-correlation between each of the five different tomographic source bins, denoted ‘S’, with the two different lens bins, denoted ‘L’ (see Table A.1 for details). Here we have constructed both E-mode (left) and B-mode (right) band powers in order to isolate any non-lensing B-mode distortions. As expected from the analysis of Giblin et al. (2021), the measured B-modes are found to be consistent with zero. The measured E-modes can be compared to the theoretical expectation given the best-fit set of cosmological parameters from our joint multi-probe analysis in Sect. 3 in the non-shaded regions.

Fig. 3. KiDS-1000 galaxy-galaxy lensing power spectra: tomographic band powers comparing the E-modes (left block) with the best-fit cosmological model from our combined multi-probe analysis. The tomographic bin combination of BOSS and 2dFLenS lenses (L) with KiDS-1000 sources (S), is indicated in the upper right corner of each sub-panel. Data within grey regions are not included in the cosmological analysis. The null-test B-modes (right block – note the reduced ordinate scale), are consistent with zero for both the full data vector and each bin combination individually. The errors are estimated analytically (Joachimi et al. 2021).

The shaded regions in Fig. 3 are excluded from our analysis for two reasons. For overlapping lens-source bins (L1 with S1 and L2 with S1 to S3), the intrinsic alignment terms C_gI(ℓ) are expected to become significant. This raises the question of the validity of the arguably rudimentary ‘NLA’ intrinsic alignment model when used in combination with a non-linear galaxy bias model (see Blazek et al. 2019, for a self-consistent pertubative approach to both intrinsic alignment and galaxy bias modelling). As these bin combinations carry little cosmological information, we exclude this data from our cosmological inference analysis, using it instead in a redshift-scaling null test of the catalogue in Giblin et al. (2021). For separated lens-source bins we introduce a maximum ℓ-scale beyond which the contributions from scales k > 0.3 h⁻¹ Mpc become significant¹⁰. In this regime uncertainties in the extrapolation of the Sánchez et al. (2017) non-linear galaxy bias model into the non-linear regime (Eq. (6)) may well render the C_nϵ(ℓ) model invalid. The ℓ-limit depends on the redshift of the lens bin. Figure 3 therefore serves as an important illustration of the necessity of improving non-linear galaxy bias and non-linear intrinsic alignment modelling for future studies, in order to fully exploit the cosmological signal contained within the galaxy-galaxy lensing observable.

2.5. Multi-probe covariance

Joachimi et al. (2021) present the multi-probe covariance matrix adopted in this study, verified through an analysis of over 20 000 fast full-sky mock galaxy catalogues derived from lognormal random fields. Given that only 4% of the BOSS footprint overlaps with KiDS-1000, in an initial step we validate the approximation that the BOSS anisotropic galaxy clustering observations are uncorrelated with the cosmic shear and galaxy-galaxy lensing observations. By imposing realistic overlapping BOSS and KiDS-1000 footprints in our mock catalogues, we find that cross-correlation, between the projected BOSS-like angular galaxy correlation function, and the KiDS-1000-like weak lensing signals, is less than ∼5% of the auto-correlation terms along the diagonal of covariance matrix. With such a low cross-correlation, we can safely assume independence between the clustering and lensing observations, allowing us to adopt the Sánchez et al. (2017) covariance matrix¹¹ for the anisotropic galaxy clustering observations, ξ_gg(s, μ, z), setting the clustering-lensing cross-correlation terms to zero.

The covariance of the two weak lensing observations is calculated analytically, combining terms that model pure Gaussian shape noise, survey sampling variance, and the noise-mixing that occurs between these two components, in addition to higher-order terms that account for mode-mixing between the in-survey modes and between the observed in-survey and the unobserved out-of-survey modes (known as super-sample covariance, Takada & Hu 2013). The covariance also includes a contribution to account for our uncertainty on the multiplicative shear calibration correction (Kannawadi et al. 2019). Joachimi et al. (2021) demonstrate that every term in the covariance is important, each dominating in different regions of the covariance with one exception: non-Gaussian variance between the in-survey modes is always sub-dominant. We therefore review the approximations made in these analytical calculations. Whilst the complex KiDS-1000 mask is fully accounted for in the shape-noise terms, and the super-sample terms, for all other terms it is assumed that the scales we measure are much smaller than any large-scale features in the survey footprint. Furthermore, the survey is assumed to be homogeneous in its depth, which is invalid for any ground-based survey where the survey depth becomes a sensitive function of the observing conditions (Heydenreich et al. 2020). With mock catalogues, we have the freedom to impose complex masks and variable depth to quantify the impact of these effects on the derived covariance, finding differences typically ≲10%, with a maximum difference of ∼20%. The majority of the differences were found to be driven by the mix-term between the Gaussian shape noise and the Gaussian sampling variance. Through a mock multi-probe data vector inference analysis, Joachimi et al. (2021) demonstrate that these differences in the covariance are not expected to lead to any systematic bias in the recovery of the KiDS-1000 cosmological parameters, nor to any significant differences in the confidence regions of the recovered parameters. We therefore adopt an analytical covariance in our analysis, shown in Fig. A.1, and refer the reader to Joachimi et al. (2021) for further details, where their Sect. 4 presents the mocks, Sect. 5 and Appendix E presents the analytical covariance model, and Appendix D presents detailed comparisons of the mock and analytical covariance.

2.6. Parameter inference methodology

We use the KiDS Cosmology Analysis Pipeline, KCAP¹² built from the COSMOSIS analysis framework of Zuntz et al. (2015), adopting the nested sampling algorithm MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2019). The KCAP bespoke modules include: the BOSS wedges likelihood from Sánchez et al. (2017); the band power cosmic shear and galaxy-galaxy lensing likelihood based on Eqs. (1) and (6); tools to permit correlated priors on nuisance parameters; and tools to sample over the clustering parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, a parameter which is typically only derived. Scripts are also provided to derive the best-fit parameter values at the maximum multivariate posterior, denoted MAP (maximum a posteriori), and an associated credible region given by the projected joint highest posterior density region, which we denote PJ-HPD (Joachimi et al. 2021). This concise list of new modules reflect the primary updates in the KiDS-1000 parameter inference methodology compared to previous KiDS analyses, which we discuss in more detail below.

Our 3 × 2 pt model has 20 free parameters, with five to describe flat ΛCDM in addition to fifteen nuisance parameters. Eight of these nuisance parameters describe the galaxy bias model, with four in each lens redshift bin. The remaining seven allow us to marginalise over our uncertainty on the impact of baryon feedback (one parameter), intrinsic galaxy alignment (one parameter), and the mean of the source redshift distribution in each tomographic bin (five correlated parameters). The priors adopted for each parameter are listed in Appendix B.

Adopted priors are usually survey-specific, with the intention to be uninformative on the parameter that lensing studies are most sensitive to, S₈. Different prior choices, particularly on the amplitude of the primordial power spectrum of scalar density fluctuations, A_s, have however been shown to lead to non-negligible changes in the derived S₈ parameter (Joudaki et al. 2017, 2020; Chang et al. 2019; Asgari et al. 2020). Joachimi et al. (2021) show that even with wide priors on A_s, the sampling region in the σ₈ − Ω_m plane is significantly truncated at low values of σ₈ and Ω_m, with the potential to introduce a subtle bias towards low values of σ₈. In this analysis, we address this important issue of implicit informative priors by sampling directly in S₈. By adopting a very wide S₈ prior, our constraints on S₈ are therefore not impacted by our choice of prior. We note, however, that this approach is expected to lead to a more conservative constraint on S₈, compared to an analysis that adopts a uniform prior¹³ on lnA_s.

We account for the uncertainty in our source redshift distributions using nuisance parameters, ${\delta }_z^i$, which modify the mean redshift of each tomographic bin i. By analysing mock KiDS catalogues, Wright et al. (2020a) determined the mean bias per redshift bin, μⁱ, and also the covariance between the different redshift bins, C_δz. This covariance arises from sampling variance in the spectroscopic training sample, which impacts, to some degree, the redshift calibration of all bins. We therefore adopt the multivariate Gaussian prior 𝒩(μ; C_δz), for the vector δ_z (see Sect. 3 of Hildebrandt et al. 2021, for details).

Adopting the Bayesian paradigm for inference, we provide our constraints in the form of a series of samples that describe the full posterior distribution¹⁴. In Sect. 3 we explore this multi-dimensional posterior in the traditional way, visualising the 2D and 1D marginal posterior distributions for a selection of parameters. In cosmological parameter inference it is standard to also report a point estimate of the one-dimensional marginal posterior distribution with an associated 68% credible interval. It is not always stated, however, how these point estimates and intervals are defined.

We provide two different point estimates for our cosmological parameters, with the first reporting the standard maximum of the marginal distribution, along with a credible interval that encompasses 68% of the marginal highest posterior density, which we denote by M-HPD. For the high-dimensional parameter space of a multi-probe weak lensing analysis, we find that this standard marginalised point estimate leads to a value for S₈ that is lower than the maximum of the multivariate joint posterior, with an offset of up to ∼1σ, dependent on which probes are combined (see Sect. 7 of Joachimi et al. 2021). This is not a result of an error in the KCAP inference pipeline. Rather, it is a generic feature of projecting high-dimensional asymmetric distributions into one dimension, prompting the development of an alternative approach to reporting point estimates for cosmological parameters.

Our fiducial S₈ constraints follow this alternative, reporting the parameter value at the maximum of the joint posterior (MAP), along with a 68% credible interval based on the joint, multi-dimensional highest posterior density region, projected onto the marginal posterior of the S₈ parameter (PJ-HPD). Here we step through the posterior MULTINEST samples, ordered by their decreasing posterior density. For each model parameter we determine the extrema within the n highest posterior samples, and the posterior mass contained within the marginal distribution of each parameter, limited by the extrema values. We iterate, increasing the number of samples analysed, n, until the posterior mass reaches the desired 68% level. The PJ-HPD credible interval is then reported as the parameter extrema at this point n in the sample list, and we repeat the process for each model parameter of interest (see Sect. 6.4 of Joachimi et al. 2021, for further details).

We note that the MAP reported by MULTINEST provides a noisy estimate of the true MAP due to the finite number of samples, and we therefore conduct an optimisation step using the 36 samples with the highest posterior values as starting points. We use both the Nelder & Mead (1965) and Powell (1964) optimisation algorithms, as well as two-step optimisation using both algorithms. While the MAP estimates found in this optimisation step increase the posterior probability by a factor of 2 to 4 compared to the MAP estimated from the MULTINEST samples, they exhibit scatter in parameters constrained by the galaxy clustering likelihood. We suspect this is due to numerical noise in the galaxy clustering likelihood which results in many local minima that inhibit the convergence of the optimisation step. This suspicion is strengthened by the fact that the MAP estimates do not exhibit this scatter for probe combinations that exclude the clustering observable. For this reason, we report the median of the MAP estimates, weighted by their posterior probability, for probes that include the galaxy clustering likelihood, since a global optimisation would be computationally prohibitively expensive. For the other probes the reported MAP is given by the parameter set at the maximum posterior found amongst all estimates (see also Muir et al. 2020, who adopt a similar approach).

We note that the presence of offsets between marginal S₈ constraints and those derived from the full multivariate joint posterior highlights how efforts to accurately quantify tension based solely on one-point estimates should be undertaken with some level of caution. Tension can also be assessed in terms of the overlap between the full posterior distributions (see for example Handley & Lemos 2019; Lemos et al. 2020; Raveri & Hu 2019), which we discuss further in Sect. 3.2 and Appendix G.

3. Results

We present our multi-probe constraints on the cosmological parameters of the flat ΛCDM model in Fig. 4, showing the marginalised posterior distributions for matter fluctuation amplitude parameter, σ₈, the matter density parameter, Ω_m, and the dimensionless Hubble parameter, h, where the BOSS galaxy clustering constraints (shown blue), break the σ₈ − Ω_m degeneracy in the KiDS-1000 cosmic shear constraints (shown pink), resulting in tight constraints on σ₈ in the combined 3 × 2 pt analysis (shown red). Reporting the MAP values with PJ-HPD credible intervals for the parameters that we are most sensitive to, we find

$$\begin{array}{cc}& {\sigma }_8=0.{76}_{-0.020}^{+0.025}\hfill \\ \hfill & {\mathrm{\Omega }}_{\mathrm{m}}=0.{305}_{-0.015}^{+0.010}\hfill \\ \hfill & S_8=0.{766}_{-0.014}^{+0.020}.\hfill \end{array}$$(7)

Fig. 4. Marginal multi-probe constraints on the flat ΛCDM cosmological model, for the matter fluctuation amplitude parameter, σ8, the matter density parameter, Ωm, and the dimensionless Hubble parameter, h. The BOSS galaxy clustering constraints (blue), can be compared to the KiDS-1000 cosmic shear constraints (pink), the combined 3 × 2 pt analysis (red), and CMB constraints from Planck Collaboration VI (2020, grey).

Our constraints can be compared to the marginalised posterior distributions from Planck (shown grey in Fig. 4), finding consistency between the marginalised constraints on Ω_m and h, but an offset in σ₈, which we discuss in detail in Sect. 3.2.

Tabulated constraints for the full set of cosmological parameters are presented in Appendix C, quoting our fiducial MAP with PJ-HPD credible intervals along with the marginal posterior mode with M-HPD credible intervals. As discussed in Joachimi et al. (2021), the marginal mode estimate is known to yield systematically low values of S₈ in mock data analyses. This effect can be seen in Fig. 5 which compares the joint posterior constraints (solid) with the marginal posterior constraints (dashed).

Fig. 5. Constraints on the structure growth parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ for different probe combinations: 3 × 2 pt, KiDS-1000 cosmic shear, BOSS galaxy clustering, cosmic shear with galaxy-galaxy lensing (GGL), and cosmic shear with galaxy clustering. Our fiducial and preferred MAP with PJ-HPD credible interval (solid) can be compared to the standard, but shifted, marginal posterior mode with M-HPD credible intervals (dotted). Our results can also be compared to weak lensing measurements from the literature, which typically quote the mean of the marginal posterior mode with tail credible intervals (dashed).

We find good agreement between the different probe combinations and single-probe S₈ constraints, demonstrating internal consistency between the different cosmological probes, in Fig. 5. As forecast by Joachimi et al. (2021), the addition of the galaxy-galaxy lensing observable adds very little constraining power, with similar results found for the full 3 × 2 pt analysis and the combined cosmic shear and clustering analysis. This is primarily a result of the significant full area of BOSS in comparison to the size of the BOSS-KiDS overlap region. The lack of an accurate galaxy bias model on the deeply non-linear scales that weak lensing probes also prohibits the inclusion of large sections of our galaxy-galaxy lensing data vector, shown in Fig. 3. The addition of the galaxy-galaxy lensing does, however, serve to moderately tighten constraints on the amplitude of the intrinsic alignment model A_IA, as seen in Fig. 6.

Fig. 6. Marginalised posterior distributions for an extended set of cosmological parameters covering the matter density parameter, Ωm, the matter fluctuation amplitude parameter, σ8, the structure growth parameter, S8, the spectral index, ns, the dimensionless Hubble parameter, h, the baryon feedback amplitude parameter, Abary, the intrinsic alignment amplitude, AIA, and the linear bias parameters for the low and high BOSS redshift bins, b1. The KiDS-1000 cosmic shear results (pink), can be compared to the BOSS galaxy clustering results (blue), the combination of cosmic shear with BOSS and 2dFLenS galaxy-galaxy lensing (GGL, purple), and the full 3 × 2 pt analysis (red). The combination of cosmic shear with galaxy clustering (orange) is only distinguishable from the 3 × 2 pt result in the Abary and AIA panels. For parameters constrained by the CMB, we also include constraints from Planck Collaboration VI (2020, grey).

Figure 5 also demonstrates the good agreement between our constraints and weak lensing results from the literature, comparing to cosmic shear-only results from the Hyper Suprime-Cam Strategic Programme (HSC, Hikage et al. 2019; Hamana et al. 2020), DES Y1 (Troxel et al. 2018) and an earlier KiDS analysis (KV450, Hildebrandt et al. 2020), in addition to the previous KV450-BOSS ‘2 × 2 pt’ analysis of Tröster et al. (2020b) and the DES Y1 3 × 2 pt analysis from Abbott et al. (2018). We refer the reader to Asgari et al. (2021a) for a discussion and comparison of different cosmic shear results. In Sect. 3.1 we present a more detailed comparison of our results with 3 × 2 pt results in the literature.

Figure 6 displays the marginal posterior distributions for an extended set of cosmological parameters. We find that the allowed range for the linear galaxy bias, b₁, in each redshift bin (lower two rows), is almost halved with the addition of the weak lensing data. This constraint does not arise, however, from the sensitivity of the galaxy-galaxy lensing observable to galaxy bias (shown to be relatively weak in the purple cosmic shear + GGL contours). Instead, in this analysis, it is a result of the degeneracy breaking in the σ₈ − Ω_m plane, tightening constraints on σ₈ which, for galaxy clustering, is degenerate with galaxy bias. The improved constraints on galaxy bias do not, however, fold through to improved constraints on h, which the weak lensing data adds very little information to.

For our primary cosmological parameter, S₈, our constraints are uninformed by our choice of priors. This statement cannot be made for the other ΛCDM parameters, however, as shown in Fig. 6. The most informative prior that we have introduced to our 3 × 2 pt analysis is on the spectral index, n_s. As noted by Tröster et al. (2020b), the BOSS galaxy clustering constraints favour a low value for n_s, where they find n_s = 0.815 ± 0.085. From the Tröster et al. (2020b) sensitivity analysis to the adopted maximum clustering scale we observe that this preference appears to be driven by the amplitude of the large scale clustering signal with s > 100 h⁻¹ Mpc. We note that spurious excess power in this regime could plausibly arise from variations in the stellar density impacting the BOSS galaxy selection function (Ross et al. 2017). Our choice to impose a theoretically motivated informative prior for n_s, as listed in Table B.1, helps to negate this potential systematic effect without degrading the overall goodness-of-fit to the galaxy clustering measurements. Our prior choice is certainly no more informative than the n_s priors that are typically used in weak lensing and clustering analyses (see for example Abbott et al. 2018; eBOSS Collaboration 2020). We recognise, however, that this well-motivated prior choice acts to improve the BOSS-only error on Ω_m by roughly a third, and decrease the BOSS-only best-fitting value for Ω_m and h by ∼0.5σ (see Fig. B.1). With < 10% differences on the constraints on S₈ and h, however, and only a ∼0.1σ difference in the BOSS-only best-fitting value for S₈, which is consistent with the typical variation between different MULTINEST analyses, we conclude that our prior choice does not impact on our primary S₈ constraints (see Appendix B). With the informative or uninformative n_s prior, our constraints on h remain consistent with the Hubble parameter constraints from both Planck Collaboration VI (2020) and Riess et al. (2019).

Figure 7 illustrates the results of a series of sensitivity tests, where we explore how our 3 × 2 pt constraints on S₈ change when: we ignore the impact of baryon feedback (the ‘No baryon’ case) by fixing A_bary = 3.13, corresponding to the non-linear matter power spectrum for a dark-matter only cosmology; we limit the analysis to a linear galaxy bias model, setting all higher-order bias terms in Eq. (3) to zero, as well as restricting the redshift-space distortion model to a Gaussian velocity distribution; and we remove individual tomographic bins from our weak lensing observables. The systematic offset that arises from the use of a linear-bias model highlights the importance of accurate non-linear galaxy bias modelling in 3 × 2 pt analyses. This series of tests is dissected further in Appendix E, and complements the detailed KiDS-1000 internal consistency analysis of Asgari et al. (2021a, Appendix B), which demonstrates that the change seen with the removal of tomographic bin 4 is fully consistent with expected statistical fluctuations.

Fig. 7. 3 × 2 pt constraints on S8 for a series of sensitivity tests; when we ignore the impact of baryon feedback (the ‘No baryon’ case), limit the analysis to a linear galaxy bias model (the ‘No higher order GC’ case), and remove individual tomographic bins from our weak lensing observables.

Table 1 records the goodness-of-fit for each component in our 3 × 2 pt analysis, where we report the χ² value at the maximum posterior, ${\chi }_{\text{MAP}}^2$ (see Sect. 2.6 for a discussion of our optimised MAP-finder). The effective number of degrees of freedom (d.o.f.) does not equate to the standard difference between the total number of data points (Data d.o.f.) and the total number of model parameters (20 in the case of our 3 × 2 pt analysis), as a result of the adopted priors and the non-linear dependencies that exist between the model parameters. For some probe combinations we calculate the effective number of degrees of freedom in the model (Model d.o.f.), using the estimator described in Sect. 6.3 of Joachimi et al. (2021). As this approach is computationally expensive, however, we also estimate the Model d.o.f. following Raveri & Hu (2019), recognising that, for the cases explored in Joachimi et al. (2021), this approach results in a slightly lower model d.o.f.

We find that the goodness-of-fit is excellent for the BOSS galaxy clustering. For all other cases, the goodness-of-fit is certainly acceptable¹⁵, with the probability to exceed the measured χ² given by p ≳ 0.01. We note that the cosmic shear analysis of Asgari et al. (2021a) shows no significant changes in the inferred cosmological parameters when using different two-point statistics which exhibit an excellent goodness-of-fit. As such, we could be subject to an unlucky noise fluctuation that particularly impacts the band power estimator in Eq. (1). Cautiously inspecting Fig. 1, as ‘χ-by-eye’ is particularly dangerous with correlated data points, we nevertheless note a handful of outlying points, for example the low ℓ-scales in the fifth tomographic bin. We also note that Giblin et al. (2021) document a significant but low-level PSF residual systematic in the KiDS-1000 fourth and fifth tomographic bins that was shown to reduce the overall goodness-of-fit in a cosmic shear analysis, but not bias the recovered cosmological parameters (see the discussion in Amara & Réfrégier 2008). Future work to remove these low-level residual distortions is therefore expected to further improve the goodness-of-fit.

3.1. Comparison with weak lensing surveys

Our results are consistent with weak lensing constraints in the literature. We limit our discussion in this section to published 3 × 2 pt analyses, referring the reader to Asgari et al. (2021a) who discuss how the KiDS-1000 cosmic shear results compare with other weak lensing surveys. We note that direct comparisons of cosmological parameters should be approached with some caution, as the priors adopted by different surveys and analyses are often informative (see Sect. 6.1 in Joachimi et al. 2021). Homogenising priors for cosmic shear analyses, for example, has been shown to lead to different conclusions when assessing inter-survey consistency (Chang et al. 2019; Joudaki et al. 2020; Asgari et al. 2020).

Abbott et al. (2018) present the first year 3 × 2 pt DES analysis (DES Y1), finding $S_8=0.{773}_{-0.020}^{+0.026}$, where they report the marginal posterior maximum and the tail credible intervals. This is in excellent agreement with our equivalent result, differing by 0.3σ, with the DES-Y1 error being 40% larger than the KiDS-1000-BOSS 3 × 2 pt results. The inclusion of BOSS to our 3 × 2 pt analysis results in tight constraints on Ω_m. This leads to joint KiDS-1000-BOSS constraints on ${\sigma }_8=0.{760}_{-0.023}^{+0.021}$ that are more than twice as constraining compared to the DES Y1-alone 3 × 2 pt analysis which found ${\sigma }_8=0.{817}_{-0.056}^{+0.045}$, as shown in Fig. 8. This comparison serves to highlight the additional power that can be extracted through the combination of spectroscopic and photometric surveys, and the promising future for the planned overlap between the Dark Energy Spectroscopic Instrument survey (DESI Collaboration 2016) and the 4 m Multi-Object Spectroscopic Telescope (4MOST, Richard et al. 2019), with Euclid and the Vera C. Rubin Observatory Legacy Survey of Space and Time (Laureijs et al. 2011; LSST Science Collaboration 2009), in addition to the nearer-term ∼1400 deg² of overlap between BOSS and HSC (Aihara et al. 2019).

Fig. 8. Marginalised posterior distribution in the σ8 − Ωm plane, comparing the 3 × 2 pt analyses from KiDS-1000 with BOSS and 2dFLenS, with the 3 × 2 pt analysis from DES Y1 (Abbott et al. 2018), and the CMB constraints from Planck Collaboration VI (2020). The KiDS-1000 3 × 2 pt result can also be compared to our previous KV450-BOSS analysis from Tröster et al. (2020b).

van Uitert et al. (2018) and Joudaki et al. (2018) present 3 × 2 pt analyses for the second KiDS weak lensing release (KiDS-450), finding, respectively, $S_8=0.{800}_{-0.027}^{+0.029}$ (KiDS with GAMA) and S₈ = 0.742 ± 0.035 (KiDS with BOSS and 2dFLenS limited to the overlap region). Both results are consistent with our KiDS-1000 results, noting that the increase in our S₈ constraining power, by a factor of ∼2 in this analysis, is driven by increases in both the KiDS survey area, and the analysed BOSS survey area.

The impact of doubling the KiDS area can be seen by comparing to Tröster et al. (2020b), in Fig. 8, who present a joint cosmic shear and galaxy clustering analysis of the KV450 KiDS release with the full BOSS area, finding S₈ = 0.728 ± 0.026. The ∼40% improvement in constraining power is consistent with expectations from the increased survey area, but a straightforward area-scaling comparison is inappropriate given that KiDS-1000 features improvements in the accuracy of the shear and photometric redshift calibrations, albeit at the expense of a decrease in the effective number density (see Giblin et al. 2021; Hildebrandt et al. 2021, for details).

The offset in S₈ between the KiDS-1000-BOSS and KV450-BOSS S₈ constraints reflects a number of differences between the two analyses. First, as the S₈ constraints from the 3 × 2 pt analysis are primarily driven by KiDS (see Fig. 6), we expect a reasonable statistical fluctuation in this parameter given the sampling variance arising from the significant increase in the KiDS survey area. Using a simple model analysis in Appendix F, we conclude that we should expect differences, on average, of |ΔS₈| = 0.016, and as such the increase that we find in S₈ between KV450 and KiDS-1000 is consistent with the expectation from simple statistical fluctuations. BOSS primarily constrains Ω_m which is impacted by the choice of prior on n_s. The wider n_s prior adopted in Tröster et al. (2020b), favours a slightly higher but less well-constrained value for Ω_m, leading to a slightly lower but less well-constrained value for σ₈, when combined with cosmic shear (see Appendix B). If we had also chosen an uninformative prior on n_s for our KiDS-1000-BOSS analysis, a decision that we cannot revise post unblinding, this would have likely served to exacerbate any tension with the Planck CMB constraints.

3.2. Comparison with Planck

In our KiDS-1000-BOSS 3 × 2 pt analysis we find good agreement with Planck for the matter density parameter, Ω_m, and the Hubble parameter, h, (see Fig. 4). The amplitude of matter fluctuations, σ₈, that we infer from the clustering of galaxies within, and lensing by, the large-scale structure of the low-redshift Universe is lower, however, than that inferred by Planck¹⁶ from the CMB.

To quantify the level of discrepancy in the amplitude of matter fluctuations, we first concentrate on the parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ as it is tightly constrained and only exhibits negligible degeneracies, if at all, with the other cosmological parameters, Ω_m, h, and n_s, as illustrated in Fig. 6. Comparing the reported marginal S₈ constraints, we find S₈ to be 8.3 ± 2.6% lower than the CMB constraint from Planck Collaboration VI (2020).

We define the widely used S₈-difference measure

$$\begin{array}{c}\hfill \tau =\frac{|{\overline{S_8}}^{3\times 2pt}-{\overline{S_8}}^{Planck}|}{\sqrt{\mathrm{Var}\left[S_8^{3\times 2\mathrm{pt}}\right]+\mathrm{Var}\left[S_8^{\mathit{Planck}}\right]}},\end{array}$$(8)

where $\overline{S_8}$ and Var[S₈] denote the means and variances of the Planck and 3 × 2 pt S₈ posterior distributions. If both distributions are Gaussian, τ can be used to measure how likely it is that the mean of the difference between the distributions is consistent with zero.

Comparing the S₈ posterior distributions between our 3 × 2 pt analysis and the Planck plik_lite_TTTEEE + lowl + lowE likelihood, we find τ = 3.1, meaning there is a 3.1σ difference between the KiDS-1000 and Planck constraints. Adopting two tension measures that do not assume Gaussianity of the marginal posterior distributions, both the ‘Hellinger’ distance and the distribution of the S₈ parameter shifts indicate 3.1σ difference between our 3 × 2 pt analysis and Planck (see Appendices G.1 and G.2 for details). Our result thus continues the general trend of low-redshift probes preferring low amplitudes of matter fluctuations¹⁷ (Heymans et al. 2013; Alam et al. 2017; Leauthaud et al. 2017; Abbott et al. 2018, 2020; Hikage et al. 2019; Bocquet et al. 2019; Lange et al. 2019; Palanque-Delabrouille et al. 2020; Wright et al. 2020b; Singh et al. 2020). In these cases the reported low S₈, or σ₈, constraints are formally statistically consistent with Planck, and well below the detection of any anomalies at the 5σ-level. Considering, however, the ∼3σ difference that we have reported, and the overall trend in the literature, we would argue that we are reaching an uncomfortable point when it comes to regarding the S₈ offset as a simple statistical fluke.

Sanchez (2020) pointed out that comparing S₈ between different experiments can be misleading due to the implicit dependence of σ₈ on h. Besides the intrinsic dependence of the amplitude of the matter power spectrum on h, measurements of σ₈ also depend on h through the value of 8 h⁻¹ Mpc, the radius of the sphere within which the matter fluctuations are measured. In this way, constraints on σ₈ derived from data sets with different posterior distributions on h represent the average of σ(R) over different ranges of scales. We therefore also consider S₁₂ = σ₁₂(Ω_m h²/0.14)^0.4 (Sanchez 2020), where ${\sigma }_{12}^2$ is the variance of the linear matter field at redshift zero in spheres of radius 12 Mpc. We find $S_{12}=0.{754}_{-0.018}^{+0.015}$, with the value inferred by Planck being $S_{12}=0.{817}_{-0.015}^{+0.011}$, a difference of τ = 3.0, in agreement with the S₈ results.

In light of the large parameter spaces that are being considered, we recognise that focussing on a single parameter can paint a simplistic picture of the agreement, or disagreement, between probes. On a fundamental level, the question we wish to answer is whether a single model of the Universe can describe both the CMB as well as the low-redshift large-scale structure of the Universe. Within our Bayesian inference framework, the Bayes factor provides a natural approach to model selection. The two models under consideration are

• M₁:

  Both our 3 × 2 pt data and Planck’s measurements of the CMB are described by a single flat ΛCDM cosmology.

• M₂:

  The two data sets are described by different cosmologies for the low- and high- redshift Universe, respectively.

The Bayes factor is then

$$\begin{array}{c}\hfill R=\frac{P(\mathit{d}|{\mathrm{M}}_1)P({\mathrm{M}}_1)}{P(\mathit{d}|{\mathrm{M}}_2)P(M_2)},\end{array}$$(9)

where P(d|M_i) is the probability of the data d under model M_i – the Bayesian evidence.

We assume the model priors P(M₁) and P(M₂) to be equal, that is, we make no a-priori assumption on the likelihood of M₁ or M₂. We use anesthetic¹⁸ (Handley 2019) to compute R and find lnR = 3.1 ± 0.3, which can be interpreted as odds of 23 ± 6 in favour of model M₁, and consistency between our 3 × 2 pt measurement and Planck. Given the dependence of R on the parameter priors (Handley & Lemos 2019), we consider this result with some caution and review a series of alternative metrics that seek to quantify the tension between the full KiDS-1000 and Planck multi-dimensional cosmological parameter constraints. These metrics are summarised in Table 2, and below, with further details provided in Appendix G.

Handley & Lemos (2019) propose the ‘suspiciousness’ statistic S that is based on the Bayes factor, R, but hardened against prior dependences. We find that the probability of observing our measured suspiciousness statistic is 0.08 ± 0.02, which corresponds to a KiDS-Planck tension at the level of 1.8 ± 0.1σ (see Appendix G.3 for details).

Raveri & Hu (2019) introduce a number of metrics to quantify the consistency of data sets. We consider their Q_DMAP metric, which explores the change in the goodness-of-fit when two data sets are combined. A reduction in the goodness-of-fit can then be translated into a tension metric, finding a KiDS-Planck tension at the level of 2.1σ (see Appendix G.4 for details). Their Q_UDM metric generalises the notion of parameter differences between posteriors to multiple dimensions, finding KiDS-Planck tension at the level of 1.9σ (see Appendix G.5 for details).

Comparing tension metrics, we find a fairly consistent picture. Reviewing tension in terms of the single parameter that our 3 × 2 pt analysis is most sensitive to, S₈, leads to a ∼3σ tension with Planck. Including additional parameters into the tension analysis, parameters which KiDS is mainly insensitive to, serves to effectively dilute the tension, reducing the measure to the ∼2σ level.

4. Conclusions

In this analysis we have presented constraints on the flat ΛCDM cosmological model by combining observations of gravitational lensing and galaxy clustering to directly probe the evolution and distribution of the large-scale structures in the Universe. Our survey of the z ≲ 1 low-redshift Universe finds a matter distribution that is less clustered, compared to predictions from the best-fitting ΛCDM model to early-Universe CMB observations (Planck Collaboration VI 2020). This tendency for low-redshift probes to favour a smoother matter distribution compared to the CMB expectation has persisted since the first large-scale weak lensing survey (CFHTLenS, Heymans et al. 2013), but the significance of this effect has always been tantalisingly around, or below, the ∼3σ level. It is therefore unclear if these differences are merely a statistical fluctuation, unaccounted for systematic errors, or a sign of interesting new physics.

Our new result does not lead to a resolution in the matter of statistical fluctuations, finding a 3.1σ offset in the structure growth parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ with $S_8=0.{766}_{-0.014}^{+0.020}$. Comparing the marginal S₈ constraints, we find S₈ to be 8.3 ± 2.6% lower than the CMB constraint from Planck Collaboration VI (2020). For a series of ‘tension’ metrics that quantify differences in terms of the full posterior distributions, we find that the KiDS-1000 and Planck results agree at the ∼2σ level. Through our series of image simulation analyses (Kannawadi et al. 2019), catalogue null-tests (Giblin et al. 2021), variable depth mock galaxy survey analyses (Joachimi et al. 2021), optical-to-near-infrared photometric-spectroscopic redshift calibration, validated with mocks (Wright et al. 2020a; van den Busch et al. 2020; Hildebrandt et al. 2021), internal consistency tests (Asgari et al. 2021a, Fig. 6 and Appendix E), and marginalisation over a series of nuisance parameters that encompass our theoretical and calibration uncertainties (Appendix B), we argue that we have, however, addressed the question of 3 × 2 pt systematic errors, robustly assessing and accounting for all sources of systematics that are known about in the literature.

The KiDS-1000 cosmic shear constraints are highly complementary to the BOSS galaxy clustering constraints, leading to tight constraints in our joint 3 × 2 pt analysis that are more than twice as constraining for the matter fluctuation amplitude parameter, ${\sigma }_8=0.{760}_{-0.023}^{+0.021}$, compared to previous 3 × 2 pt analyses. In the future, analysis of the clustering and galaxy-galaxy lensing of photometric samples with very accurate photometric redshifts (see for example Vakili et al. 2019), presents an opportunity for a future alternative KiDS-only 3 × 2 pt photometric analysis, similar to the approach taken in Abbott et al. (2018).

In the next few years, two weak lensing surveys will see first light, with the launch of the Euclid satellite and the opening of the Vera C. Rubin Observatory. These observatories will build the first two ‘full-sky’ weak lensing surveys, which are highly complementary in terms of their differing strengths in depth and spatial resolution¹⁹. Combined with complementary overlapping redshift spectroscopy from DESI, 4MOST and Euclid, the multi-probe weak lensing and spectroscopic galaxy clustering methodology, which we have implemented in this analysis, provides a promising route forward for these next generation surveys. We view this 3 × 2 pt approach as just the start of the story, however, looking forward to a future combined analysis of weak lensing and galaxy clustering with both photometric and spectroscopic lenses, a combination which we call a ‘6 × 2 pt’ approach (Bernstein 2009). This would allow for the optimal combination of information from the clustering cross-correlation of spectroscopic and photometric galaxies (Newman 2008), an observable that we currently only use as an independent tool to validate our photometric redshift calibration (Hildebrandt et al. 2021). Developments in the area of highly non-linear galaxy bias, baryon feedback and intrinsic alignment modelling, along with a sufficiently flexible but tractable redshift distribution model and an accurate ‘6 × 2 pt’ covariance estimate, will all be required in order to realise this long-term goal. The effort will, however, be worthwhile allowing for the implementation of arguably the most robust methodology available to mitigate systematic errors, whilst simultaneously enhancing cosmological parameter constraints.

The ESO-KiDS public survey completed observations in July 2019, spanning 1350 deg². We therefore look forward to the fifth and final KiDS data release, ‘KiDS-Legacy’, along with new results from the concurrent ‘Stage-III’ surveys, DES and HSC, whilst the community prepares for the next exciting chapter of ‘full-sky’ weak lensing surveys.

Acknowledgments

We are completely indebted to Eric Tittley at the IfA for going well beyond the call of duty to save the KiDS-1000 data products after an explosion in our server room destroyed the RAID. We thank our anonymous referee for their positive feedback and useful comments, in addition to our external blinder Matthias Bartelmann who revealed the key for which of the three catalogues analysed was the true unblinded catalogue on the 9th July 2020, right at the end of the KiDS-1000 study which was submitted to A&A on the 30th July 2020. We also wish to thank the Vera C. Rubin Observatory LSST-DESC Software Review Policy Committee (Camille Avestruz, Matt Becker, Celine Combet, Mike Jarvis, David Kirkby, Joe Zuntz with CH) for their draft Software Policy document which we followed, to the best of our abilities, during the KiDS-1000 project. Following this draft policy, the software used to carry out the various analyses presented in this paper is open source at github.com/KiDS-WL/Cat_to_Obs_K1000_P1 and github.com/KiDS-WL/KCAP. The figures in this work were created with matplotlib (Hunter 2007) and getdist (Lewis 2019), making use of the numpy (Oliphant 2006) and scipy (Jones et al. 2001) software packages. We also made extensive use of the TREECORR and COSMOSIS software packages and thank Mike Jarvis and Joe Zuntz for their continuing enhancements and clear documentation. This project has received significant funding from the European Union’s Horizon 2020 research and innovation programme. We thank and acknowledge support from: the European Research Council under grant agreement No. 647112 (CH, TT, MA, CL and BG), No. 770935 (HHi, AHW, JLvdB and AD) and No. 693024 (SJ) in addition to the Marie Skłodowska-Curie grant agreement No. 797794 (TT). We also acknowledge support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research (CH, FK); the Deutsche Forschungsgemeinschaft Heisenberg grant Hi 1495/5-1, (HHi); the Netherlands Organisation for Scientific Research Vici grant 639.043.512 (AK, HHo) and grant 621.016.402 (JdJ); the Alexander von Humboldt Foundation (KK); the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12, and the Polish National Science Center through grants no. 2018/30/E/ST9/00698 and 2018/31/G/ST9/03388 (MB); the Royal Society through an Enhancement Award RGF/EA/181006 (BG); the Australian Research Council grants DP160102235 and CE17010013 (KG); the Beecroft Trust (SJ); STFC grant ST/N000919/1 (LM); the Netherlands Research School for Astronomy and Target (GVK); the NSFC of China under grant 11973070, the Shanghai Committee of Science and Technology grant No. 19ZR1466600, and the Key Research Programme of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013 (HYS). The results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). Data processing for VIKING has been contributed by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh. The BOSS-related results in this paper have been made possible thanks to SDSS-III. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the US Department of Energy Office of Science. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The 2dFLenS-related results are based on data acquired through the Australian Astronomical Observatory, under programme A/2014B/008. It would not have been possible without the dedicated work of the staff of the AAO in the development and support of the 2dF-AAOmega system, and the running of the AAT. Author contributions. All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (CH & TT) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have either made a significant contribution to the data products, or to the scientific analysis.

References

Appendix A: Galaxy properties and the 3  ×  2 pt covariance

In this appendix we tabulate the properties of the KiDS-1000 tomographic source samples, along with the properties of the BOSS and 2dFLenS lens samples, in Table A.1. We list the spectroscopic redshift selection for the lenses (z_min < z_s ≤ z_max), and the photometric redshift selection for the sources (z_min < z_B ≤ z_max), along with the mean redshift of each sample. For the source sample the true redshift distributions are estimated in Hildebrandt et al. (2021), using the SOM methodology from Wright et al. (2020a). The shear calibration correction, m, which can also be referred to in the literature as the responsivity, R = 1 + m, is listed for each source bin (Kannawadi et al. 2019). The effective number density of lenses and sources defines the number of galaxies per square arcminute in the case of unit weights and, for the sources, unit responsivity (see Eqs. (C.11) and (C.13) in Joachimi et al. 2021). We also list the effective ellipticity dispersion, σ_ϵ, i, per ellipticity component, i, for each the weighted and calibrated source galaxy samples (Eq. (C.8) in Joachimi et al. 2021).

Figure A.1 displays the correlation coefficients of the 3 × 2 pt covariance matrix for the three observables; cosmic shear E-mode power spectra, 𝒞_E, galaxy-galaxy lensing E-mode power spectra, 𝒞_nϵ, and the anisotropic galaxy clustering in low and high redshift bins, ξ_gg (see Sect. 2 for details). The cross-correlation between the two lensing and the clustering observables is set to zero, as mock data analyses showed these correlations to be negligible for the KiDS and BOSS footprints (Joachimi et al. 2021). The Fourier-space lensing observables are shown to be significantly less correlated between ℓ-scales, in comparison to the physical-scale clustering observables.

Fig. A.1. Correlation coefficients of the 3 × 2 pt covariance matrix for the cosmic shear, 𝒞E, galaxy-galaxy lensing, 𝒞nϵ and galaxy clustering observables, ξgg (see Sect. 2 for details). Here the band powers 𝒞, are related to the angular power spectrum C, in Eqs. (1) and (3), as 𝒞 = ℓ2C/(2π).

Appendix B: Parameter priors

In this appendix we tabulate the adopted KiDS-1000 priors and sampling parameters in Table B.1. The uniform prior on the dimensionless Hubble constant, h, reflects distance-ladder ±5σ constraints from Riess et al. (2016), which encompasses the value of h favoured by Planck Collaboration VI (2020). The uniform prior on the baryon density, ω_b = Ω_b h², reflects big bang nucleosynthesis ±5σ constraints from Olive & Particle Data Group (2014). The uniform prior on the CDM density, ω_c = Ω_c h², reflects Supernova Type Ia ±5σ constraints on Ω_m from Scolnic et al. (2018) combined with the most extreme allowed values of h and ω_b, given their priors.

As discussed in Sect. 2.6 we choose to sample with an uninformative uniform prior on S₈ to avoid implicit informative priors from a uniform prior on the primordial power spectrum amplitude A_s. We choose a fixed model for the properties of neutrinos, adopting the normal hierarchy at the minimum sum of masses, Σm_ν = 0.06 eV, following Planck Collaboration VI (2020). We consider extended models, including variations on neutrino mass in Tröster et al. (2020a).

The uniform prior on the scalar spectral index, n_s, reflects a restriction in our likelihood implementation, where the Sánchez et al. (2017) galaxy clustering likelihood becomes prohibitively slow for n_s > 1.1. With the upper limit of the top-hat prior fixed by this computational limitation, we choose to symmetrise the prior around the theoretical expectation of n_s = 0.97. In Fig. B.1 we demonstrate the impact of this informative prior on the BOSS-only galaxy clustering constraints, highlighting how informative this choice of prior is. We argue that adopting an informative prior is justified however, given our theoretical prior knowledge of the Harrison-Zel’dovich spectrum. Figure B.1 also helps to illustrate that had we chosen an uninformative prior on n_s for our KiDS-1000-BOSS analysis, a decision taken more than a year before unblinding our analysis, this would have likely served to exacerbate any tension with the Planck CMB constraints.

Fig. B.1. Impact of the ns prior: Comparing marginalised posterior distributions for the BOSS galaxy clustering analysis for our fiducial analysis (blue) with the constraints when adopting an uninformative prior on ns (orange). Opening the parameter space to arguably unphysical values of ns favours higher values and weaker constraints on Ωm. Constraints on S8 are, however, fairly insensitive to the choice of ns prior.

Turning to astrophysical priors, the galaxy bias parameter top-hat priors on b₁, b₂, ${\gamma }_3^{-}$, (Eq. (3)), and on a_vir (see the ‘fingers of god’ model in Eqs. (6)–(9) of Joachimi et al. 2021), match those adopted in Tröster et al. (2020b), which cover a wider range than those used in Sánchez et al. (2017). The two BOSS redshift slices have independent sets of parameters. Wide uniform priors for the intrinsic alignment parameter A_IA are chosen to be uninformative. Uniform priors on the baryon feedback parameter A_bary are chosen such that the resulting Mead et al. (2016) model of the non-linear matter power spectrum encompasses both the most aggressive feedback model from the van Daalen et al. (2011) suite of hydrodynamical simulations, along with the dark matter-only case where A_bary = 3.13.

There are five additional correlated nuisance parameters, ${\delta }_z^i$, that model uncertainty in the mean of the source redshift distributions. We adopt a multivariate Gaussian prior for the vector δ_z with a mean μ = (0.0001,  0.0021,  0.0129,  0.0110,   − 0.0060), and a covariance, C_δz, as calibrated using mock galaxy catalogues in Wright et al. (2020a). The diagonal terms of C_δz are typically at the level of ∼(0.01)², with off-diagonal correlation coefficients ranging between ∼0.1 and ∼0.3 (see Sect. 3 and Fig. 2 of Hildebrandt et al. 2021, for details).

Appendix C: Parameter constraints

In this appendix we tabulate the maximum posterior (MAP) and marginalised constraints on the flat ΛCDM cosmological parameters, in Table C.1, for the different combinations of the three large-scale structure probes considered in this work. For constraints from KiDS-1000 cosmic shear alone, we refer the reader to Asgari et al. (2021a).

Appendix D: Modelling intrinsic galaxy alignment

In this analysis we adopt the Bridle & King (2007) NLA model in order to marginalise over our uncertainty in the contribution to the observed two-point shear correlation function from the intrinsic alignment (IA) of galaxies within their surrounding density field. More sophisticated models exist, however, and in this appendix we briefly discuss these alternatives. We then provide justification for our choice by summarising the analysis of Fortuna et al. (2021) which demonstrates that for the statistical power of KiDS-1000, the use of the somewhat adhoc NLA model is sufficiently flexible so as not to introduce any biases in the cosmological parameter constraints.

There are two advanced methods to determine a non-linear IA model. The first uses a perturbative approach to model the non-linear behaviour of the tidal alignment (where the galaxy is preferentially aligned with the ‘stretching axis’ of the tidal quadrupole) and tidal torquing (where the galaxy disc forms perpendicular to the angular momentum axis which is dependent on the tidal field) (Blazek et al. 2019; Vlah et al. 2020). The second uses a halo model approach, which introduces a model for the small-scale alignment of satellite galaxies within central haloes (Schneider & Bridle 2010), and allows for different alignment strengths to be included for the evolving red and blue galaxy population (Fortuna et al. 2021). On large physical scales both techniques, along with the NLA model, recover the linear alignment model of Hirata & Seljak (2004). On small physical scales the accuracy of the perturbative approach is limited by the order of the corrections adopted (Blas et al. 2013). The accuracy of the halo model approach is limited by the challenge of modelling the transition between the properties of galaxies within single halos and across the full density field, commonly referred to as the one-to-two halo transition (see for example the discussion in Mead et al. 2021). Both models are nevertheless a significant improvement on the NLA model, which developed an adhoc solution to a small-scale mismatch of the Hirata & Seljak (2004) model with IA numerical simulations (Heymans et al. 2004), by replacing the linear matter power spectrum in the Hirata & Seljak (2004) formalism with the non-linear matter power spectrum.

The IA observations of red and blue galaxies have been both direct (Joachimi et al. 2011; Mandelbaum et al. 2011; Singh et al. 2015; Tonegawa et al. 2018; Johnston et al. 2019) and indirect (Heymans et al. 2013; Samuroff et al. 2019), recovering a common conclusion of significant alignment for red galaxies, and, as yet, no detection of alignment for blue galaxies. Direct observations are limited by the depth of the spectroscopic, or high-accuracy photometric redshift, galaxy samples that can be studied. Indirect observations, where an IA model is constrained simultaneously with a cosmological model, are limited by degeneracies with nuisance parameters. As the IA and cosmological signal scale differently with redshift, a flexible IA model can absorb any systematic errors in the shape of the source redshift distributions, such that the resulting IA parameter constraints are not a true reflection of the underlying IA model. This point is nicely illustrated in Fig. 5 of Efstathiou & Lemos (2018) where indirect constraints on the amplitude of the NLA model for the full galaxy sample are shown to vary from −6 < A_IA < 6 across a wide range of published cosmic shear surveys. Wright et al. (2020b) present another example, analysing mock galaxy catalogues to improve the accuracy of the priors on the redshift uncertainty for the previous KiDS data release. In this case the fiducial IA constraint A_IA = 0.95 ± 0.67 reduces to A_IA = 0.28 ± 0.59 with the inclusion of a more accurate redshift uncertainty model. As the intrinsic alignment model depends on the cosmology, any tendency for the model to incorrectly absorb redshift errors would inadvertently lead to unexpected biases in the cosmological parameter constraints. The more freedom afforded to the IA model, the more opportunity there is for such biases to occur. Given this concern, we choose to adopt the minimal IA model freedom afforded by current direct observational IA constraints. In this way any unaccounted errors in our redshift distributions can be detected through a poor goodness-of-fit of the combined cosmological and IA model, as seen, for our second tomographic redshift bin, in the KiDS-1000 internal consistency analysis of Asgari et al. (2021b, Appendix B.2).

Our choice of the one-parameter NLA model is motivated by the IA halo model analysis of Fortuna et al. (2021). In this analysis central red galaxies are modelled using the Hirata & Seljak (2004) model, with an amplitude A_red and a luminosity dependent scaling ∝L^β. Combining the observed constraints from Joachimi et al. (2011), Singh et al. (2015) they model a simple power-law scaling model with A_red = 5.3 ± 0.6 and β = 1.2 ± 0.4. The additional constraints from Johnston et al. (2019) motivate a broken power-law model with A_red = 5.1 ± 1.0 and β_L ≥ L0 = 1.2 ± 0.4, which they also consider. Central blue galaxies follow the Hirata & Seljak (2004) model, with a Johnston et al. (2019) constrained amplitude of A_blue = 0.2 ± 0.4 and no luminosity dependence, as there is no observational evidence to support this. The small-scale satellite galaxy alignments are modelled following Schneider & Bridle (2010). The amplitude, radial and luminosity dependence of the alignment of red and blue satellite galaxies within their host halo is independently constrained following Georgiou et al. (2019). The large-scale galaxy alignment signal is considered to be sourced solely by central galaxies, as there is currently no observational evidence to support large-scale (two-halo) alignments between the different satellite populations.

The direct IA observational constraints adopted by Fortuna et al. (2021) are determined from low redshift galaxy samples, where the properties of the full galaxy population differ significantly from the high redshift galaxy population analysed in KiDS-1000. Given that the fundamental physical processes that underpin the intrinsic alignment mechanisms are unlikely to significantly evolve out to z ∼ 1, however, these galaxy-type-specific constraints are relevant at the redshifts sampled by KiDS, provided the evolution in the relative fractions, luminosities and distribution of the red and blue, central and satellite, galaxy populations are accurately modelled.

Fortuna et al. (2021) determine a sample of plausible IA models for a KiDS-like survey. They combine the different galaxy-type IA contributions using a halo occupation distribution model taken from the MICE mock galaxy catalogues (Fosalba et al. 2015), with a magnitude limit corresponding to KiDS-depth. The range of full-population models encompasses the observational uncertainty on the IA model parameters for each galaxy type. Analysing these IA halo-models with the NLA model, they find the highest NLA model amplitude of A_IA = 0.44 ± 0.13 when adopting a broken power-law for the red-central luminosity scaling. This is consistent with the Asgari et al. (2021b) KiDS-1000 COSEBIs cosmic shear constraints with $A_{\mathrm{IA}}=0.{26}_{-0.34}^{+0.42}$, the KiDS-1000 band power cosmic shear constraints with $A_{\mathrm{IA}}=0.{97}_{-0.38}^{+0.29}$ and the KiDS-1000 3 × 2 pt band power constraints with $A_{\mathrm{IA}}=1.{07}_{-0.31}^{+0.27}$, with the largest 1.9σ difference found for the 3 × 2 pt results.

Adopting the NLA model in a cosmological parameter inference of a mock KiDS-like data vector of the cosmic shear signal contaminated by the range of different IA halo models allowed by observations, Fortuna et al. (2021) conclude that the redshift dependence of the true IA halo model is not large enough to bias the cosmological parameters in a KiDS-like cosmic shear analysis with the NLA model. Asgari et al. (2021b) nevertheless explore one extension of the NLA model, with the inclusion of a redshift-dependent scaling term, which does not change the recovered KiDS-1000 cosmic shear MAP S₈ value. It only serves to increase the marginal credible region of S₈ by ∼10%. Given these analyses, we choose to adopt a minimal one-parameter NLA model in our 3 × 2 pt analysis, but recognise that in future surveys such as LSST and Euclid, this adhoc model will no longer be sufficient given the expected statistical power of these next-generation surveys (Blazek et al. 2019; Fortuna et al. 2021).

Appendix E: Sensitivity tests

In this appendix we present, with Fig. E.1, the marginalised posterior distributions for a selection of the 3 × 2 pt sensitivity tests explored in Sect. 3. These tests were all shown to recover consistent constraints on S₈, in Fig. 5. Here we explore these tests in more detail.

Fig. E.1. Marginalised posterior distributions for the extended set of cosmological parameters shown in Fig. 6, comparing the fiducial 3 × 2 pt analysis (red) to a selection of our sensitivity test analyses where we ignore the impact of baryon feedback (the ‘No baryon’ case, sea-green), limit the analysis to a linear galaxy bias model (the ‘No higher order GC’ case, lime-green), and remove individual tomographic bins from our weak lensing observables (orange, purple and pink).

We compare our analysis which fully marginalises over our uncertainty in the baryon feedback parameter (red), with our ‘No baryon’ case (sea-green), where A_bary = 3.13, corresponding to the non-linear matter power spectrum for a dark-matter only cosmology. Here we find very little difference, as our 3 × 2 pt analysis already favours high values of A_bary. Our choice of scales and adoption of the band-power cosmic shear statistic for this analysis also makes us less sensitive to uncertainties in the baryon feedback parameter, compared to a standard two-point correlation function analysis (Asgari et al. 2020).

The removal of the two highest photometric redshift bins (blue and purple) primarily impacts S₈. These two bins carry the majority of the signal-to-noise in our analysis, and so it is not surprising that the removal of nearly half the constraining data, in each case, can result in ∼1σ changes in the recovered S₈ (see for example the discussion in Appendix F). We refer the reader to Asgari et al. (2021a) where we present a detailed internal consistency analysis of the cosmic shear signal, following Köhlinger et al. (2019), (see also Efstathiou & Lemos 2018), concluding that the two highest photometric redshift bins are consistent with the full data set. A potential ∼3σ-level flag is, however, raised in Asgari et al. (2021a) over the internal consistency of the second tomographic bin. In our analysis where we remove the two lowest photometric redshift bins (orange), we find that these bins contribute very little to the S₈ constraint and only serve to tighten the constraints on the intrinsic alignment parameter A_IA.

Finally we turn to the galaxy bias test (lime green), where we limit the analysis to a linear galaxy bias model, b₁, setting all higher-order bias terms in Eq. (3) to zero²⁰, as well as imposing a Gaussian galaxy velocity distribution by setting a_vir to zero. We can see that this biases the recovered cosmological constraints, as the amplitude of the linear galaxy bias increases in an attempt to model the enhanced power on small scales that the non-linear bias induces. The erroneous increase in b₁ leads to a decrease in the recovered S₈, and also an overall reduction in the goodness-of-fit of the model with ${\chi }_{\text{MAP}}^2$ = 379.8 for ∼300 degrees of freedom. This result should serve as a point of caution for 3 × 2 pt analyses that adopt an effective linear bias model (see also the discussion in Asgari et al. 2021b), although we note that our analysis is particularly sensitive to the galaxy bias model given the high signal-to-noise BOSS clustering observations that probe physical scales as low as s_min = 20 h⁻¹ Mpc.

Appendix F: Expected S₈ differences between partially overlapping weak lensing surveys

In this appendix we construct a simple model to estimate the expected statistical fluctuation in S₈ constraints from partially overlapping weak lensing surveys, specifically our previous (KV450, Wright et al. 2020b), and current KiDS analyses. Assuming that we only wish to constrain the overall amplitude of the measurement, the inference can be approximated by a linear least-squares problem, where S₈, the parameter of interest, is Gaussian distributed.

In the following, measurements derived from the KV450 footprint are denoted by an X in the subscript, while measurements on the new area added for KiDS-1000 are denoted with a Y. Measurements derived from the full KiDS-1000 footprint are denoted with a Z. Now, let $S_X^{\mathrm{KV}450}$ denote the value of S₈ estimated from the KV450 footprint, using the KV450 methodology; $S_Y^{\mathrm{KiDS}\text{-}1000}$ the value estimated from the newly added area using the KiDS-1000 methodology; and $S_Z^{\mathrm{KiDS}\text{-}1000}$ the KiDS-1000 value estimated from the full footprint. Approximating the estimates from disjoint footprints as independent, the S₈ measurement from the full KiDS-1000 area is then given by the area-weighted average²¹

$$\begin{array}{c}\hfill S_Z^{\mathrm{KiDS}\text{-}1000}=\frac{A_X{S}_X^{\mathrm{KiDS}\text{-}1000}+A_Y{S}_Y^{\mathrm{KiDS}\text{-}1000}}{A_Z}·\end{array}$$(F.1)

Here A_X is the effective area of the KV450 footprint, A_Z is the effective area of the KiDS-1000 footprint, and A_Y = A_Z − A_X, is the additional area added between the two data releases. The uncertainty ${\sigma }_X^{\mathrm{KiDS}\text{-}1000}$ on the measurement $S_X^{\text{KiDS01000}}$ is then related to the uncertainty ${\sigma }_Z^{\text{KiDS-1000}}$ on measurements $S_Z^{\text{KiDS-1000}}$, as

$$\begin{array}{c}\hfill {\sigma }_X^{\mathrm{KiDS}\text{-}1000}=\sqrt{\frac{A_Z}{A_X}}{\sigma }_Z^{\mathrm{KiDS}\text{-}1000},\end{array}$$(F.2)

and analogously for ${\sigma }_Y^{\mathrm{KiDS}\text{-}1000}$. It is important to note that the uncertainty ${\sigma }_X^{\mathrm{KV}450}$ differs from ${\sigma }_X^{\mathrm{KiDS}\text{-}1000}$ due to differences in the KiDS-1000 and KV450 methodologies, such as the adopted m-calibration uncertainty.

We define $\mathrm{\Delta }=S_Z^{\text{KiDS-1000}}-S_X^{\text{KV450}}$, the offset between the KiDS-1000 and KV450 S₈ measurements. The uncertainty on Δ is then given by

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{\Delta }}^2={\left({\sigma }_X^{\mathrm{KV}450}\right)}^2+{\left({\sigma }_Z^{\mathrm{KiDS}\text{-}1000}\right)}^2-2\sqrt{\frac{A_X}{A_Z}}{\sigma }_X^{\mathrm{KV}450}{\sigma }_Z^{\mathrm{KiDS}\text{-}1000},\end{array}$$(F.3)

where we have assumed a 100% correlation between the KiDS-1000 and KV450 measurements within the KV450 footprint area X, but approximate the areas X and Y as fully uncorrelated. This approximation neglects the large-scale correlations between X and Y. Given that the majority of the error budget for KiDS-1000 stems from the random shape noise component, which is independent between X and Y, this approximation is sufficient for this toy model. The absolute difference, |Δ| has the cumulative distribution function, CDF,

$$\begin{array}{c}\hfill \mathrm{CDF}(x)=\mathrm{erf}\left(\frac{x}{{\sigma }_{\mathrm{\Delta }}\sqrt{2}}\right),\end{array}$$(F.4)

and an expectation value $E[|\mathrm{\Delta }|]=\sqrt{2/\pi }{\sigma }_{\mathrm{\Delta }}$.

We now compare the actual KV450 and KiDS-1000 S₈ constraints, given the effective areas A_KV450 = A_X = 341.3 deg² and A_KiDS-1000 = A_Z = 777.4 deg². Using the marginal S₈ constraints from the two-point shear correlation function analysis, comparing Wright et al. (2020b, KV450: $S_8=0.{716}_{-0.038}^{+0.043}$) with Asgari et al. (2021a, KiDS-1000: $S_8=0.{768}_{-0.020}^{+0.016}$), we find |Δ| = 0.052 = 1.6σ_Δ, which can be compared against the expected offset of E[|Δ|] = 0.026. We expect to find an offset of this, or a larger magnitude, 10% of the time.

This simple model analysis is sufficient to conclude that the increase in S₈ that we find between KV450 and KiDS-1000 is consistent with the expectation from simple statistical fluctuations. A complete assessment could be conducted by analysing a reasonable fraction of the 20 000 KiDS mock catalogues from Joachimi et al. (2021), with the KiDS-1000 and KV450 footprints imposed. Unfortunately this is, however, out-of-scope for this analysis, given the significant compute power that would be incurred.

Appendix G: Tension estimators

In this appendix we review the range of tension estimators employed in Sect. 3.2 to quantify the tension between the cosmological parameter constraints from Planck and our KiDS-1000-BOSS 3 × 2 pt analysis. These are in addition to the standard Gaussian offset measure, τ, in Eq. (8), and the Bayes factor, R, in Eq. (9).

G.1. Hellinger tension measure: d_H

To assess the tension between the marginal S₈ distributions of two experiments, we have employed the widely used Gaussian-approximation measure, τ, given in Eq. (8). As an alternative that does not depend on the Gaussian assumption, we construct a sample-based, one-dimensional tension estimator adopting the Hellinger distance d_H. The Hellinger distance is widely used in statistics as a stable metric in optimisation problems that require the comparison of distributions (see for example Beran 1977). It is given by

$$\begin{array}{c}\hfill d_{\mathrm{H}}^2\left[p;q\right]=\frac{1}{2}{\displaystyle \int }\mathrm{d}x{\left[\sqrt{p(x)}-\sqrt{q(x)}\right]}^2,\end{array}$$(G.1)

where the integral runs over the support of the PDFs of the distributions p and q. We choose this distance definition as it is symmetric, avoids PDFs in the denominator and the associated instability for sample estimates, and has a closed-form expression for the case of two Gaussians:

$$\begin{array}{c}\hfill d_{\mathrm{H}}^2\left[\mathcal{N}({\mu }_1,{\sigma }_1^2);\mathcal{N}({\mu }_2,{\sigma }_2^2)\right]=1-\sqrt{\frac{2{\sigma }_1{\sigma }_2}{{\sigma }_1^2+{\sigma }_2^2}}exp\{-\frac{{({\mu }_1-{\mu }_2)}^2}{4({\sigma }_1^2+{\sigma }_2^2)}\},\end{array}$$(G.2)

where the normal distributions have means μ_1, 2, and variances ${\sigma }_{1,2}^2$.

We estimate the Hellinger distance between our samples via a discretised version of Eq. (G.1). The PDFs are built from the samples in two ways: through binning into a histogram (with the number of bins chosen to scale with the square root of the number of samples), and through kernel density estimation (using a Gaussian kernel with the width chosen by Scott’s rule). The former approach will in general over-estimate tension as tails are underpopulated due to the discrete sampling, whereas the latter approach under-estimates tension due to the smoothing of the PDFs by the kernel density estimator. Empirically, we find that a simple average of the two approaches yields good accuracy over the few-σ tension range that we are interested in.

We invert Eq. (G.2) to determine the expected difference in the means, |μ₁ − μ₂| if the two distributions were Gaussian, for a given value of d_H. To do so, we use the standard deviations estimated from the samples. The resulting difference in the means is inserted into Eq. (8) to obtain a measure of tension in the familiar σ units²². Comparing the 3 × 2 pt and PlanckS₈ distributions, we find a Hellinger distance d_H = 0.95, which corresponds to a tension of 3.1σ.

G.2. Distribution of parameter shifts: p_S

Following Raveri et al. (2020, see also Köhlinger et al. 2019), we consider the distribution P(Δθ) of the parameter shifts Δθ = θ_Planck − θ_{3 × 2 pt + Planck}. The significance of any shift can then be assessed by calculating how much of P(Δθ) is enclosed within the iso-probability surface at P(0):

$$\begin{array}{c}\hfill p_{\mathrm{S}}={\displaystyle {\int }_{P(\mathrm{\Delta }\mathit{\theta })>P(0)}P(\mathrm{\Delta }\mathit{\theta })\mathrm{d}\mathrm{\Delta }\mathit{\theta }.}\end{array}$$(G.3)

In other words, how much of the distribution of parameter shifts P(Δθ) is more likely than the probability of no shifts. Accurately evaluating the integral in Eq. (G.3) when P(Δθ) is not centred on zero requires very large numbers of samples to sufficiently cover the tails of the distribution. For a single parameter, S₈, we find p_S = 0.9981, with a PTE of 0.0019, or 3.1σ. This is in agreement with the other single-parameter tension estimators Eqs. (8) and (G.1). For the five shared parameters between our 3 × 2 pt analysis and Planck we find lower significances, around 2σ, albeit with very large uncertainties owing to the insufficient sampling of the tails of the parameter shift distribution.

G.3. Suspiciousness: S

Handley & Lemos (2019) propose the ‘suspiciousness’ statistic S based on the Bayes factor, R, but hardened against prior dependences. They define

$$\begin{array}{c}\hfill lnS=lnR-lnI,\end{array}$$(G.4)

where lnR is the logarithm of the Bayes factor in Eq. (9):

$$\begin{array}{c}\hfill lnR=ln{Z}_{3\times 2\mathrm{pt}+\mathit{Planck}}-ln{Z}_{3\times 2\mathrm{pt}}-ln{Z}_{\mathit{Planck}},\end{array}$$(G.5)

with the evidence Z_i given by

$$\begin{array}{c}\hfill Z={\displaystyle \int }\mathcal{L}\pi \mathrm{d}\mathit{\theta },\end{array}$$(G.6)

the integral of the likelihood ℒ and prior π over the parameters θ. The information ratio I in Eq. (G.4) is defined as

$$\begin{array}{c}\hfill lnI={\mathcal{D}}_{3\times 2\mathrm{pt}}+{\mathcal{D}}_{\mathit{Planck}}-{\mathcal{D}}_{3\times 2\mathrm{pt}+\mathit{Planck}},\end{array}$$(G.7)

with 𝒟_i being the Kullback-Leibler divergence between the posterior P and prior for probe i:

$$\begin{array}{cc}\hfill \mathcal{D}& ={\displaystyle \int Pln\frac{P}{\pi }\mathrm{d}\mathit{\theta }=\int Pln\mathcal{L}\mathrm{d}\mathit{\theta }-lnZ}\hfill \\ \hfill & ={⟨ln\mathcal{L}⟩}_{\mathrm{P}}-lnZ.\hfill \end{array}$$(G.8)

The second equality follows from Bayes theorem: P = ℒπ/Z. Using this definition of 𝒟 allows us to rephrase the suspiciousness solely in terms of the expectation values of the log-likelihoods:

$$\begin{array}{cc}\hfill lnS=& {⟨ln{\mathcal{L}}_{3\times 2\mathrm{pt}+\mathit{Planck}}⟩}_{{\mathrm{P}}_{3\times 2\mathrm{pt}+\mathit{Planck}}}-{⟨ln{\mathcal{L}}_{3\times 2\mathrm{pt}}⟩}_{{\mathrm{P}}_{3\times 2\mathrm{pt}}}\hfill \\ \hfill & -{⟨ln{\mathcal{L}}_{\mathit{Planck}}⟩}_{{\mathrm{P}}_{\mathit{Planck}}}.\hfill \end{array}$$(G.9)

We find lnS = −2.0 ± 0.1. For Gaussian posteriors, the quantity d − 2lnS is distributed as ${\chi }_{\text{d}}^2$, where d = d_{3 × 2 pt} + d_Planck − d_{3 × 2 pt + Planck} is the difference in the Bayesian model dimensionalities. This allows us to assign a probability of the observed suspiciousness under the assumption that the two data sets are in concordance. We calculate d, following Handley & Lemos (2019), from the variances of the log-likelihoods. Other estimates of the Bayesian model dimensionality, such as those introduced in Raveri & Hu (2019) yield very similar results. We calculate d = 3.3 ± 0.7 and therefore conclude that the probability of observing our measured suspiciousness statistic is 0.08 ± 0.02, or 1.8 ± 0.1σ.

G.4. Goodness-of-fit change: Q_DMAP

Raveri & Hu (2019) introduce the Q_DMAP statistic that quantifies how the goodness-of-fit changes when two data sets are combined, based on the difference of the log-likelihoods, $ln\mathcal{L}({\mathit{\theta }}_i^{\mathrm{MAP}})$, at the MAP of the respective posteriors²³,

$$\begin{array}{c}\hfill Q_{\mathrm{DMAP}}=2ln\mathcal{L}\left({\mathit{\theta }}_{3\times 2\mathrm{pt}}^{\mathrm{MAP}}\right)+2ln\mathcal{L}\left({\mathit{\theta }}_{\mathit{Planck}}^{\mathrm{MAP}}\right)-2ln\mathcal{L}\left({\mathit{\theta }}_{3\times 2\mathrm{pt}+\mathit{Planck}}^{\mathrm{MAP}}\right).\end{array}$$(G.10)

For Gaussian posteriors, the sampling distribution of Q_DMAP is χ²-distributed with d_{3 × 2 pt} + d_Planck − d_{3 × 2 pt + Planck} degrees of freedom. Here we follow Raveri & Hu (2019) and calculate the model dimensionalities by $d_i=N-\mathrm{tr}[{\mathcal{C}}_{\pi }^{-1}{\mathcal{C}}_{\mathrm{P}}]$, where N is the total number of varied parameters, 𝒞_π is the prior parameter covariance, and 𝒞_P is the posterior parameter covariance. We find consistent results when using the variance of the log-likelihoods, the approach taken for the suspiciousness calculation. We find Q_DMAP = 9.5, with d = 3.6, corresponding to a probability of 0.037, or 2.1σ.

G.5. Parameter difference in multiple dimensions: Q_UDM

Raveri & Hu (2019) introduce the Q_UDM statistic that quantifies the difference between cosmological parameters across multiple dimensions. Here

$$\begin{array}{c}\hfill Q_{\mathrm{UDM}}=\mathrm{\Delta }{\overline{\mathit{\theta }}}^{\mathrm{T}}{({\mathcal{C}}_{\mathit{Planck}}-{\mathcal{C}}_{3\times 2\mathrm{pt}+\mathit{Planck}})}^{-1}\mathrm{\Delta }\overline{\mathit{\theta }},\end{array}$$(G.11)

where $\mathrm{\Delta }\overline{\mathit{\theta }}={\overline{\mathit{\theta }}}_{\mathit{Planck}}-{\overline{\mathit{\theta }}}_{3\times 2\mathrm{pt}+\mathit{Planck}}$ is the shift in the mean of the shared parameters, and 𝒞_Planck and 𝒞_{3 × 2 pt + Planck} are the posterior parameter covariances for Planck and 3 × 2 pt + Planck, respectively. The estimator Eq. (G.11) is asymmetric in the first data set if the posteriors are non-Gaussian. In this case, the data set with the more Gaussian posterior should be used, which in our case is Planck. We use tensiometer²⁴ to compute Q_UDM and its associated number of degrees of freedom. We find Q_UDM = 7.7 with d = 3. The sampling distribution is approximated by ${\chi }_{\mathrm{d}}^2$, such that the probability of the observed value of Q_UDM is 0.054, or 1.9σ.

Appendix H: Redundancy, pipeline validation and software review

In this appendix we briefly review the redundancy in the KiDS-1000 analysis, from pixels through to parameters. Multi-band pixel processing and photometry in the optical has been carried out on the full KiDS-1000 data set using two independent pipelines, ASTROWISE and THELI (Begeman et al. 2013; Erben et al. 2013). This approach allowed us to resolve a range of different issues, primarily in the astrometric solution, for the handful of problem fields that resisted automated processing. We adopt the ASTROWISE reduction for multi-band photometry and the THELIr-band reduction for object detection and weak lensing shape measurement (Kuijken et al. 2019). We consider only one shape measurement technique, lensfit (Miller et al. 2013), but this is calibrated using a series of different image simulations that vary the input galaxy properties to assess the sensitivity of the shear estimator to: the presence of blending with unresolved and undetected objects, blending due to enhanced galaxy clustering, photometric redshift selection bias, size-ellipticity correlations in the galaxy properties, varying stellar density, and the choice of smooth or realistic galaxy profiles (see Kannawadi et al. 2019; Giblin et al. 2021, for details).

Our fiducial photometric redshift calibration has been compared to two additional independent calibration approaches, and has been validated on mock photometry catalogues (see Wright et al. 2020a; van den Busch et al. 2020; Hildebrandt et al. 2021, for details). Using a series of null-tests we have validated the resulting shear-redshift catalogues in Giblin et al. (2021). Our catalogue-to-observables pipeline, based on the two-point correlation function code TREECORR (Jarvis 2015), has been validated through an independent analysis using the alternative ATHENA package (Kilbinger et al. 2014), and through mock catalogue analysis (Joachimi et al. 2021). Our cosmological inference code KCAP with CAMB (Lewis et al. 2000), has been validated against the Core Cosmology Library (CCL, Chisari et al. 2019) and through the recovery of input parameters in mock data analysis (Joachimi et al. 2021). In Asgari et al. (2021a) we also verify that we produce the same results in our cosmic-shear only analysis using a completely independent inference code based on MONTEPYTHON with CLASS (Lesgourgues 2011; Brinckmann & Lesgourgues 2019; Köhlinger et al. 2019; Hildebrandt et al. 2020). Our analysis of BOSS follows Sánchez et al. (2017), which is in good agreement with the numerous independent parallel analyses of the same DR12 data set presented in Alam et al. (2017).

In regards to software review, for the majority of cases we have adopted the ‘four-eye’ approach, collaboratively building code through a git repository, with major updates reviewed through pull requests. This approach has applied throughout the full pixels-to-parameters process, and to all types of software, both our significant tools, analysis code and, for the most part, simple paper support scripts. Our repositories are open source²⁵ for others to use and build upon, but with a caveat for users to recognise that we are not software engineers. Adopting this style of software review for KiDS-1000 has been an extremely beneficial exercise for the team, with many lessons learnt for how to improve our software engineering skills and our approach to open source collaborative coding for future projects.

Appendix I: Post-unblinding analyses

This analysis was carried out ‘blind’, such that our final key result, our constraint on S₈, was unknown until all analysis choices were fixed. Our blinding strategy creates three versions of the catalogue, where one is the truth, and the other two are modified to introduce up to a ∼ ± 2σ deviation in the recovered value of S₈ (Kuijken et al. 2015; Giblin et al. 2021). In contrast to previous KiDS analyses, we set ourselves a challenge to only run our blind data analysis inference once, after fully developing and road-testing the pipeline using mock catalogues and data vectors (Joachimi et al. 2021). Unexpectedly, the only analysis where we failed to meet this challenge were for data vectors that included galaxy clustering. Here a bug in a naming convention in an updated COSMOSIS-CAMB interface resulted in ∼0.5σ errors in the recovery of σ₈ in our initial BOSS re-analysis. This was corrected and updated for our fiducial analysis before unblinding.

Our fiducial results, presented in Fig. 4 and tabulated in Appendix C, were carried out for all three blinds using a covariance matrix derived assuming a fiducial cosmology given by the best-fit parameters from Tröster et al. (2020b). We reserved our iterated-covariance analysis, however, for post-unblinding, re-analysing the true data vector with an updated covariance matrix derived adopting the best-fit parameters from our initial blinded run of the true catalogue. This iterative step changed our value for the full 3 × 2 pt analysis S₈ by 0.2σ for the MAP constraint, and by 0.1σ for the marginal constraint. Our tension-consistency analysis with the CMB constraints from Planck was also conducted after unblinding, as we were unblind to Planck throughout the process.

As we are only interested in relative shifts in the sensitivity tests presented in Fig. E.1, these were only conducted for a single blind, and re-calculated post-unblinding for the true catalogue. Wishing to be fully transparent in this appendix, we note that owing to limited resources the sensitivity test was not updated during the blinded phase of the project in order to correct for the COSMOSIS-CAMB interface error which only impacted the BOSS constraints. We argue that this was appropriate given the nature of the test and given that this error has been corrected post-unblinding. Our cosmic shear and galaxy clustering analysis, featured in our internal consistency test in Fig. 5, was also not re-processed blind after the BOSS error was corrected, as both our forecast in Joachimi et al. (2021), and our initial analysis, confirmed that the constraints from our cosmic shear and galaxy clustering analysis were almost identical to constraints from the full 3 × 2 pt analysis, which was correctly analysed for all three blinds.

As reported in Giblin et al. (2021), co-author Kannawadi was unblinded early in the process to permit accurate calibration of the shear measurements in his re-analysis of KiDS-1000-like image simulations. Co-authors Wright and Heymans also wish to record that they independently had a suspicion over which blind was the truth based on the changes that the SOM gold selection made to the blinded KiDS-1000 catalogue effective number density and ellipticity dispersion, compared to the impact the SOM selection made to the KiDS-VIKING-450 catalogue analysed in Wright et al. (2020b). These suspicions were never shared with the rest of the team, nor discussed with each other, and on unblinding were found to both be false, and different! This issue nevertheless highlighted to us the challenges of catalogue-level blinding for successive data releases where little has changed in the core data reduction process. This is particularly relevant if the people working at the coal-face of catalogue production are the same people working on the cosmological analysis. Future KiDS blinding will therefore likely adopt the approach advocated by Sellentin (2020), where the covariance matrix is modified in the analysis.

All Tables

All Figures

Fig. 1. KiDS-1000 cosmic shear power spectra: tomographic band powers comparing the E-modes (upper left block) with the best-fit cosmological model from our combined multi-probe analysis. The tomographic bin combination is indicated in the upper right corner of each sub-panel. The null-test B-modes (lower right block – note the reduced ordinate scale), are consistent with zero for both the full data vector and each bin combination individually. The errors are estimated analytically (Joachimi et al. 2021). See Sect. 3 for a discussion on the goodness-of-fit.

In the text

	Fig. 2. BOSS-DR12 anisotropic clustering from Sánchez et al. (2017): the transverse (pink), intermediate (blue) and parallel (black) clustering wedges in two redshift bins, compared with the best-fit cosmological model from our combined multi-probe analysis. The errors, estimated from mock BOSS catalogues (Kitaura et al. 2016), are highly correlated, particularly at large scales (see Fig. A.1).
In the text

Fig. 3. KiDS-1000 galaxy-galaxy lensing power spectra: tomographic band powers comparing the E-modes (left block) with the best-fit cosmological model from our combined multi-probe analysis. The tomographic bin combination of BOSS and 2dFLenS lenses (L) with KiDS-1000 sources (S), is indicated in the upper right corner of each sub-panel. Data within grey regions are not included in the cosmological analysis. The null-test B-modes (right block – note the reduced ordinate scale), are consistent with zero for both the full data vector and each bin combination individually. The errors are estimated analytically (Joachimi et al. 2021).

In the text

Fig. 4. Marginal multi-probe constraints on the flat ΛCDM cosmological model, for the matter fluctuation amplitude parameter, σ8, the matter density parameter, Ωm, and the dimensionless Hubble parameter, h. The BOSS galaxy clustering constraints (blue), can be compared to the KiDS-1000 cosmic shear constraints (pink), the combined 3 × 2 pt analysis (red), and CMB constraints from Planck Collaboration VI (2020, grey).

In the text

Fig. 5. Constraints on the structure growth parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ for different probe combinations: 3 × 2 pt, KiDS-1000 cosmic shear, BOSS galaxy clustering, cosmic shear with galaxy-galaxy lensing (GGL), and cosmic shear with galaxy clustering. Our fiducial and preferred MAP with PJ-HPD credible interval (solid) can be compared to the standard, but shifted, marginal posterior mode with M-HPD credible intervals (dotted). Our results can also be compared to weak lensing measurements from the literature, which typically quote the mean of the marginal posterior mode with tail credible intervals (dashed).

In the text

Fig. 6. Marginalised posterior distributions for an extended set of cosmological parameters covering the matter density parameter, Ωm, the matter fluctuation amplitude parameter, σ8, the structure growth parameter, S8, the spectral index, ns, the dimensionless Hubble parameter, h, the baryon feedback amplitude parameter, Abary, the intrinsic alignment amplitude, AIA, and the linear bias parameters for the low and high BOSS redshift bins, b1. The KiDS-1000 cosmic shear results (pink), can be compared to the BOSS galaxy clustering results (blue), the combination of cosmic shear with BOSS and 2dFLenS galaxy-galaxy lensing (GGL, purple), and the full 3 × 2 pt analysis (red). The combination of cosmic shear with galaxy clustering (orange) is only distinguishable from the 3 × 2 pt result in the Abary and AIA panels. For parameters constrained by the CMB, we also include constraints from Planck Collaboration VI (2020, grey).

In the text

	Fig. 7. 3 × 2 pt constraints on S8 for a series of sensitivity tests; when we ignore the impact of baryon feedback (the ‘No baryon’ case), limit the analysis to a linear galaxy bias model (the ‘No higher order GC’ case), and remove individual tomographic bins from our weak lensing observables.
In the text

	Fig. 8. Marginalised posterior distribution in the σ8 − Ωm plane, comparing the 3 × 2 pt analyses from KiDS-1000 with BOSS and 2dFLenS, with the 3 × 2 pt analysis from DES Y1 (Abbott et al. 2018), and the CMB constraints from Planck Collaboration VI (2020). The KiDS-1000 3 × 2 pt result can also be compared to our previous KV450-BOSS analysis from Tröster et al. (2020b).
In the text

	Fig. A.1. Correlation coefficients of the 3 × 2 pt covariance matrix for the cosmic shear, 𝒞E, galaxy-galaxy lensing, 𝒞nϵ and galaxy clustering observables, ξgg (see Sect. 2 for details). Here the band powers 𝒞, are related to the angular power spectrum C, in Eqs. (1) and (3), as 𝒞 = ℓ2C/(2π).
In the text

Fig. B.1. Impact of the ns prior: Comparing marginalised posterior distributions for the BOSS galaxy clustering analysis for our fiducial analysis (blue) with the constraints when adopting an uninformative prior on ns (orange). Opening the parameter space to arguably unphysical values of ns favours higher values and weaker constraints on Ωm. Constraints on S8 are, however, fairly insensitive to the choice of ns prior.

In the text

Fig. E.1. Marginalised posterior distributions for the extended set of cosmological parameters shown in Fig. 6, comparing the fiducial 3 × 2 pt analysis (red) to a selection of our sensitivity test analyses where we ignore the impact of baryon feedback (the ‘No baryon’ case, sea-green), limit the analysis to a linear galaxy bias model (the ‘No higher order GC’ case, lime-green), and remove individual tomographic bins from our weak lensing observables (orange, purple and pink).

In the text