© The Authors 2023

1. Introduction

Gravitational lensing, the theory that describes the deflection of light by massive objects, reveals a wealth of information about the evolution of matter structure in the Universe (see, e.g. Hamana et al. 2020; Asgari et al. 2021; Amon et al. 2022, for recent cosmic shear analyses). Due to the accurate theoretical description and control over systematic inaccuracies, the most commonly used methods focus on two-point statistics, namely the two-point correlation functions and their Fourier counterparts called power spectra. These statistics are excellent for capturing the Gaussian information contained in the data and are complete if the data are Gaussian-distributed, such as the cosmic microwave background (CMB; e.g. Planck Collaboration V 2020). In the late Universe, however, non-linear gravitational instabilities generate a significant amount of non-Gaussian features, whose information can only be accessed with higher-order statistics. Furthermore, since higher-order statistics scale differently with cosmology and are affected differently by residual systematic effects, the constraining power on cosmological parameters increases by jointly investigating second- and higher-order statistics (see, e.g. Kilbinger & Schneider 2005; Bergé et al. 2010; Pires et al. 2012; Fu et al. 2014; Pyne & Joachimi 2021).

As the current analysis of the estimation of cosmological parameters reaches the per cent level, tensions arise between observations of the early and late or local Universe. A famous tension is the one for the Hubble parameter H₀ (Di Valentino et al. 2021a) but it is not the subject of this work. More interesting for us is the tension in the matter clustering parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, where it seems that the local Universe is less clustered than observations of the CMB suggest (Hildebrandt et al. 2017; Joudaki et al. 2020; Heymans et al. 2021; Di Valentino et al. 2021b).

A recent development in analysis methods has enabled the joint investigation of weak lensing and galaxy clustering data (van Uitert et al. 2018; Joudaki et al. 2018; Abbott et al. 2018; DES Collaboration 2022), yielding significantly better constraints, especially along the σ₈ − Ω_m degeneracy axis. Even though foreground clustering data introduces largely uncertain astrophysical parameters, such as the galaxy bias, that complicate the analysis, these joint analyses inform us better about the correlation between galaxies and the underlying matter distribution (Sánchez et al. 2017). Here again, two-point statistics have been favoured so far for the reasons mentioned above, such that the combination of all measurements (cosmic shear, galaxy clustering, and galaxy-galaxy lensing) are generally referred to as ‘3×2pts’ statistics.

To access the additional information contained in the non-Gaussian features, a competitive statistic to the 3×2pts method was recently proposed in Gruen et al. (2016), coined the ‘density-split statistics’ (DSS hereafter). This technique measures the tangential shear on the full pixelated survey footprint and bins the resulting shear profiles as a function of the foreground mass density. For example, high galaxy density regions generally trace large matter over-density regions, in which the tangential shear is expected to be larger, and this varies with cosmology. The DSS, therefore, captures information both from the shape and amplitude of the shear profiles and from the number of foreground galaxies in each density bin, with the latter helping to measure the galaxy bias significantly.

The first ingredient needed is a prediction model to interpret the measurements and constrain cosmological and astrophysical parameters. This can be constructed either from simulations (see, e.g., Harnois-Déraps et al. 2021; Zürcher et al. 2022, for examples of simulation-based inference using the lensing peak count) or from analytical calculations, where for instance Reimberg & Bernardeau (2018) and Barthelemy et al. (2021) made use of large deviation theory (LDT) to model the reduced-shear correction to the aperture mass probability distribution function (PDF). Another approach to access higher-order moments is discussed in Halder et al. (2021), Halder & Barreira (2022), and Heydenreich et al. (2022), where third-order cosmic shear statistics were modelled directly from the bispectrum. Although, the simulation-based approach has advantages regarding the numerical incorporation of critical systematic effects such as the intrinsic alignment (IA) of galaxies (see, e.g. Harnois-Déraps et al. 2022, hereafter HD22) and baryonic processes extracted from hydrodynamical simulations. However, it typically requires large simulation suites that jointly vary all the parameters under consideration. On the other hand, analytical modelling of the DSS can better dissect the basic underlying properties of the LSS, and it can be computed sufficiently fast enough at any point in the cosmological parameter space. Such a model was derived in Friedrich et al. (2018, hereafter F18), based on non-perturbative modelling of the matter density PDF. For a given cosmology, mean foreground galaxy density, and redshift distributions of the foreground and background galaxies, the F18 model predicts the mean tangential shear profiles and the PDF of the galaxy counts in each mass density bin. In Gruen et al. (2018, hereafter G18), the F18 model was used to constrain cosmological parameters from measurements of the Dark Energy Survey (DES) First Year and Sloan Digital Sky Survey (SDSS) data, yielding results competitive with the main DES 3×2 pt analysis (Abbott et al. 2018).

To date, no cosmological constraints from DSS exist, except that of G18. However, the methods have been improved significantly. In particular, Brouwer et al. (2018) have presented a contemporary measurement of the DSS extracted from the third data release of the Kilo-Degree Survey data (KiDS), wherein the foreground galaxies were selected to mimic the spectroscopic Galaxy And Mass Assembly survey (hereafter GAMA; Driver et al. 2011). They developed an optimal methodology in their work, notably showing how the resulting signal-to-noise ratio (S/N) depends on the smoothing scale for the density map of foreground galaxies.

Burger et al. (2022, hereafter B22) modified the analytical model by F18 for an application to galaxy density fields smoothed with general filters. As discussed in Burger et al. (2020), compensated filter functions outperform the previously used top-hat filter functions in terms of the overall S/N of the shear signals and in recovering the correlation between the galaxy and matter density contrast. B22 mention another advantage of compensated filter functions: they are more compact in Fourier space and, therefore, can better suppress large-ℓ modes where baryonic effects play an important role, as studied in Asgari et al. (2020). On the downside, compensated filters complicate the LDT-like calculations (Barthelemy et al. 2021). Nevertheless, B22 show that the density split statistics with compensated filters can still be accurately modelled in a computationally tractable manner after calibrating residual inaccuracies at large and small scales on the simulations of Takahashi et al. (2017).

The current paper presents the first cosmological inference based on a DSS analysis of the KiDS data. We exploited the model advances presented in B22, using the dense sample of bright galaxies presented in Bilicki et al. (2021), to construct our foreground density maps and compute the tangential shear from the lensing catalogue constructed from the fourth KiDS data release. Our inference includes a marginalisation over several residual systematic uncertainties. We verified with numerical N-body and hydrodynamical simulations that our measurements are robust against the IA of galaxies and baryonic feedback.

This work is structured as follows. In Sect. 2 we review the basics of the DSS and introduce small modifications made to our model compared to the one from B22. In Sect. 3 we present the observed data used in our analysis, and then we describe in Sect. 4 the simulations needed for the validation of our inference pipeline which is described in Sect. 5. In Sect. 6 we perform our validation of the model together with an investigation on IA and baryonic physics which could potentially contaminate our results. In Sect. 7 we finally present our main results and conclude with a discussion and summary in Sect. 8.

2. Theoretical background

The DSS essentially measures the tangential shear around sub-areas of the sky that are assigned according to the galaxy foreground density. It is therefore closely related to aperture statistics, which we introduce here first. Given a convergence field κ(θ), the aperture mass map is defined as

$$\begin{array}{c}\hfill M_{\mathrm{ap}}\left(\mathit{\theta }\right){\displaystyle \int {\mathrm{d}}^2}{\theta }^{\prime }\kappa (\mathit{\theta }+{\mathit{\theta }}^{\prime })U(|\mathit{\theta }\mathbf{\prime }|),\end{array}$$(1)

where θ is the position on the flat sky, and U(ϑ) is a compensated, axisymmetric filter function, such that ∫ϑ U(ϑ) dϑ = 0. The aperture mass, M_ap, can also be expressed in terms of the tangential shear γ_t (Schneider 1996) and a second filter function Q as

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\theta })={\displaystyle \int {\mathrm{d}}^2}{\theta }^{\prime }{\gamma }_{\mathrm{t}}(\mathit{\theta }+{\mathit{\theta }}^{\prime })Q(|{\mathit{\theta }}^{\prime }|),\end{array}$$(2)

where

$$\begin{array}{c}\hfill Q(\vartheta )=\frac{2}{{\vartheta }^2}{\displaystyle \underset{0}{\overset{\vartheta }{\int }}}\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }U({\vartheta }^{\prime })-U(\vartheta ).\end{array}$$(3)

The above relation between the two filters U and Q can be inverted,

$$\begin{array}{c}\hfill U(\vartheta )=2{\displaystyle \underset{\vartheta }{\overset{\infty }{\int }}}\mathrm{d}{\vartheta }^{\prime }\frac{Q({\vartheta }^{\prime })}{{\vartheta }^{\prime }}-Q(\vartheta ),\end{array}$$(4)

allowing us to work either with convergence maps or shear catalogues. Replacing the convergence by the foreground galaxy number counts n(θ) in Eq. (1), we define the aperture number counts, or simply aperture number, as (Schneider 1998)

$$\begin{array}{c}\hfill N_{\mathrm{ap}}(\mathit{\theta }){\displaystyle \int {\mathrm{d}}^2}{\theta }^{\prime }n(\mathit{\theta }+{\mathit{\theta }}^{\prime })U(|{\mathit{\theta }}^{\prime }|).\end{array}$$(5)

This definition is equivalent to the ‘Counts-in-Cell’ (CiC) statistics mentioned in Gruen et al. (2016) if the filter U is defined as a top-hat. In that case, however, U is not compensated; hence, one can not relate the filters U and Q.

The general idea of the DSS is to divide the survey area into quantiles 𝒬 according to the aperture number N_ap and then measure the mean tangential shear in the corresponding quantiles ⟨γ_t|𝒬⟩. We detail in Sect. 2.1 how we achieve this in the data and in Sect. 2.2 how we predict it analytically.

2.1. Measuring the DSS vector

We followed several ordered pipeline steps to extract the DSS data vector for the cosmological inference. First, we distributed the foreground (lens) galaxies onto a HEALPix (Górski et al. 2005) grid n(θ) of nside = 4096, which resulted in a pixel area of A_HP ≈ 0.74 arcmin². Second, we determined the aperture number field N_ap with a filter function U, with a finite filter radius Θ and maximal one transition from positive to negative values at θ_tr. This was achieved with the healpy function smoothing, with a beam window function that is the U-filter in the spherical harmonic space determined with healpy function beam2bl. Since Eq. (5) assumes full knowledge of n(θ) on the sky which has to be modified in the presence of a mask m(θ) as

$$\begin{array}{cc}\hfill N_{\mathrm{ap}}(\mathit{\theta })& =\frac{{\displaystyle {\int }_0^{{\theta }_{\mathrm{tr}}}}U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime }}{{\displaystyle {\int }_0^{{\theta }_{\mathrm{tr}}}}m(\mathit{\theta }+{\mathit{\theta }}^{\prime })U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime }}{\displaystyle {\int }_0^{{\theta }_{\mathrm{tr}}}}n(\mathit{\theta }+{\mathit{\theta }}^{\prime })U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime }\hfill \\ \hfill & +\frac{{\int }_{{\theta }_{\mathrm{tr}}}^{\mathrm{\Theta }}U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime }}{{\int }_{{\theta }_{\mathrm{tr}}}^{\mathrm{\Theta }}m(\mathit{\theta }+{\mathit{\theta }}^{\prime })U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime }}{\displaystyle {\int }_{{\theta }_{\mathrm{tr}}}^{\mathrm{\Theta }}}n(\mathit{\theta }+{\mathit{\theta }}^{\prime })U({\mathit{\theta }}^{\prime }){\mathrm{d}}^2{\theta }^{\prime },\hfill \end{array}$$(6)

where in this work m(θ) was the KiDS-1000 mask. For example, the second part of this equation vanishes for a top-hat filter. We divided the filter in this way to prevent the masked area from entering the positive part to decrease N_ap systematically and artificially increase the aperture number when entering the negative filter region. By separating the compensated filter into its positive and negative parts, we could correct both individually. Furthermore, since this correction became less accurate in heavily masked regions, we included only those pixels where the number of unmasked pixels within the given filter radius (effective area) was greater than 50% of the total number of pixels inside the same circle (maximal area), which we treated as our ‘good’ pixels. This can change from pixel to pixel because our HEALPix map originated from a flat sky mask. To avoid a pixel being considered good, yet for more than 50% of pixels to be missing in the negative part, we included for compensated filters only those pixels where the effective area for the positive and the negative part was greater than 50% of their individual maximal area. With our choice of the 50% threshold, we attempted to achieve a compromise between statistical power and falsely measured N_ap values. G18 considered only regions with at least 80% coverage, but since the KiDS footprint is very narrow, we had to relax that threshold to avoid shot noise-dominated data vectors. Therefore, we decided to use the highest threshold that yielded shear profiles that do not deviate significantly from shear profiles with smaller threshold values. The result is shown in Fig. A.1, where it is seen that the shear profiles with threshold values of 50% or smaller are quite similar and start to deviate for higher threshold values. Next, we allocated those good pixels to five quantiles 𝒬 according to their N_ap value. The pixels from each quantile are then correlated with the tangential shear information from the source catalogues using the treecorr (Jarvis et al. 2004) software in ten log-spaced bins with angular separation 10 arcmin <  ϑ <  120 arcmin. This resulted in measurements of the five tangential shear profiles ⟨γ_t|𝒬⟩, that is, one per quantile. For all measured profiles, the shear around random points is subtracted, which ensures that the average overall quantiles vanish by definition. Finally, we constructed our data vector, which consists of the shear profiles from the highest two and lowest two quantiles, plus the mean of the aperture number values in the same four quantiles. We must exclude the information of one quantile 𝒬 since the other four quantiles fully determine it by construction, for the reason explained above. The same is true for mean aperture number values in those quantiles, whose average is fixed by the total galaxy number density measured in the data. This results in ten θ bins, four quantiles and two source bins in a data vector of 80 + 4 = 84 elements. Although the information content is independent of which quantile is removed, we exclude the middle quantile for the whole analysis since it would have the least cosmological information, if the quantiles were analysed separately.

2.2. Modelling the DSS vector

Our modelling of the DSS signal is inspired by the LDT approach and builds from the original F18 model and the subsequent improvements presented in B22. We refer the reader to these two references for the complete details of the model calculations and highlight here only the broad principles and the minor modifications we have made. Briefly, the model consists of three key ingredients: (i) the PDF of the matter density contrast, smoothed with the filter function U, labelled δ_m, U; (ii) the expectation value of the convergence inside a radius ϑ given the smoothed matter density contrast defined above; (iii) the distribution of N_ap values given δ_m, U. B22 shows how these are computed for arbitrary filter functions and quantile counts. We, however, focus here on the ‘adapted compensated’ filter case, introduced in B20 and shown in Fig. 3 in B22, and five quantiles. As in B22, the model was calibrated using the full-sky simulations described in Takahashi et al. (2017) to suppress the residual differences of the modelled and measured data vector. The KiDS-1000 lens distribution used in the current paper peaks at a lower redshift than that in B22, and we know that the DSS model is slightly less accurate in that case, but we subsequently show that the calibration is accurate enough to yield unbiased results.

Since real galaxies are not expected to be perfectly Poisson-distributed, we modified the distribution of the aperture number computed in the B22 model in a way that allows for super-Poissonian shot noise. Inspired by F18, we achieved this by scaling the galaxy number density n₀ with a free parameter α >  0, such that n₀ α⁻¹ can be interpreted as an effective number density of Poissonian tracers. This implies that the quantity p(N_ap α⁻¹|δ_m, U) follows a log-normal distribution, instead of p(N_ap|δ_m, U) as in B22. Consequently, the characteristic function Ψ, determining the parameters of the log-normal distribution, must be modified to (see Eq. (36) of B22)

$$\begin{array}{cc}\hfill \mathrm{\Psi }(t)=& exp\left(2\pi \frac{n_0}{\alpha }{\displaystyle {\int }_0^{\infty }}\mathrm{d}\vartheta \vartheta (1+b⟨w_{\vartheta }|{\delta }_{\mathrm{m},U}⟩)[{\mathrm{e}}^{\mathrm{i}tU(\vartheta )}-1]\right),\hfill \end{array}$$(7)

where w_ϑ is the mean 2D density contrast on a circle at ϑ (see Eq. (37) in B22), and b is the linear galaxy bias. Moreover, to ensure that the mean aperture number remains constant, we further modified the calculation of p(N_ap|δ_m, U) as (see Eq. (A39) of B22):

$$\begin{array}{c}\hfill p\left(N_{\mathrm{ap}}|{\delta }_{\mathrm{m},U}\right)\to \frac{1}{\alpha }p\left(N_{\mathrm{ap}}{\alpha }^{-1}|{\delta }_{\mathrm{m},U}\right).\end{array}$$(8)

Since the expectation value ⟨N_ap|δ_m, U⟩∝α and the variance ⟨(N_ap − ⟨N_ap|δ_m, U⟩)²|δ_m, U⟩∝α² we see that the ratio of variance to expectation value is proportional to α as required to describe deviations from Poissonian samples. Similar to Friedrich et al. (2018), we require α >  0.1 in our parameter sampling for numerical reasons. We also compared our definition of this α parameter to the one implemented in Friedrich et al. (2018) and found no differences in the predictions.

Compared with numerical simulations, this model has been shown in B22 to be accurate, with residual inaccuracies to be everywhere significantly smaller than the statistical noise of the KiDS-1000 data. We, therefore, do not need to include a modelling error in our uncertainty budget. Furthermore, we also tested a non-linear galaxy bias model, where we exchanged the constant galaxy bias b with b = b₁ + b₂δ_m, U >  0. However, b₂ was highly correlated with other parameters such as Ω_m which prevented our parameter estimation from converging and therefore had to be excluded for this analysis. We also tested if the assumption of a linear galaxy bias is satisfied (see Fig. A.2 and its description), and the results of that test can be summarised as follows: a linear galaxy bias model is sufficient if an analysis using shear and N_ap information gives similar cosmological results as using only shear information since the shear profiles are basically insensitive to the galaxy bias model.

3. Observational data

In our analysis, we exploit the fourth data release of the KiDS (Kuijken et al. 2015, 2019; de Jong et al. 2015, 2017), which is a public survey carried out at the European Southern Observatory¹. KiDS was designed for weak lensing applications, producing high-quality images with VST-OmegaCAM camera. Thanks to the infrared data from its overlapping partner survey VIKING (VISTA Kilo-degree Infrared Galaxy survey, Edge et al. 2013), galaxies are observed in nine optical and near-infrared bands, u, g, r, i, Z, Y, J, H, K_s, allowing for better control over redshift uncertainties (Hildebrandt et al. 2021, hereafter H21) to earlier releases. The weak lensing data in KiDS DR4 are collectively called ‘KiDS-1000’ as they cover ∼1000 deg² of images; this reduces to 777.4 deg² of the effective area after masking. These galaxies are further split into lens and source samples, which we discuss in more detail in the following sections, with properties summarised in Table 1.

3.1. Lens catalogues

Our primary lens catalogue is the ‘KiDS-bright’ sample described in Bilicki et al. (2021, hereafter Bi21), a flux-limited galaxy catalogue with accurate and precise photometric redshifts, z_ph, derived using the nine photometric bands available in the KiDS-1000 data. This highly pure and complete² galaxy dataset was selected to match the properties of the partly overlapping Galaxy And Mass Assembly (GAMA, Driver et al. 2011) spectroscopic redshift, z_sp, dataset. KiDS-bright is limited to r <  20 mag, covers ∼1000 deg² and contains about one million galaxies after artefact masking. To obtain photometric redshift estimates, Bi21 took advantage of the large amount of spectroscopic calibration data measured by GAMA and trained a supervised machine-learning neural network algorithm implemented in the ANNz2 software (Sadeh et al. 2016) to map an input space of 9-band magnitudes to an output redshift. The ANNz2 training sample consists of matched KiDS galaxies with spectroscopic redshifts from the GAMA equatorial fields, where that survey is the most complete and provides representative training data. The trained model was subsequently applied to the entire inference dataset, the photometrically selected KiDS galaxies with the same r <  20 cut and magnitudes detected in the same nine bands. This sample spans the redshift range of 0 <  z ≲ 0.6; however, since our analytical DSS model is less accurate for very small redshifts, we further exclude all galaxies with z_ph <  0.1. This cut only slightly lowers the number of lenses and results in a projected number density of 0.325 arcmin⁻², as summarised in the first row of Table 1.

The main properties of interest to us are the galaxy bias, which has not been measured before for the KiDS-bright sample, the galaxy number density, and the redshift distribution, which is needed in the modelling. As shown in Bi21, the photometric redshift distribution of the full KiDS-bright sample is measured with high precision: jackknife sub-sampling reveals negligible mean bias, with a small overall scatter of σ_z ≈ 0.018(1 + z). However, for our theoretical model, we need to estimate our newly selected foreground sample’s true n(z) distribution. We take advantage of the very good match between the GAMA spectroscopic sample and the KiDS-bright dataset, allowing us to build an accurate model of the photometric redshift error distribution, as discussed in Bi21. Following the description in Peacock & Bilicki (2018), the best-estimated n_be(z) of the true n(z) can be obtained from a convolution between the normalised photometric redshift distribution of all the galaxies in our selected sample, n_ph(z), and a photo-z error model p_δz(Δz),

$$\begin{array}{c}\hfill n_{\mathrm{be}}(z)={\displaystyle \int n_{\mathrm{ph}}}(z_{\mathrm{ph}})p_{\delta z}(z_{\mathrm{ph}}-z)\mathrm{d}{z}_{\mathrm{ph}}.\end{array}$$(9)

Following Bilicki et al. (2014), we adopted a ‘modified Lorentzian’,

$$\begin{array}{c}\hfill p_{\delta z}(\mathrm{\Delta }z)\propto {(1+\frac{\mathrm{\Delta }{z}^2}{2a{s}^2})}^{-a},\end{array}$$(10)

which was shown in Bi21 to reproduce the photo-z errors better than a Gaussian error model. To estimate the parameters a and s, Eq. (10) is fitted to the KiDS-bright galaxies that also have GAMA spectroscopic redshifts, where Δz = z_ph − z_sp. The best-fit a and s values for our selection are provided in the top row of Table 1, and the resulting n_be(z) is shown in black in Fig. 1.

Fig. 1. Redshift distributions, n(z), of the galaxy samples. The lens sample is obtained from the KiDS-bright galaxies described in Bi21. The blue line shows nsp(z), i.e. the redshift distribution of KiDS galaxies for which we have spectra from the GAMA survey. The red line shows nph(z), the distribution of the full KiDS-bright sample as estimated by ANNz2 with a photometric redshift cut of zph <  0.1. The black line shows our best-estimated nbe(z): a smoothed version of nph(z) that better accounts for photometric redshift errors. The cyan, orange, and brown lines show the third, fourth and fifth redshift bins of the KiDS-1000 data, as estimated in H21. As the third bin strongly overlaps the sources, it is excluded from the analysis.

In order to estimate the uncertainty on our n(z) estimate, we notice that Eq. (9) has two ingredients: the photometric redshift distribution n_ph(z), which is ‘exact’ in the sense that they are directly measured, and the photo-z error model p_δz. Putting aside possible systematic effects related to adopting this machine learning approach in this framework, we only need to quantify the uncertainty associated with our choice of the error model. To account for this, we measured how much the output n(z) changes when the a and s parameters are determined from different sub-samples on the sky. For this, we split the KiDS-bright × GAMA matched sample into ten sub-samples along the right ascension, where each sub-sample has the same number of objects. We fitted a and s to each sub-sample with Eq. (10) and convolve the resulting p_δz with the full p(z_ph) of the KiDS-bright sample. The resulting ten n(z) distributions are almost indistinguishable from our best estimate, as displayed in Fig. A.3. We, therefore, conclude that we can safely neglect the error coming from the Lorentzian fit.

It is difficult to estimate all uncertainties accurately on the n(z) estimate since, for this, we would need to test the ANN2z algorithm on another spectroscopic survey with the same selection. We investigate this further and study the impact of changing the shape and the mean of the n(z). Therefore, besides the best-estimated n_be(z), we also consider using the n_ph(z) directly, as well as the spectroscopic redshift distribution n_sp(z) itself, coming from the matched KiDS-bright × GAMA galaxies; both also shown in Fig. 1. We further allow the n(z) to shift along the redshift direction to give the analysis some flexibility, where the shift value δ⟨z⟩ is drawn from a Gaussian with a standard deviation of 0.01 and vanishing mean, motivated by the uncertainty on the mean of the source redshift distribution.

In addition to the full lens sample described above, we take advantage of the colour information contained in the KiDS-bright data to construct colour-selected sub-samples. This allows us to constrain the bias of blue and red galaxy populations separately from the DSS signal. Following Bi21, we use an empirical split between red and blue galaxies based on their location on the absolute r-band magnitude M_r and the rest-frame u − g colour diagram. The rest-frame quantities are based on LePhare (Arnouts et al. 1999) with the derivations presented in Bi21. We apply a cut through the green valley in the colour-magnitude diagram, which results in a line that delimits the red and blue samples that satisfy

$$\begin{array}{c}\hfill u-g=0.825-0.025M_r.\end{array}$$(11)

We identify those galaxies that are at least 0.05 mag above (below ) the cut line as red (blue) galaxies. We estimate their underlying redshift distributions following the same approach as for the full sample, and the resulting n(z) are shown in Fig. B.1. The effective number densities and best-fit parameters of the modified Lorentzian redshift error model are also listed in the second and third rows of Table 1. We finally note that, while these colour-selected sub-samples are particularly interesting from a galaxy formation perspective, our main cosmological results are obtained from the full lens sample, which has the highest signal-to-noise.

3.2. Source catalogues

The fiducial KiDS-1000 cosmic shear catalogue consists of five tomographic bins, whose redshifts are calibrated using the self-organising map (SOM) method³ of Wright et al. (2020) and presented in Hildebrandt et al. (2021). Although the five bins can be exploited in a cosmic shear analysis as in Asgari et al. (2021), for the DSS analyses, one must also be cautious about source-lens coupling, which arises if sources and lenses belong to the same gravitational potential. This can significantly affect our signal and bias our cosmological inference if left unmodelled. A significant redshift overlap between the source and lens distribution can result in further contamination by the IA of source galaxies that are tidally connected with the foreground lenses. We measure this effect in Sect. 4.2 from IA-infused weak lensing simulations and show that, given the KiDS-bright n(z), this can be avoided by excluding the first three tomographic bins of the KiDS-1000 sources from the analysis. An additional lens-source coupling complication is the so-called boost factor, which arises due to the clustering of sources with over-dense and the anti-correlation of sources with under-dense regions. This can be taken into account by modifying the source n(z) depending on clustering properties (Gruen et al. 2018). Since we exclude the third bin, we do not need to consider this further complication. The fourth and fifth redshift bins, shown as the orange and brown lines in Fig. 1, are separate enough from the lens n(z) to avoid any appreciable lens-source coupling and are therefore used in our cosmological analysis. The uncertainty on the redshift distribution and the residual systematic offsets are very small, as listed in the last two rows of Table 1. The galaxy shear estimates are provided by the lensfit tool (Miller et al. 2013; Fenech Conti et al. 2017) and are described in more detail in Giblin et al. (2021), where it is shown that shear-related systematic effects do not cause more than a 0.1σ shift in S₈ ≡ σ₈(Ω_m/0.3)^0.5 when measured by cosmic shear two-point functions.

4. Simulated data

Besides the real KiDS-1000 data, we validate our inference pipeline on several simulated data sets, study the impact of key systematic uncertainties and carry out the cosmological inference. We use the publicly available FLASK tool (Full-sky Log-normal Astro-fields Simulation Kit) described in Xavier et al. (2016) to estimate the covariance of errors in the DSS data vector of the KiDS-1000; the cosmo-SLICS+IA simulations, described in HD22, to quantify the impact of IA on our measurements and to validate the new N_ap segment in our pipeline (see the end of Sect. 2) that was not present in the B22 model; and the Magneticum lensing simulations, first introduced in Hirschmann et al. (2014), to investigate the impact of stellar and AGN feedback. More details are provided in the following sections.

4.1. FLASK log-normal simulations

Our cosmological inference analysis requires an estimate of the error covariance of the DSS data vector. Since an analytical covariance matrix for the DSS is challenging to compute, we make use instead of an ensemble of log-normal simulations produced with the publicly available FLASK tool⁴ (Xavier et al. 2016). In Hilbert et al. (2011) it is shown that log-normal random fields are a good approximation to the 1-point PDF of the weak lensing convergence and shear field, and Friedrich et al. (2020) show that they are in fact accurate enough to estimate the covariance matrix for higher-order statistics in Stage-III lensing surveys (see their Fig. 4)⁵. Compared to full N-body simulations, FLASK log-normal random fields are computationally cheap to create. The fact that FLASK outputs full-sky maps has the advantage that it can easily be masked to match the footprint of the data, making area re-scaling unnecessary. For the creation of our mock catalogues, we use the cosmological parameters that approximately match current cosmological analyses and fixed the matter density parameter to Ω_m = 0.3, the normalisation of the matter power spectrum to σ₈ = 0.74, the dimensionless Hubble parameter to h = 0.7, the dark energy equation-of-state parameter to w₀ = −1 and the power spectrum power-law index to n_s = 0.97. Furthermore, we provide FLASK with the angular power spectrum of the projected matter density field, the convergence power spectrum for both source bins, and the two matter-lensing cross-spectra. By assuming a flat universe throughout this paper, given the n(z) shown in Fig. 1, and using the PYCCL software package⁶ (Chisari et al. 2019) to get the 3D matter density contrast power spectrum P_δ(ℓ/χ, χ), we calculate the angular power spectrum by use of the Limber-approximated projection (Kaiser 1992) as

$$\begin{array}{c}\hfill C_{i,j}(\ell )={\displaystyle \underset{0}{\overset{\infty }{\int }}}\mathrm{d}\chi \frac{W_i(\chi )W_j(\chi )}{{\chi }^2}{P}_{\delta }(\ell /\chi ,\chi ),\end{array}$$(12)

where i, j are placeholder for either the galaxy or convergence projection, such that $W_g(\chi )=n_{\mathrm{l}}(z[\chi ])\frac{\mathrm{d}z(\chi )}{\mathrm{d}\chi }$ for the lenses with redshift distribution n_l, while for the source with redshift distribution n_s we have instead

$$\begin{array}{c}\hfill W_{\mathrm{s}}(\chi )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}{\displaystyle {\int }_{\chi }^{\infty }}\mathrm{d}{\chi }^{\prime }\frac{\chi ({\chi }^{\prime }-\chi )}{{\chi }^{\prime }a(\chi )}{n}_{\mathrm{s}}(z[{\chi }^{\prime }])\frac{\mathrm{d}z[{\chi }^{\prime }]}{\mathrm{d}{\chi }^{\prime }}.\end{array}$$(13)

Besides these angular power spectra, FLASK needs the log-normal shift parameters κ₀ and δ₀, where −κ₀ and −δ₀ defining the lower limits of the log-normal random variable of the convergence and matter density fields, respectively. Whereas the shift parameter κ₀ = {0.02, 0.03} for the convergence power spectra for the two source bins can be determined directly from the fitting formula Eq. (38) in Hilbert et al. (2011), we estimated the shift parameter δ₀ = {0.57, 0.59, 0.56} for the three lens samples (full, red, blue), as described in Gruen et al. (2016), by assuming that it can be approximated by the shift parameter of the smoothed density contrast, which in turn is calculated from our model for a top-hat filter function (see Eq. (23) in B22). Given this setup, FLASK returns a foreground density map δ_m,2D(θ) and two sets of correlated shear and convergence grids γ_1,2(θ),κ(θ), one per tomographic source bin. We populate our mock KiDS-bright galaxies on the density map by sampling, for each pixel θ, a Poisson distribution with mean parameter λ = n_eff[1+b δ_m,2D(θ)], where b = 1.4 is the (constant) linear galaxy bias estimated from preliminary analyses with only some realisations⁷ and n_eff = 0.325 arcmin⁻² is the mean galaxy density of the KiDS-bright sample. Similarly, we populate the two source planes by Poisson-sampling for each pixel a number of source galaxies n_pix with parameter λ = n_effA_pix, where A_pix is the area of the pixel under consideration⁸, and the effective number density n_eff is taken from Table 1. We finally combine the two shear components of each object with their convergence to construct reduced shear components g_1, 2, and further combine these with a shape noise contribution ϵ^s taken from sampling a Gaussian distribution with vanishing mean and deviation σ_ϵ also taken from Table 1. This results in catalogues containing observed ellipticities ϵ_obs transformed as (Seitz & Schneider 1997)

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{obs}}=\frac{{\mathit{ϵ}}^{\mathrm{s}}+\mathit{g}}{1+{\mathit{ϵ}}^{\mathrm{s},\ast }\mathit{g}}.\end{array}$$(14)

The quantities in bold here are all complex numbers, and the asterisk ‘’ indicates complex conjugation. This procedure ensures that we match the number of foreground and background galaxies in the data and the associated shape noise level.

4.2. cosmo-SLICS+IA

As mentioned earlier, the second suite of simulations is used to validate the inference pipeline and study the impact of IA on our DSS measurements. We use for this the fiducial suite of the cosmo-SLICS presented in Harnois-Déraps et al. (2019), which consists of a set of 50 simulated light cones of 100 deg² each, run in a ΛCDM universe with Ω_m = 0.2905, Ω_Λ = 0.7095, Ω_b = 0.0473, h = 0.6898, σ₈ = 0.836 and n_s = 0.969. The mocks follow the non-linear evolution of 1536³ particles up to z = 0, computed by the cubep³mN-body code (Harnois-Déraps et al. 2013). For Fourier modes of comoving wave number k <  2.0 h Mpc⁻¹, the cosmo-SLICS three-dimensional dark matter power spectrum P(k) agrees within 2% with the predictions from the Extended Cosmic Emulator (Heitmann et al. 2014), followed by a progressive deviation for higher k-modes (Harnois-Déraps et al. 2019), offering a sufficient resolution to model Stage-III galaxy surveys. The particle data were assigned onto mass sheets at 18 redshifts and then post-processed into 10 × 10 deg² light cones. Lensing maps were produced at 18 source redshift planes for each cosmo-SLICS light cone and used to interpolate lensing information onto galaxy catalogues.

Similar to B22, we construct cosmo-SLICS mock source samples that reproduce a number of key data properties, including the tomographic n(z), the galaxy number density n_eff and the shape noise levels. As for the FLASK simulations, we use Eq. (14) to add shape noise to the reduced shear signal. Source galaxies are placed at random positions on the light cones, and the shear quantities (γ_1/2, κ) are interpolated at these positions from the enclosing lensing maps.

We also construct mock KiDS-bright samples by populating the light cone mass maps with galaxies that trace the underlying dark matter field linearly, following the method presented in Harnois-Déraps et al. (2018). We here again fix the galaxy bias to 1.4 and an effective number density of n_eff = 0.325 arcmin⁻².

4.2.1. IA infusion

The impact of galaxy IA is a known secondary signal to the cosmic shear measurements that have been neglected in past DSS studies. In this paper, we verify the validity of this assumption by measuring our statistics in simulated source data that are infused with IA. We find that IA influences our data vector only if the lenses’ n(z) overlap with that of the sources. Following the methods described in HD22, the IA properties of these galaxies are computed as

$$\begin{array}{c}\hfill {ϵ}_1^{\mathrm{IA}}=-\frac{A_{\mathrm{IA}}\overline{C_1}\overline{\rho }(z)}{D_{+}(z)}(s_{\mathit{xx}}-s_{\mathit{yy}}),{ϵ}_2^{\mathrm{IA}}=-\frac{2A_{\mathrm{I}A}\overline{C_1}\overline{\rho }(z)}{D_{+}(z)}{s}_{\mathit{xy}},\end{array}$$(15)

where s_ij = ∂_ijϕ are the Cartesian components of the projected tidal field tensors interpolated at their positions, with ϕ being the gravitational potential. In the above expression, A_IA captures the strength of the coupling between the ellipticities and the tidal field, $\overline{\rho }(z)$ is the matter density, D₊(z) is the linear growth factor, $\overline{C_1}=5\times {10}^{-14}{M}_{\odot }^{-1}{h}^{-2}$ Mpc³, as calibrated in Brown et al. (2002). These intrinsic ellipticity components ${ϵ}_{1/2}^{\mathrm{IA}}$ are then combined with the cosmic shear signal by Eq. (14), resulting in an IA-contaminated weak lensing sample that is consistent with the NLA model of Bridle & King (2007). We refer the reader to HD22 for full details about the IA infusion method. We test several values of A_IA, more precisely, we infused A_IA = {1, 1.5, 2}, and inspect in each case the impact on the DSS data vector.

4.3. Magneticum

Baryon feedback is also known to affect the distribution of the large-scale structure significantly, as the sustained outflows of energy arising from stellar winds, supernovae and AGN reduce the clustering on intra-cluster scales by up to tens of per cent (van Daalen et al. 2011). The exact strength of this suppression is still largely uncertain, with different hydro-dynamical simulations predicting different redshift and scale dependencies (see, e.g. Chisari et al. 2015, for a review of recent results). Without consensus, we opted to measure the DSS in one of these hydro-dynamical simulations for which the impact is quite high and inspect how an extreme baryon impact would affect our data vector.

The Magneticum lensing simulations were first introduced in Hirschmann et al. (2014) and used to mock up KiDS-450 and Stage-IV cosmic shear data (Martinet et al. 2021), and subsequently in Harnois-Déraps et al. (2021) to study the impact of baryons in the peak count analysis of the Dark Energy Survey Y1 data. The underlying matter field is constructed from the Magneticum Pathfinder simulations⁹, more specifically by the Run-2 and Run-2b data described in Hirschmann et al. (2014) and Ragagnin et al. (2017). These are based on the GADGET3 smoothed particle hydrodynamical code (Springel 2005) and are able to reproduce a large number of observations (see Castro et al. 2021, for more details). These both co-evolve dark matter particles of mass 6.9 × 10⁸h⁻¹ M_⊙ and gas particles with mass 1.4 × 10⁸h⁻¹ M_⊙, in comoving volumes of side 352 and 640 h⁻¹ Mpc, respectively. Included key mechanisms are radiative cooling, star formation, supernovae, AGN, and their associated feedback on the matter density field. From sequences of projected mass planes, we use the procedure outlined above for the cosmo-SLICS simulations to generate KiDS-1000 sources and KiDS-bright lenses for ten pseudo-independent light cones, each covering 100 deg². We repeat the same procedure on dark matter-only light cones, such that any difference is caused by the presence of baryons.

The cosmo-SLICS+IA and Magneticum light cones are square-shaped, a geometry that accentuates the edge effects when the aperture filter overlaps with the light cone boundaries. One could, in principle, weight the outer rims for each N_ap map, such that the whole map can be used; although this would increase our statistical power, it could also introduce a systematic offset. We opted instead to exclude the outer rim for each realisation, resulting in an effective area of 36 deg², where a 2 deg band has been removed, matching the size of the adapted filter. This procedure also ensures that roughly the same number of background galaxies are used to calculate the shear profile around each pixel.

5. Cosmological parameter inference

Before performing several Markov chain Monte Carlo (MCMC) samplings in the following two sections, we describe here the pipeline of our Monte-Carlo sampler. In our different MCMC runs, the model vector, the data vector and the covariance matrix are varied, but the overall pipeline stays the same.

Our statistical analysis has the following two free cosmological parameters that we fitted for: the matter density parameter Ω_m and the normalisation of the power spectrum σ₈. We additionally vary the galaxy bias term b and the super-Poisonnian shot-noise parameter α (see Eq. (7)). We detail the prior ranges of all parameters in Table 2, where we also show the Gaussian priors for the nuisance parameters used in the data analysis (but not in the simulation-based validation runs).

For the estimated covariance matrix $\stackrel{\sim }{C}$, which itself is a random variable, Percival et al. (2022) suggested a procedure that uses a more general joint prior of the mean and covariance matrix as the Jeffreys prior proposed in Eq. (6) in Sellentin & Heavens (2016). The method by Percival et al. (2022) leads to credible intervals that can also be interpreted as confidence intervals with approximately the same coverage probability. From a data vector d and a covariance matrix $\stackrel{\sim }{C}$ measured from n_r simulated survey realisations, the posterior distribution of a model vector m that depends on n_θ parameters Θ is

$$\begin{array}{c}\hfill \mathit{P}\left(\mathit{m}(\mathbf{\Theta })|\mathit{d},\stackrel{\sim }{C}\right)\propto {|\stackrel{\sim }{C}|}^{-\frac{1}{2}}{(1+\frac{{\chi }^2}{n_{\mathrm{r}}-1})}^{-m/2},\end{array}$$(16)

where

$$\begin{array}{c}\hfill {\chi }^2={[\mathit{m}(\mathbf{\Theta })-\mathit{d}]}^{\mathrm{T}}{\stackrel{\sim }{C}}^{-1}[\mathit{m}(\mathbf{\Theta })-\mathit{d}].\end{array}$$(17)

The power-law index m is

$$\begin{array}{c}\hfill m=n_{\theta }+2+\frac{n_{\mathrm{r}}-1+B(n_{\mathrm{d}}-n_{\theta })}{1+B(n_{\mathrm{d}}-n_{\theta })}\end{array}$$(18)

with n_d being the number of data points and

$$\begin{array}{c}\hfill B=\frac{n_{\mathrm{r}}-n_{\mathrm{d}}-2}{(n_{\mathrm{r}}-n_{\mathrm{d}}-1)(n_{\mathrm{r}}-n_{\mathrm{d}}-4)}.\end{array}$$(19)

By setting m = n_r the formalism of Sellentin & Heavens (2016) is recovered.

Finally, since the model prediction is too slow for our MCMC, we use the emulation tool contained in CosmoPower (Spurio Mancini et al. 2022), which was first developed to emulate power spectra but can easily be adapted for arbitrary vectors. We trained the emulator on 4000 model points in the parameter space {Ω_m, σ₈, b, α} distributed in a Latin hypercube, where we also included δ⟨z⟩ Gaussian distributed values with the mean as shown in Table 2 but twice the standard deviation¹⁰. To quantify the accuracy of the emulator, we calculated the model at 500 independent points in the same parameter space, as determined with the emulator or directly with the model and show the model vector accuracy in Fig. A.4. The fractional error is better than 2% (95% confidence level).

5.1. Reporting parameter constraints and goodness of fit

In this work, we followed the approach of Joachimi et al. (2021) to report our parameter constraints. In particular, we seek to report the global best fit to the data, that is the set of parameter values that provide the maximum a posteriori (MAP) distribution, computed as

$$\begin{array}{c}\hfill {\mathbf{\Theta }}_{\mathrm{MAP}}=\underset{\mathbf{\Theta }}{\mathrm{argmax}}[\mathit{P}(\mathit{m}(\mathbf{\Theta })|\mathit{d},\stackrel{\sim }{C})],\end{array}$$(20)

where we found the maximum by running several minimisation processes. To estimate the resulting uncertainties around the MAP, we use the suggested projected-joint-highest-posterior-density (PJ-HPD) method, which calculates the parameter ranges that encompass the 68% and 95% credible intervals.

Furthermore, with the degrees of freedom (d.o.f.), we also report the reduced χ²/d.o.f. to quantify the goodness of fit, where the χ² results from the point in the high-dimensional parameter space that has the highest posterior probability. To unbias the covariance matrix $\stackrel{\sim }{C}$, which is used to estimate the χ²-values, we instead of inverting $\stackrel{\sim }{C}\mathrm{h}$ with the known Hartlap factor (Hartlap et al. 2007) defined as h = (n_r − 1)/(n_r − n_d − 2), but rather use

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}^{\prime }=\frac{(n_{\mathrm{r}}-1)[1+B(n_{\mathrm{d}}-n_{\theta })]}{n_{\mathrm{r}}-n_{\mathrm{d}}+n_{\theta }-1}\stackrel{\sim }{C}.\end{array}$$(21)

To estimate the d.o.f. we measure for 1000 mock data vectors the best χ² and fit a χ² distribution to it. The 1000 mock data vectors are drawn from a multivariate Gaussian distribution, where the mean is the model prediction at the MAP values, and the covariance is the corresponding covariance matrix for that particular model. As we have only four free parameters and the rest are fixed by prior knowledge, we expect that the resulting d.o.f. is only slightly less than the raw number of elements in the data vector. Lastly, we report the p-value in each case, which provides the probability of finding a χ² that is more extreme for the given d.o.f., and therefore indicates the goodness of fit. For our analysis, we choose a significance level of 0.01 to be a reliable fit. We have verified with some selected cosmo-SLICS nodes that the theoretical model for the KiDS-bright sample is valid and accurate for σ₈ <  1.0, with indications that it loses accuracy for larger values of σ₈. This does not affect our results, given that the preferred values of σ₈ are well below this limit.

6. Validating the model on simulations

In this section, we validate our model on simulated measurements, several of which are infused with known and controlled systematic effects. The first (fiducial) test establishes that our model is unbiased in the simplest setup, where lens galaxies linearly trace the pure dark matter density maps, while source galaxies are given by the noise-free pure gravitational shear. The second test verifies that our results are unchanged in the presence of IA, as described in Sect. 4.2, while, lastly, we investigate the impact of baryonic physics on our statistics. The fiducial and the IA tests use 20 light cones, which have an unmasked area close to that of the KiDS-1000 footprint¹¹. The measurements on the Magneticum mocks use the ten available light cones. Furthermore, we measured the shear profiles from KiDS-like mocks with shape noise to validate our model on a realistic data vector.

6.1. Validation on intrinsic alignment

To quantify the influence of IA on our results, we performed an MCMC analysis on the mock infused with a strong IA amplitude (A_IA = 2.0) and compared it to our fiducial model (A_IA = 0). For this, we used the mocks described in Sect. 4.2, excluding the third source tomographic bin due to the lens-source coupling. The resulting profiles and mean relative number counts that are used for our pipeline validation are shown in Figs. 2 and 3, respectively. Although the aperture number is not affected by IA, it has a slight effect on the profiles, and hence we verify how it impacts the full data analysis. Although the fourth quantile in N_ap has 2 σ deviation, the p-value is approximately 0.2, which shows that this is consistent with being a statistical fluctuation. We performed this validation test for the adapted and top-hat filters but show the resulting shear signals, and the corresponding mean aperture number values only for the adapted filter since those of the top-hat filter are very similar and would not yield more insights.

Fig. 2. Shear profiles measured with the adapted filter for the KiDS-bright-like lenses and sources from the cosmo-SLICS+IA simulations for two different IA amplitudes (see the legend). The orange regions are estimated from the covariance matrix, while the black dashed lines are obtained from our IA-free analytical predictions at the input cosmology and using the n(z) shown in Fig. 1. As the differences between the shear profiles with varying IA amplitude can barely be seen, we display the relative difference between them in the bottom panels for the highest and lowest two quantiles.

Fig. 3. Relative difference between the mean ⟨Nap⟩ measured in the cosmo-SLICS+IA bright mocks for four quantiles, and the value of ⟨Nap⟩ predicted by the model at the same cosmology. The blue error bars show the KiDS-1000 statistical uncertainty measured from FLASK.

To quantify our decision to discard the third source tomographic bin from our analysis, we show in Fig. A.5 the same shear profiles as in Fig. 2 but also the ones resulting from the third redshift bin, where we clearly see the third tomographic bin is heavily affected by IA, due to a significant overlap in redshift between the lens and source populations, and therefore would need additional modelling of IA, which we disregard for this work, and thus we exclude the third bin.

The MCMC results for the two IA amplitudes (no IA and IA = 2.0) are shown in Fig. 4. First, it is clearly shown that changing the IA amplitude does not affect the posterior at all. Second, and very importantly, the input cosmology is recovered. We observe a small offset on the parameter α, but the other parameters are all recovered within the 1 σ region. This confirms that the < 6% deviations seen on the aperture number count presented in Fig. 4 do not impact our cosmological inference. Finally, we observe an anti-correlation between the S₈ or σ₈ parameter and the galaxy bias b parameter, which is expected since all three parameters are directly correlated with the amplitude of the shear signal. This correlation and the correlation of S₈ and σ₈ to the α parameter could potentially impair the robustness of the later constrained parameters. However, these parameters are particularly important for our model and, therefore, cannot be ignored. The same validation is done for the top-hat filter in Appendix C.

Fig. 4. Pipeline validation: cosmological inference with the adapted filter using the cosmo-SLICS simulations with and without IA infusion, analysed with our model that ignores IA. The posteriors are almost indistinguishable from each other.

6.2. Validation on baryonic feedback

As a last important verification, we investigate for the first time the impact of baryons on the DSS with the Magneticum simulations described in Sect. 4.3. By combining the different DM-only and Hydro mock data for the lenses and sources, we end up with four scenarios (lens-source = DM-DM, DM-Hydro, Hydro-DM and Hydro-Hydro). Figure 5 shows the residuals between the shear profiles measured from dark matter-only mocks (DM-DM) to the other three combinations for each quantile. Clearly, the deviations are well inside the expected KiDS-1000 uncertainty. The biggest differences are seen if baryonic feedback processes are included in the lens mocks: some pixels, close to the N_ap threshold between two quantiles, are shifted to another quantile by the presence of baryons. This is in concordance with the mean aperture numbers reported in Fig. 6, which shows that the mean aperture numbers with baryons are slightly lower. Different to our studies of baryonic feedback such as Heydenreich et al. (2022) or Harnois-Déraps et al. (2021), the inclusion of baryons in the sources has only a minor impact on the DSS, but as expected, becoming more important at small scales. In light of this, we can safely neglect the impact of baryons in our real data analysis.

Fig. 5. Absolute differences between the mean shear profiles for all quantiles using either dark matter-only or full hydrodynamical Magneticum simulations. The residuals are with respect to the dark matter-only mocks for the lenses and sources and are always within the expected statistical uncertainty shown as the orange bands.

Fig. 6. Relative difference between the mean ⟨Nap⟩ to compare measurements from the hydro and dark matter-only Magneticum simulations. The relative difference is always well inside the expected statistical uncertainty of KiDS-1000.

7. Results and discussion

After validating our model to simulations in B22 and our additional testing on the impact of IA and baryonic physics, we are well equipped to analyse real lensing data accurately. As described in Sect. 3, we use the KiDS-bright sample as our foreground lenses and the fourth and fifth KiDS-1000 tomographic bins as our cosmic shear data.

To recap, in our fiducial analysis, we used the n(z) shown as the black solid line in Fig. 1, we varied the two cosmological parameters Ω_m and σ₈ as well as the two astrophysical parameters b and α. We marginalised over the systematic effects parameters describing the δ⟨z⟩ and m-bias uncertainty. In Figs. 7 and 8 we display the resulting shear profiles and mean aperture number. The model shown in these figures is computed at the best-fit MAP values, listed in the first column of Table 3. In that table, the p-values indicate that the data are well fitted by the model, being all well above our threshold value fixed at p = 0.01. The d.o.f. are estimated as described in Sect. 5, where we show in Fig. A.6 the distribution of χ² values and for which a χ²-distribution with 81 d.o.f. fits well. As expected, the resulting d.o.f. is slightly lower than the raw number of elements. The reduced χ² values are slightly below the expectation of 1.0, potentially indicating that the uncertainties could be slightly overestimated, although they are well inside the expected reduced χ² scatter of $\pm \sqrt{2}$, hence do not warrant further investigation.

Fig. 7. Shear profiles measured in the data, compared to the best-fit predictions (MAP) obtained with values listed in Table 3, for the adapted filter. The shaded region shows the statistical uncertainty estimated from 1000 FLASK realisations.

Fig. 8. Relative difference between the mean ⟨Nap⟩ to compare measurements from the KiDS-bright sample and our model, evaluated at the MAP values shown in Table 3. The measured ⟨Nap⟩ are all greater than the predicted ⟨Nap⟩ at MAP just indicating that the measured p(Nap) is broader than the predicted one.

Using the approach described in Sect. 5 to estimate the uncertainty around the MAP, we find

$$\begin{array}{c}\hfill S_8^{\mathrm{DSS}}=0.{731}_{-0.018}^{+0.030},\end{array}$$(22)

which is consistent with and competitive to the KiDS-1000 cosmic shear constraints from Asgari et al. (2021),

$$\begin{array}{c}\hfill S_8^{\mathrm{COSEBIs}}=0.{759}_{-0.021}^{+0.024}.\end{array}$$(23)

We present the posterior of these two analyses in Fig. 9, where the consistency between the two probes is obvious. The DSS has a slightly lower constraining power on S₈; however, the Ω_m–σ₈ degeneracy is broken, thanks to the additional information provided by the foreground data. But even for the shear-only case shown in green in Fig. 9, the DSS has better constraining power for the Ω_m and σ₈ parameters, although the lower bound should be taken with caution as we excluded all Ω_m <  0.2, as the model does not agree with the cosmo-SLICS for smaller Ω_m values. Furthermore, although the results of Fig. A.2 show that the inferred S₈ might be smaller compared to the truth if a linear galaxy bias model is not sufficient, the consistency between the shear-only DSS analysis to the complete DSS analysis supports the robustness of our inferred parameters with respect to the galaxy bias model (see the discussion at the end of Sect. 2.2). The comparison to the COSEBIs analysis reveals competitive S/N for the S₈ parameter while using only a fraction of the lensing sources (tomographic bins four and five). Caution should be taken when comparing their respective constraining power, as the COSEBIs analysis marginalises over more cosmological parameters and samples the parameter space differently. Nevertheless, it seems that a joint COSEBIs-DSS analysis could further improve the constraints, which we leave for future work. Lastly, as the DSS estimates of S₈ are slightly lower but with higher uncertainty, we measure a similar tension to the CMB results as the COSEBIs analysis.

Fig. 9. Cosmological posteriors for the adapted filter for the best-estimated nbe(z) in black using the full data vector and in green using only shear information compared to the COSEBIs posteriors in orange presented in Asgari et al. (2021) and to the (TT, TE, EE+lowE) results of Planck Collaboration V (2020) in red. The sharp cut of the green posterior is due to the conservative prior of Ωm >  0.2, as the model does not agree with the cosmo-SLICS for smaller Ωm values. The shear-only posterior is shown in this figure to support the assumption of using a linear galaxy bias model.

In the next section, we investigate the robustness of our results with respect to the lens redshift distribution n(z), of varying the covariance matrix. We further present our galaxy colour-split analysis and additionally discuss the galaxy bias b and α results.

7.1. Impact of lens redshift distribution

To estimate the impact of the shape of the lens galaxy redshift distribution (on top of shifting the mean by δ⟨z⟩ in the sampling), we repeat the analysis for the three n(z) shown in Fig. 1. These are the smoothed version of the photometric redshift distribution n_be(z), the photometric redshift distribution n_ph(z) itself without any smoothing, and the spectroscopic redshift n_sp(z) from those GAMA galaxies that are also in the KiDS-bright sample. Although Bi21 showed that GAMA is representative and that mismatches should be rare, the results from the GAMA spectroscopic n_sp(z) should be taken with caution because the equatorial fields have a relatively small sky coverage, leading to features in the n_sp(z) that are caused by the large-scale structure present in these fields. For this investigation, we use the same setup as for the fiducial analysis, varying the two cosmological parameters Ω_m and σ₈ together with the α and the linear galaxy bias b parameter. We also marginalised over the nuisance parameters shown bottom half of Table 2. In Fig. 10 we display the different posteriors following from the three alternatives n(z). It is clearly seen that the posteriors are shifted along the Ω_m − σ₈ degeneracy axis, whereas these shifts partially cancel out for S₈. Due to the different lens n(z), the amplitude and the slope of the shear signal predictions are slightly different. In particular, higher amplitude and steeper slope of the shear profiles result in larger σ₈ and smaller Ω_m, and vice-versa. Furthermore, we notice that the linear galaxy bias b and the noise α are stable against changes in the n(z).

Fig. 10. Posteriors of the full DSS vector resulting from using the different lens n(z) shown in Fig. 1. The n(zphot) and n(z) have, by construction, the same mean redshift, while the mean redshift of the spectroscopic redshift estimate is ∼0.015 lower. The main DSS results include a marginalisation over our uncertainty in the mean redshifts of both the lens and source samples, which partially compensates for some of these differences. The posterior obtained using the GAMA spectroscopic n(z), shown in blue, should be taken with caution as it is estimated from a smaller sky coverage and, as such, contains larger statistical fluctuations.

Lastly, in order to investigate the impact of our photometric redshift cut at z_ph = 0.1 (see Sect. 3.1), we perform two additional analyses, this time modifying the photometric redshift threshold to z_ph >  0.15 and z_ph >  0.2. We find that the posteriors shifts are smaller than the 68% credibility region, indicating that removing low-redshift galaxies does not result in systematically different Ω_m or σ₈.

7.2. Cosmology scaling of the covariance matrix

The covariance matrix used in the main analysis is determined at a specific point in the parameter space (see Sect. 4.1), which is not identical to the MAP. However, assuming the MAP values are the true parameters, the most robust likelihood analysis would be achieved with a data covariance matrix estimated at the MAP values. This section, therefore, explores the impact of considering a cosmology-dependent covariance matrix.

In the related literature, there are two common approaches on whether the cosmology should be kept fixed in the covariance matrix (van Uitert et al. 2018) or varied at each point sampled by the MCMC, as in Eifler et al. (2009). It is argued in Carron (2013) that the latter could result in over-constraints, whereas Kodwani et al. (2019) argues that the effect is small. We, therefore, explore both methods here.

We achieve our cosmology rescaling by assuming that the covariance matrix scales quadratically with the signal. This is only strictly true in the Gaussian regime; nevertheless, it is a good first approximation, even though the impact on the non-linear mode coupling is neglected in this approach.

To achieve the rescaling, we compute at a new cosmology Θ the ratio between the predicted model m(Θ) and the model predicted at the FLASK cosmology m(Θ^F),

$$\begin{array}{c}\hfill r_i(\mathbf{\Theta })=\frac{m_i(\mathbf{\Theta })}{m_i\left({\mathbf{\Theta }}^{\mathrm{F}}\right)}.\end{array}$$(24)

We then multiply each element of the fiducial covariance matrix, $C_{\mathrm{ij}}^F$, by the scaling factors,

$$\begin{array}{c}\hfill C_{\mathit{ij}}(\mathbf{\Theta })=C_{\mathit{ij}}^{\mathrm{F}}{r}_i(\mathbf{\Theta })r_j(\mathbf{\Theta }),\end{array}$$(25)

and obtain a cosmology-rescaled covariance matrix.

As explained in Eifler et al. (2009), this method is only valid for a noise-free covariance matrix since it wrongly rescales the shape-noise component and possibly over-estimates the cosmology changes. Finally, we note that in Eq. (16) the determinant of the covariance matrix changes with cosmology as well and needs to be recalculated.

Using our fiducial setup, we determine the posterior distribution in two distinct ways: first, by varying the covariance matrix alongside the model vector at each step of the MCMC, and second, by scaling the covariance matrix to the MAP value. For the latter approach, we use an iterative process, where we first estimate the MAP with the fiducial covariance matrix, then use the MAP parameters to scale the covariance matrix and find new MAP parameters; we repeated that process 100 times. As seen in Fig. A.7 after approximately 20 iterations the MAP values converged. The results are shown in Fig. 11, where the red posteriors used a covariance rescaled to the converged MAP value; the blue posteriors are for a full parameter-dependent covariance matrix varied in the MCMC; the black posteriors show the fiducial covariance. The red and black posteriors are almost identical, which is not surprising given the fact that the MAP values are very close to the parameters used to determine the covariance matrix. However, the blue contours slightly differ from the other two but are still within half 1 σ; we are therefore not concerned about the impact on our constraints of this analysis choice. In Stage IV surveys, this is even less important due to the tighter posterior and the, therefore, smaller cosmology variation.

Fig. 11. Comparison between the posteriors obtained from our fiducial setup (with the covariance matrix calculated at the initial cosmology, see the black contours), with those obtained after scaling the covariance at the best-fit parameters (red) and to those obtained by varying the covariance with the parameters (blue).

7.3. Red and blue split

In this section, we present our final investigation, dividing the KiDS-bright sample into red and blue galaxies according to their colour as described in Sect. 3.1, and carry out a joint analysis. The motivation for this is to learn more about the behaviour of different galaxy types and as a cosmological robustness check. As for the main analysis, we use the best-estimated n_be(z), resulting from smoothing the photometric redshifts after applying a z_ph >  0.1 cut. In this setup, our data vector has 168 elements. But given the likelihood modelling described in Sect. 5, we are still confident in our results with respect to the remaining noise in the covariance matrix. As shown in Table 4, the reduced χ²/d.o.f. = 1.00 indicates a valid fit. The resulting posteriors are shown in Fig. 12 and the MAP values are stated in Table 4.

Fig. 12. Posterior for the full KiDS-bright sample shown in green, while the joint red+blue posteriors represent the results of the colour-selected samples. In the latter case, the red and blue samples share the same cosmology (the dark blue contours) by construction. The resulting measured and best-fit predicted shear profiles for the red and blue samples are displayed in Fig. B.2 and the corresponding mean aperture number values are seen in Fig. B.3.

First, we notice that the cosmological parameters of the full KiDS-bright sample analysis and the joint red and blue analysis are consistent and within 0.48 σ in σ₈ and almost identical in Ω_m. Of particular interest in this investigation are the results obtained for the two astrophysical parameters, where we see that, as expected, the blue (red) sample prefers a smaller (larger) linear galaxy bias b compared to the full sample. This is in line with the fact that red galaxies are known to be more strongly clustered than blue galaxies and therefore have a larger galaxy bias (Mo et al. 2010). Furthermore, we find that the α parameter for the blue sample is significantly larger than unity, whereas the red sample tends to be below one. This shows that blue galaxies follow a super-Poisson distribution and red galaxies a sub-Poisson distribution. Following the results of Friedrich et al. (2022) who found that a higher satellite fraction leads to a higher α value, the blue sample has more satellites. The full sample overlaps with the blue and red posteriors and is consistent with a normal-Poisson distribution (α = 1.0).

8. Summary and conclusions

In this work, we present an unblinded density split statistic analysis of the fourth data release of the Kilo-Degree survey (KiDS-1000). The analytical model used to infer cosmological and astrophysical parameters was first developed in Friedrich et al. (2018) and then modified in Burger et al. (2022), and we further validated it on realistic simulated data. The lenses used to construct the foreground density map are taken from the KiDS-bright sample described in Bi21, while for our sources, we used the fourth and fifth tomographic redshift bins of the KiDS-1000 data described in H21.

We investigated for the first time the impact of baryons and IA on the DSS. While the effect of the former is suppressed due to the implied smoothing of the density map, IA can have an important role if the redshift distributions of the lenses and sources overlap. We carried out a full analysis on contaminated mock data without overlapping redshift distributions and found that for our selected data, we are immune to both of the systematic effects at the level of the inferred posteriors.

We explored the uncertainty of the redshift distribution of the lenses by investigating the impact on the posterior of changing the mean and the shape of the n_ph(z). In particular, for this, we used the photometric redshift distribution n_ph(z) with and without smoothing, as well as the distribution obtained directly from the overlapping spectroscopic GAMA galaxies. We found that the posteriors vary by less than ∼0.5 σ. Notably, we observed that of all parameters, Ω_m is the most affected and is generally lower for broader n(z); in contrast, σ₈ increased, leaving S₈ values changed by ∼0.5 σ. We assigned an extra error term to this uncertainty, resulting in $S_8=0.{73}_{-0.02}^{+0.03}\pm 0.01$ for the n(z) shape after marginalising over the other systematic effects. These constraints are competitive and consistent with the KiDS-1000 cosmic shear constraints from Asgari et al. (2021).

Furthermore, we investigated the impact of varying the covariance matrix with cosmological parameters, where we used an iterative process once to scale the covariance matrix to the MAP best-fit parameters, and we varied them alongside the parameters in the MCMC process once as well; for all cases, we recorded no significant deviations.

As our final result, we divided the full KiDS-bright sample into red and blue galaxies as a cosmology robustness check and to learn more about different galaxy types. For this, we performed a joint analysis of the red and blue samples with a joint covariance matrix with the smoothed version of n_ph(z) as our best redshift estimate. The resulting posteriors of the full and joint red+blue analyses agree as to the cosmological parameters within 0.35 σ. Furthermore, this shows the expected behaviour of the linear galaxy bias, where blue (red) galaxies have a lower (higher) bias than the full sample. The α parameter, which accounts for super-Poisson or sub-Poisson shot noise, also revealed interesting results. Whereas red galaxies have an α value that tends to be smaller than one (≈2 σ), blue galaxies have an α value significantly larger than one (≈6 σ), meaning that blue galaxies are super-Poisson distributed and red galaxies are sub-Poisson distributed. According to Friedrich et al. (2020), this reveals information about the halo occupation distribution, with samples with a larger fraction of satellite galaxies tending to have larger α values.

We conclude from our results that the density split statistic is a valuable tool with a major advantage in the Ω_m-σ₈ degeneracy breaking. In addition to this, it also yields a new way to measure the galaxy bias on linear scales and the Poissanity of different galaxy types. We save the inference of the dark-energy-equation of state w for future analysis; this can be achieved by a tomographic analysis of high-precision lensing and clustering data. Other aspects for future analysis include the modelling of a more complex galaxy bias model, and lastly, the impact of the filter size, where we expect smaller filter sizes to be more constraining yet also more sensitive to non-linear scales and baryonic analysis.

Acknowledgments

We thank the anonymous referee for the very constructive and fruitful comments. This paper went through the KiDS review process, where we especially want to thank the KiDS-internal referee Benjamin Joachimi for his fruitful comments to improve this work. Further, we would like to thank Mike Jarvis for maintaining treecorr, and Alessio Mancini to develop the CosmoPower emulator, which improved our work significantly. Lastly, we thank Sven Heydenreich, Laila Linke, Patrick Simon and Jan Luca van den Busch for very valuable discussions. PAB acknowledges support by the German Academic Scholarship Foundation. JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). OF gratefully acknowledges support by the Kavli Foundation and the International Newton Trust through a Newton-Kavli-Junior Fellowship and by Churchill College Cambridge through a postdoctoral By-Fellowship. MB is supported by the Polish National Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698, 2018/31/G/ST9/03388 and 2020/39/B/ST9/03494, and by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12. HH is supported by a Heisenberg grant of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1) as well as an ERC Consolidator Grant (No. 770935). AHW is supported by a European Research Council Consolidator Grant (No. 770935). TC is supported by the INFN INDARK PD51 grant and the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA. KD acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679 and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. The calculations for the Magneticum simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr83li and with the support by M. Petkova through the Computational Center for Particle and Astrophysics (C2PAP). CH acknowledges support from the European Research Council under grant number 647112, and 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. BJ acknowledges support by STFC Consolidated Grant ST/V000780/1. KK acknowledges support from the Royal Society, and Imperial College. NM acknowledges support from the Centre National d’Études Spatiales (CNES) fellowship. HYS acknowledges the support from CMS-CSST-2021-A01 and CMS-CSST-2021-B01, NSFC of China under grant 11973070, the Shanghai Committee of Science and Technology grant No.19ZR1466600 and Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013. TT acknowledges support from the Leverhulme Trust. Author contributions: all authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (PAB, OF, JHD, PS), where we ordered (OF, JHD, PS) alphabetically. PAB led the paper, JHD provided all necessary simulations, OF and PS contributed to the modelling, and all four helped in the development of the analysis. The first author group is followed by two further 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: Additional plots

Fig. A.1. Residual shear profiles of the highest quantile filter for different threshold values that the effective area must exceed to be used with respect to the shear profiles with a 0.5 threshold. It is seen that for threshold values above 0.6 the shear profiles deviate more than one standard deviation (grey band) due to decreasing remaining area.

Fig. A.2. Posterior distributions, using different b2 values to distribute galaxies to the same density contrast, from which a data vector is measured. The lower right panel shows the posterior results from simulations, where galaxies are distributed with a HOD. The model vector uses for all four cases a linear galaxy bias model.

In this section, we show additional plots which belong to the main text. We start by visualising the shear profiles that result in different threshold values for the effective area above which pixels are used for the analysis in Fig. A.1.

In Fig. A.2 we show different posterior distributions, to test whether our analytical model with a linear galaxy bias model gives robust cosmological results even if galaxies are placed with a Poisson process with mean $\lambda =n_{\mathrm{eff}}[1+b_1{\delta }_{\mathrm{m},2\mathrm{D}}(\mathit{\theta })+b_2({\delta }_{\mathrm{m},2\mathrm{D}}^2(\mathit{\theta })-⟨{\delta }_{\mathrm{m},2\mathrm{D}}^2)⟩]$, where δ_m, 2D is the projected density contrast or with a halo occupation distribution (HOD) description (Smith et al. 2017). For the Poisson process test, we use the simulation of Takahashi et al. (2017) with sources that mimic the fourth and fifth KiDS-1000 bins and lenses that mimic the KiDS-bright sample. For the HOD analysis, we measure the data vector, and the covariance from 614 SLICS realisations with noisy KiDS-1000 sources and GAMA lens mocks (see Harnois-Déraps et al. 2018, for a detailed description), where the covariance is scaled to approximately match the KIDS-1000 footprint. It is clearly seen that the posterior for b₂ >  0 is strongly biased if shear and N_ap information are used. In contrast, using only shear information, the posterior is unbiased as they are almost insensitive to the galaxy bias model. Although the HOD analysis already shows that the linear galaxy bias assumption is sufficient, if the posterior using shear and N_ap information is consistent with the posterior using only shear information gives additional confidence that b₂ ≈ 0.

Fig. A.3. Redshift distribution of the full KiDS-bright sample resulting from Lorentzian fitting parameters (a and s), which in turn are determined from different patches on the sky. In the bottom panel, the absolute differences to the best-estimated nbe(z) are shown.

Fig. A.4. Difference between the true and emulated model data vector m(Θ) scaled by the standard deviation of the measured data estimated with the FLASK.

Fig. A.5. Shear profiles of the cosmo-SLICS for the bright lens and all three source mocks with the n(z) shown in Fig. 1 for two different IA amplitudes for the adapted filter. The third source tomographic bin is strongly contaminated by IA due to the significant overlap with the lens n(z) and is therefore excluded from this analysis.

In Fig. A.3 we show the negligible difference of the lens n(z) if the Lorentzian fitting parameters (a and s) are estimated from different patches of the sky. In Fig. A.4 we display the accuracy of the emulator, and in Fig. A.5, that the third tomographic bin is indeed contaminated by IA, in Fig. A.6 the χ² distribution of mock data vectors around the MAP for the adapted filter with the best-estimated n_be(z), and lastly in Fig. A.7 the iterative process to find the optimal MAP values by scaling the covariance matrix to the previously found MAP.

Fig. A.6. Distribution of χ2 values of mock data vectors that follow from multivariate Gaussian distributions, where the mean is the model prediction at MAP and the covariance is the corresponding covariance for that particular model. The red line shows the χ2 values using the real data vector. The orange line is a χ2 distribution with d.o.f. = 81, which is slightly smaller than the 84 elements of the data vector.

Fig. A.7. Change of MAP values due to scaling of the covariance to the previously measured MAP values. Roughly after ten iterations, the process converged, where the occasional outliers happened if the minimisation process stopped to early by coincidence.

Appendix B: Additional material for the red and blue analysis

Here in this chapter, we show the complementary plots for the joint red+blue analysis. In Fig. B.1 the redshift distribution for the red and blue samples is shown, resulting from the smoothing method described in Sect. 3.1. In Fig. B.2 we show the shear profiles; in Fig. B.3 the mean aperture number values for both samples which are used as the model and data vectors.

Fig. B.1. Best estimated redshift distributions of the red and blue KiDS-bright samples in red and blue, with the full KiDS-bright sample in green.

Fig. B.2. Measured and predicted shear profiles at the MAP values for the adapted filter for the red and blue samples. The shaded region is the expected KiDS-1000 uncertainty estimated from the 1000 FLASK realisations with shape noise.

Fig. B.3. Mean aperture number of the red and blue sample, where the predictions follow from a joint minimisation process. Different to the full bright sample the measured p(Nap) is smaller than the predicted Map, resulting in the measured ⟨Nap⟩ being smaller.

Appendix C: Top-hat filter analysis

Fig. C.1. MCMC results for the top-hat filter using the model and simulations where we infused the pure shear signal with different amplitudes of IA.

This section shows the corresponding plots that belong to the analysis with the top-hat filter. We start by showing in Fig. C.1 the validation for the top-hat with the cosmo-SLICS, which shows that the true parameters are always inside 1 σ. The discussion from Sect. 7.1 about the impact of different redshifts distributions is summarised for the top-hat filter in Table C.1, which shows with the given reduced χ² and corresponding p-value that the model is well fitted to the data given the covariance matrix, and result in parameters that are consistent to the ones constrained with the adapted filter.

Finally we show in Fig. C.2 and Table C.2 the analogous results for the top-hat filter to the adapted filter as shown in Fig. 12 and Table 4. The p-value for the joint red+blue indicates that the given d.o.f. has a significant tension between the measured data and the best fit model. Although this reduced χ² is not ideal for the given d.o.f., we still perform the analysis, but the posteriors should be taken with caution, which is already true because we are uncertain about the true n(z) of the sub-samples. Overall the results show the same trends as for the adapted filter, where the blue galaxies result in smaller bias b and larger α than the red galaxies. The reason for the top-hat filter performing worse than the adapted filter is unclear. Nevertheless, besides the fact that the n(z) of the red and blue sample is quite uncertain, both filters are probing on fundamentally different scales so that different behaviours are not surprising.

Fig. C.2. Posterior for the full KiDS-bright sample shown in green, and the joint red+blue posteriors that represent the results of the individual colour-selected samples. In the latter case, by construction, the red and blue samples share the same cosmology (the dark blue contours).

All Tables

All Figures

Fig. 1. Redshift distributions, n(z), of the galaxy samples. The lens sample is obtained from the KiDS-bright galaxies described in Bi21. The blue line shows nsp(z), i.e. the redshift distribution of KiDS galaxies for which we have spectra from the GAMA survey. The red line shows nph(z), the distribution of the full KiDS-bright sample as estimated by ANNz2 with a photometric redshift cut of zph <  0.1. The black line shows our best-estimated nbe(z): a smoothed version of nph(z) that better accounts for photometric redshift errors. The cyan, orange, and brown lines show the third, fourth and fifth redshift bins of the KiDS-1000 data, as estimated in H21. As the third bin strongly overlaps the sources, it is excluded from the analysis.

In the text

Fig. 2. Shear profiles measured with the adapted filter for the KiDS-bright-like lenses and sources from the cosmo-SLICS+IA simulations for two different IA amplitudes (see the legend). The orange regions are estimated from the covariance matrix, while the black dashed lines are obtained from our IA-free analytical predictions at the input cosmology and using the n(z) shown in Fig. 1. As the differences between the shear profiles with varying IA amplitude can barely be seen, we display the relative difference between them in the bottom panels for the highest and lowest two quantiles.

In the text

	Fig. 3. Relative difference between the mean ⟨Nap⟩ measured in the cosmo-SLICS+IA bright mocks for four quantiles, and the value of ⟨Nap⟩ predicted by the model at the same cosmology. The blue error bars show the KiDS-1000 statistical uncertainty measured from FLASK.
In the text

	Fig. 4. Pipeline validation: cosmological inference with the adapted filter using the cosmo-SLICS simulations with and without IA infusion, analysed with our model that ignores IA. The posteriors are almost indistinguishable from each other.
In the text

	Fig. 5. Absolute differences between the mean shear profiles for all quantiles using either dark matter-only or full hydrodynamical Magneticum simulations. The residuals are with respect to the dark matter-only mocks for the lenses and sources and are always within the expected statistical uncertainty shown as the orange bands.
In the text

	Fig. 6. Relative difference between the mean ⟨Nap⟩ to compare measurements from the hydro and dark matter-only Magneticum simulations. The relative difference is always well inside the expected statistical uncertainty of KiDS-1000.
In the text

	Fig. 7. Shear profiles measured in the data, compared to the best-fit predictions (MAP) obtained with values listed in Table 3, for the adapted filter. The shaded region shows the statistical uncertainty estimated from 1000 FLASK realisations.
In the text

	Fig. 8. Relative difference between the mean ⟨Nap⟩ to compare measurements from the KiDS-bright sample and our model, evaluated at the MAP values shown in Table 3. The measured ⟨Nap⟩ are all greater than the predicted ⟨Nap⟩ at MAP just indicating that the measured p(Nap) is broader than the predicted one.
In the text

Fig. 9. Cosmological posteriors for the adapted filter for the best-estimated nbe(z) in black using the full data vector and in green using only shear information compared to the COSEBIs posteriors in orange presented in Asgari et al. (2021) and to the (TT, TE, EE+lowE) results of Planck Collaboration V (2020) in red. The sharp cut of the green posterior is due to the conservative prior of Ωm >  0.2, as the model does not agree with the cosmo-SLICS for smaller Ωm values. The shear-only posterior is shown in this figure to support the assumption of using a linear galaxy bias model.

In the text

Fig. 10. Posteriors of the full DSS vector resulting from using the different lens n(z) shown in Fig. 1. The n(zphot) and n(z) have, by construction, the same mean redshift, while the mean redshift of the spectroscopic redshift estimate is ∼0.015 lower. The main DSS results include a marginalisation over our uncertainty in the mean redshifts of both the lens and source samples, which partially compensates for some of these differences. The posterior obtained using the GAMA spectroscopic n(z), shown in blue, should be taken with caution as it is estimated from a smaller sky coverage and, as such, contains larger statistical fluctuations.

In the text

	Fig. 11. Comparison between the posteriors obtained from our fiducial setup (with the covariance matrix calculated at the initial cosmology, see the black contours), with those obtained after scaling the covariance at the best-fit parameters (red) and to those obtained by varying the covariance with the parameters (blue).
In the text

Fig. 12. Posterior for the full KiDS-bright sample shown in green, while the joint red+blue posteriors represent the results of the colour-selected samples. In the latter case, the red and blue samples share the same cosmology (the dark blue contours) by construction. The resulting measured and best-fit predicted shear profiles for the red and blue samples are displayed in Fig. B.2 and the corresponding mean aperture number values are seen in Fig. B.3.

In the text

	Fig. A.1. Residual shear profiles of the highest quantile filter for different threshold values that the effective area must exceed to be used with respect to the shear profiles with a 0.5 threshold. It is seen that for threshold values above 0.6 the shear profiles deviate more than one standard deviation (grey band) due to decreasing remaining area.
In the text

	Fig. A.2. Posterior distributions, using different b2 values to distribute galaxies to the same density contrast, from which a data vector is measured. The lower right panel shows the posterior results from simulations, where galaxies are distributed with a HOD. The model vector uses for all four cases a linear galaxy bias model.
In the text

	Fig. A.3. Redshift distribution of the full KiDS-bright sample resulting from Lorentzian fitting parameters (a and s), which in turn are determined from different patches on the sky. In the bottom panel, the absolute differences to the best-estimated nbe(z) are shown.
In the text

	Fig. A.4. Difference between the true and emulated model data vector m(Θ) scaled by the standard deviation of the measured data estimated with the FLASK.
In the text

	Fig. A.5. Shear profiles of the cosmo-SLICS for the bright lens and all three source mocks with the n(z) shown in Fig. 1 for two different IA amplitudes for the adapted filter. The third source tomographic bin is strongly contaminated by IA due to the significant overlap with the lens n(z) and is therefore excluded from this analysis.
In the text

Fig. A.6. Distribution of χ2 values of mock data vectors that follow from multivariate Gaussian distributions, where the mean is the model prediction at MAP and the covariance is the corresponding covariance for that particular model. The red line shows the χ2 values using the real data vector. The orange line is a χ2 distribution with d.o.f. = 81, which is slightly smaller than the 84 elements of the data vector.

In the text

	Fig. A.7. Change of MAP values due to scaling of the covariance to the previously measured MAP values. Roughly after ten iterations, the process converged, where the occasional outliers happened if the minimisation process stopped to early by coincidence.
In the text

	Fig. B.1. Best estimated redshift distributions of the red and blue KiDS-bright samples in red and blue, with the full KiDS-bright sample in green.
In the text

	Fig. B.2. Measured and predicted shear profiles at the MAP values for the adapted filter for the red and blue samples. The shaded region is the expected KiDS-1000 uncertainty estimated from the 1000 FLASK realisations with shape noise.
In the text

	Fig. B.3. Mean aperture number of the red and blue sample, where the predictions follow from a joint minimisation process. Different to the full bright sample the measured p(Nap) is smaller than the predicted Map, resulting in the measured ⟨Nap⟩ being smaller.
In the text

	Fig. C.1. MCMC results for the top-hat filter using the model and simulations where we infused the pure shear signal with different amplitudes of IA.
In the text

	Fig. C.2. Posterior for the full KiDS-bright sample shown in green, and the joint red+blue posteriors that represent the results of the individual colour-selected samples. In the latter case, by construction, the red and blue samples share the same cosmology (the dark blue contours).
In the text