Issue 
A&A
Volume 669, January 2023



Article Number  A69  
Number of page(s)  18  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202244673  
Published online  10 January 2023 
KiDS1000 cosmology: Constraints from density split statistics
^{1}
ArgelanderInstitut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
email: pburger@astro.unibonn.de
^{2}
UniversitätsSternwarte, Fakultät für Physik, LudwigMaximiliansUniversität München, Scheinerstr.1, 81679 München, Germany
^{3}
School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne NE1 7RU, UK
^{4}
E.A Milne Centre, University of Hull, Cottingham Road, Hull HU6 7RX, UK
^{5}
Center for Theoretical Physics, Polish Academy of Sciences, al. Lotników 32/46, 02668 Warsaw, Poland
^{6}
Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany
^{7}
INAF – Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34143 Trieste, Italy
^{8}
INFN – Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy
^{9}
Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy
^{10}
Germany MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStraße 1, 85741 Garching, Germany
^{11}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
^{12}
Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
^{13}
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
^{14}
Laboratoire d’Astrophysique de Marseille, AixMarseille Univ., Rue Frédéric JoliotCurie 38, 13388 Marseille, France
^{15}
Shanghai Astronomical Observatory (SHAO), Nandan Road 80, Shanghai 200030, PR China
^{16}
University of Chinese Academy of Sciences, Yuquan Road 19A, Beijing 100049, PR China
^{17}
Institute for Particle Physics and Astrophysics, ETH Zürich, WolfgangPauliStrasse 27, 8093 Zürich, Switzerland
Received:
3
August
2022
Accepted:
28
October
2022
Context. Weak lensing and clustering statistics beyond twopoint functions can capture nonGaussian information about the matter density field, thereby improving the constraints on cosmological parameters relative to the mainstream methods based on correlation functions and power spectra.
Aims. This paper presents a cosmological analysis of the fourth data release of the Kilo Degree Survey based on the density split statistics, which measures the mean shear profiles around regions classified according to foreground densities. The latter is constructed from a bright galaxy sample, which we further split into red and blue samples, allowing us to probe their respective connection to the underlying dark matter density.
Methods. We used the stateoftheart model of the density splitting statistics and validated its robustness against mock data infused with known systematic effects such as intrinsic galaxy alignment and baryonic feedback.
Results. After marginalising over the photometric redshift uncertainty and the residual shear calibration bias, we measured for the full KiDSbright sample a structure growth parameter of that is competitive and consistent with twopoint cosmic shear results, a matter density of Ω_{m} = 0.27 ± 0.02, and a constant galaxy bias of b = 1.37 ± 0.10.
Key words: cosmological parameters / largescale structure of Universe / gravitational lensing: weak / methods: statistical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
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 twopoint statistics, namely the twopoint 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 Gaussiandistributed, such as the cosmic microwave background (CMB; e.g. Planck Collaboration V 2020). In the late Universe, however, nonlinear gravitational instabilities generate a significant amount of nonGaussian features, whose information can only be accessed with higherorder statistics. Furthermore, since higherorder 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 higherorder 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_{0} (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 , 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 σ_{8} − Ω_{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, twopoint statistics have been favoured so far for the reasons mentioned above, such that the combination of all measurements (cosmic shear, galaxy clustering, and galaxygalaxy lensing) are generally referred to as ‘3×2pts’ statistics.
To access the additional information contained in the nonGaussian features, a competitive statistic to the 3×2pts method was recently proposed in Gruen et al. (2016), coined the ‘densitysplit 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 overdensity 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., HarnoisDéraps et al. 2021; Zürcher et al. 2022, for examples of simulationbased 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 reducedshear correction to the aperture mass probability distribution function (PDF). Another approach to access higherorder moments is discussed in Halder et al. (2021), Halder & Barreira (2022), and Heydenreich et al. (2022), where thirdorder cosmic shear statistics were modelled directly from the bispectrum. Although, the simulationbased approach has advantages regarding the numerical incorporation of critical systematic effects such as the intrinsic alignment (IA) of galaxies (see, e.g. HarnoisDé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 nonperturbative 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 KiloDegree 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 signaltonoise 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 tophat 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 LDTlike 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 Nbody 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 subareas 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
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
where
The above relation between the two filters U and Q can be inverted,
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)
This definition is equivalent to the ‘CountsinCell’ (CiC) statistics mentioned in Gruen et al. (2016) if the filter U is defined as a tophat. 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^{2}. 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 Ufilter 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
where in this work m(θ) was the KiDS1000 mask. For example, the second part of this equation vanishes for a tophat 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 noisedominated 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 logspaced 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 fullsky simulations described in Takahashi et al. (2017) to suppress the residual differences of the modelled and measured data vector. The KiDS1000 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 Poissondistributed, we modified the distribution of the aperture number computed in the B22 model in a way that allows for superPoissonian shot noise. Inspired by F18, we achieved this by scaling the galaxy number density n_{0} with a free parameter α > 0, such that n_{0} α^{−1} can be interpreted as an effective number density of Poissonian tracers. This implies that the quantity p(N_{ap} α^{−1}δ_{m, U}) follows a lognormal distribution, instead of p(N_{ap}δ_{m, U}) as in B22. Consequently, the characteristic function Ψ, determining the parameters of the lognormal distribution, must be modified to (see Eq. (36) of B22)
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):
Since the expectation value ⟨N_{ap}δ_{m, U}⟩∝α and the variance ⟨(N_{ap} − ⟨N_{ap}δ_{m, U}⟩)^{2}δ_{m, U}⟩∝α^{2} 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 KiDS1000 data. We, therefore, do not need to include a modelling error in our uncertainty budget. Furthermore, we also tested a nonlinear galaxy bias model, where we exchanged the constant galaxy bias b with b = b_{1} + b_{2}δ_{m, U} > 0. However, b_{2} 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^{1}. KiDS was designed for weak lensing applications, producing highquality images with VSTOmegaCAM camera. Thanks to the infrared data from its overlapping partner survey VIKING (VISTA Kilodegree Infrared Galaxy survey, Edge et al. 2013), galaxies are observed in nine optical and nearinfrared 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 ‘KiDS1000’ as they cover ∼1000 deg^{2} of images; this reduces to 777.4 deg^{2} 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.
Overview of the observational KiDS1000 data.
3.1. Lens catalogues
Our primary lens catalogue is the ‘KiDSbright’ sample described in Bilicki et al. (2021, hereafter Bi21), a fluxlimited galaxy catalogue with accurate and precise photometric redshifts, z_{ph}, derived using the nine photometric bands available in the KiDS1000 data. This highly pure and complete^{2} 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. KiDSbright is limited to r < 20 mag, covers ∼1000 deg^{2} 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 machinelearning neural network algorithm implemented in the ANNz2 software (Sadeh et al. 2016) to map an input space of 9band 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^{−2}, 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 KiDSbright 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 KiDSbright sample is measured with high precision: jackknife subsampling 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 KiDSbright 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 bestestimated 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 photoz error model p_{δz}(Δz),
Following Bilicki et al. (2014), we adopted a ‘modified Lorentzian’,
which was shown in Bi21 to reproduce the photoz errors better than a Gaussian error model. To estimate the parameters a and s, Eq. (10) is fitted to the KiDSbright galaxies that also have GAMA spectroscopic redshifts, where Δz = z_{ph} − z_{sp}. The bestfit 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 KiDSbright galaxies described in Bi21. The blue line shows n_{sp}(z), i.e. the redshift distribution of KiDS galaxies for which we have spectra from the GAMA survey. The red line shows n_{ph}(z), the distribution of the full KiDSbright sample as estimated by ANNz2 with a photometric redshift cut of z_{ph} < 0.1. The black line shows our bestestimated n_{be}(z): a smoothed version of n_{ph}(z) that better accounts for photometric redshift errors. The cyan, orange, and brown lines show the third, fourth and fifth redshift bins of the KiDS1000 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 photoz 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 subsamples on the sky. For this, we split the KiDSbright × GAMA matched sample into ten subsamples along the right ascension, where each subsample has the same number of objects. We fitted a and s to each subsample with Eq. (10) and convolve the resulting p_{δz} with the full p(z_{ph}) of the KiDSbright 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 bestestimated 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 KiDSbright × 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 KiDSbright data to construct colourselected subsamples. 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 rband magnitude M_{r} and the restframe u − g colour diagram. The restframe 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 colourmagnitude diagram, which results in a line that delimits the red and blue samples that satisfy
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 bestfit 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 colourselected subsamples are particularly interesting from a galaxy formation perspective, our main cosmological results are obtained from the full lens sample, which has the highest signaltonoise.
3.2. Source catalogues
The fiducial KiDS1000 cosmic shear catalogue consists of five tomographic bins, whose redshifts are calibrated using the selforganising map (SOM) method^{3} 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 sourcelens 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 IAinfused weak lensing simulations and show that, given the KiDSbright n(z), this can be avoided by excluding the first three tomographic bins of the KiDS1000 sources from the analysis. An additional lenssource coupling complication is the socalled boost factor, which arises due to the clustering of sources with overdense and the anticorrelation of sources with underdense 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 lenssource 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 shearrelated systematic effects do not cause more than a 0.1σ shift in S_{8} ≡ σ_{8}(Ω_{m}/0.3)^{0.5} when measured by cosmic shear twopoint functions.
4. Simulated data
Besides the real KiDS1000 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 (Fullsky Lognormal Astrofields Simulation Kit) described in Xavier et al. (2016) to estimate the covariance of errors in the DSS data vector of the KiDS1000; the cosmoSLICS+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 lognormal 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 lognormal simulations produced with the publicly available FLASK tool^{4} (Xavier et al. 2016). In Hilbert et al. (2011) it is shown that lognormal random fields are a good approximation to the 1point 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 higherorder statistics in StageIII lensing surveys (see their Fig. 4)^{5}. Compared to full Nbody simulations, FLASK lognormal random fields are computationally cheap to create. The fact that FLASK outputs fullsky maps has the advantage that it can easily be masked to match the footprint of the data, making area rescaling 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 σ_{8} = 0.74, the dimensionless Hubble parameter to h = 0.7, the dark energy equationofstate parameter to w_{0} = −1 and the power spectrum powerlaw 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 matterlensing crossspectra. By assuming a flat universe throughout this paper, given the n(z) shown in Fig. 1, and using the PYCCL software package^{6} (Chisari et al. 2019) to get the 3D matter density contrast power spectrum P_{δ}(ℓ/χ, χ), we calculate the angular power spectrum by use of the Limberapproximated projection (Kaiser 1992) as
where i, j are placeholder for either the galaxy or convergence projection, such that for the lenses with redshift distribution n_{l}, while for the source with redshift distribution n_{s} we have instead
Besides these angular power spectra, FLASK needs the lognormal shift parameters κ_{0} and δ_{0}, where −κ_{0} and −δ_{0} defining the lower limits of the lognormal random variable of the convergence and matter density fields, respectively. Whereas the shift parameter κ_{0} = {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} = {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 tophat 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 KiDSbright 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^{7} and n_{eff} = 0.325 arcmin^{−2} is the mean galaxy density of the KiDSbright sample. Similarly, we populate the two source planes by Poissonsampling for each pixel a number of source galaxies n_{pix} with parameter λ = n_{eff}A_{pix}, where A_{pix} is the area of the pixel under consideration^{8}, 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)
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. cosmoSLICS+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 cosmoSLICS presented in HarnoisDéraps et al. (2019), which consists of a set of 50 simulated light cones of 100 deg^{2} each, run in a ΛCDM universe with Ω_{m} = 0.2905, Ω_{Λ} = 0.7095, Ω_{b} = 0.0473, h = 0.6898, σ_{8} = 0.836 and n_{s} = 0.969. The mocks follow the nonlinear evolution of 1536^{3} particles up to z = 0, computed by the cubep^{3}mNbody code (HarnoisDéraps et al. 2013). For Fourier modes of comoving wave number k < 2.0 h Mpc^{−1}, the cosmoSLICS threedimensional 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 kmodes (HarnoisDéraps et al. 2019), offering a sufficient resolution to model StageIII galaxy surveys. The particle data were assigned onto mass sheets at 18 redshifts and then postprocessed into 10 × 10 deg^{2} light cones. Lensing maps were produced at 18 source redshift planes for each cosmoSLICS light cone and used to interpolate lensing information onto galaxy catalogues.
Similar to B22, we construct cosmoSLICS 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 KiDSbright samples by populating the light cone mass maps with galaxies that trace the underlying dark matter field linearly, following the method presented in HarnoisDé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^{−2}.
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
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, is the matter density, D_{+}(z) is the linear growth factor, Mpc^{3}, as calibrated in Brown et al. (2002). These intrinsic ellipticity components are then combined with the cosmic shear signal by Eq. (14), resulting in an IAcontaminated 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 largescale structure significantly, as the sustained outflows of energy arising from stellar winds, supernovae and AGN reduce the clustering on intracluster 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 hydrodynamical 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 hydrodynamical 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 KiDS450 and StageIV cosmic shear data (Martinet et al. 2021), and subsequently in HarnoisDé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^{9}, more specifically by the Run2 and Run2b 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 coevolve dark matter particles of mass 6.9 × 10^{8}h^{−1} M_{⊙} and gas particles with mass 1.4 × 10^{8}h^{−1} M_{⊙}, in comoving volumes of side 352 and 640 h^{−1} 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 cosmoSLICS simulations to generate KiDS1000 sources and KiDSbright lenses for ten pseudoindependent light cones, each covering 100 deg^{2}. We repeat the same procedure on dark matteronly light cones, such that any difference is caused by the presence of baryons.
The cosmoSLICS+IA and Magneticum light cones are squareshaped, 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^{2}, 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 MonteCarlo 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 σ_{8}. We additionally vary the galaxy bias term b and the superPoisonnian shotnoise 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 simulationbased validation runs).
All varied parameters and their prior knowledge.
For the estimated covariance matrix , 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 measured from n_{r} simulated survey realisations, the posterior distribution of a model vector m that depends on n_{θ} parameters Θ is
where
The powerlaw index m is
with n_{d} being the number of data points and
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}, σ_{8}, 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^{10}. 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
where we found the maximum by running several minimisation processes. To estimate the resulting uncertainties around the MAP, we use the suggested projectedjointhighestposteriordensity (PJHPD) 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 χ^{2}/d.o.f. to quantify the goodness of fit, where the χ^{2} results from the point in the highdimensional parameter space that has the highest posterior probability. To unbias the covariance matrix , which is used to estimate the χ^{2}values, we instead of inverting with the known Hartlap factor (Hartlap et al. 2007) defined as h = (n_{r} − 1)/(n_{r} − n_{d} − 2), but rather use
To estimate the d.o.f. we measure for 1000 mock data vectors the best χ^{2} and fit a χ^{2} 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 pvalue in each case, which provides the probability of finding a χ^{2} 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 cosmoSLICS nodes that the theoretical model for the KiDSbright sample is valid and accurate for σ_{8} < 1.0, with indications that it loses accuracy for larger values of σ_{8}. This does not affect our results, given that the preferred values of σ_{8} 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 noisefree 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 KiDS1000 footprint^{11}. The measurements on the Magneticum mocks use the ten available light cones. Furthermore, we measured the shear profiles from KiDSlike 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 lenssource 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 pvalue is approximately 0.2, which shows that this is consistent with being a statistical fluctuation. We performed this validation test for the adapted and tophat filters but show the resulting shear signals, and the corresponding mean aperture number values only for the adapted filter since those of the tophat filter are very similar and would not yield more insights.
Fig. 2. Shear profiles measured with the adapted filter for the KiDSbrightlike lenses and sources from the cosmoSLICS+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 IAfree 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 ⟨N_{ap}⟩ measured in the cosmoSLICS+IA bright mocks for four quantiles, and the value of ⟨N_{ap}⟩ predicted by the model at the same cosmology. The blue error bars show the KiDS1000 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 anticorrelation between the S_{8} or σ_{8} 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_{8} and σ_{8} 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 tophat filter in Appendix C.
Fig. 4. Pipeline validation: cosmological inference with the adapted filter using the cosmoSLICS 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 DMonly and Hydro mock data for the lenses and sources, we end up with four scenarios (lenssource = DMDM, DMHydro, HydroDM and HydroHydro). Figure 5 shows the residuals between the shear profiles measured from dark matteronly mocks (DMDM) to the other three combinations for each quantile. Clearly, the deviations are well inside the expected KiDS1000 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 HarnoisDé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 matteronly or full hydrodynamical Magneticum simulations. The residuals are with respect to the dark matteronly 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 ⟨N_{ap}⟩ to compare measurements from the hydro and dark matteronly Magneticum simulations. The relative difference is always well inside the expected statistical uncertainty of KiDS1000. 
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 KiDSbright sample as our foreground lenses and the fourth and fifth KiDS1000 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 σ_{8} as well as the two astrophysical parameters b and α. We marginalised over the systematic effects parameters describing the δ⟨z⟩ and mbias 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 bestfit MAP values, listed in the first column of Table 3. In that table, the pvalues 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 χ^{2} values and for which a χ^{2}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 χ^{2} 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 χ^{2} scatter of , hence do not warrant further investigation.
Fig. 7. Shear profiles measured in the data, compared to the bestfit 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 ⟨N_{ap}⟩ to compare measurements from the KiDSbright sample and our model, evaluated at the MAP values shown in Table 3. The measured ⟨N_{ap}⟩ are all greater than the predicted ⟨N_{ap}⟩ at MAP just indicating that the measured p(N_{ap}) is broader than the predicted one. 
Marginalised MAP values and their 68% confidence intervals for the different lens n(z) of the full sample.
Using the approach described in Sect. 5 to estimate the uncertainty around the MAP, we find
which is consistent with and competitive to the KiDS1000 cosmic shear constraints from Asgari et al. (2021),
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_{8}; however, the Ω_{m}–σ_{8} degeneracy is broken, thanks to the additional information provided by the foreground data. But even for the shearonly case shown in green in Fig. 9, the DSS has better constraining power for the Ω_{m} and σ_{8} 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 cosmoSLICS for smaller Ω_{m} values. Furthermore, although the results of Fig. A.2 show that the inferred S_{8} might be smaller compared to the truth if a linear galaxy bias model is not sufficient, the consistency between the shearonly 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_{8} 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 COSEBIsDSS analysis could further improve the constraints, which we leave for future work. Lastly, as the DSS estimates of S_{8} 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 bestestimated n_{be}(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 cosmoSLICS for smaller Ω_{m} values. The shearonly 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 coloursplit 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 KiDSbright 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 largescale 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 σ_{8} 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} − σ_{8} degeneracy axis, whereas these shifts partially cancel out for S_{8}. 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 σ_{8} and smaller Ω_{m}, and viceversa. 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(z_{phot}) 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 lowredshift galaxies does not result in systematically different Ω_{m} or σ_{8}.
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 cosmologydependent 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 overconstraints, 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 nonlinear 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}),
We then multiply each element of the fiducial covariance matrix, , by the scaling factors,
and obtain a cosmologyrescaled covariance matrix.
As explained in Eifler et al. (2009), this method is only valid for a noisefree covariance matrix since it wrongly rescales the shapenoise component and possibly overestimates 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 parameterdependent 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 bestfit 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 KiDSbright 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 bestestimated 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 χ^{2}/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 KiDSbright sample shown in green, while the joint red+blue posteriors represent the results of the colourselected samples. In the latter case, the red and blue samples share the same cosmology (the dark blue contours) by construction. The resulting measured and bestfit 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. 
Marginalised MAP values and their 68% confidence intervals for the different lens samples.
First, we notice that the cosmological parameters of the full KiDSbright sample analysis and the joint red and blue analysis are consistent and within 0.48 σ in σ_{8} 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 superPoisson distribution and red galaxies a subPoisson 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 normalPoisson 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 KiloDegree survey (KiDS1000). 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 KiDSbright sample described in Bi21, while for our sources, we used the fourth and fifth tomographic redshift bins of the KiDS1000 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, σ_{8} increased, leaving S_{8} values changed by ∼0.5 σ. We assigned an extra error term to this uncertainty, resulting in for the n(z) shape after marginalising over the other systematic effects. These constraints are competitive and consistent with the KiDS1000 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 bestfit 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 KiDSbright 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 superPoisson or subPoisson 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 superPoisson distributed and red galaxies are subPoisson 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}σ_{8} 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 darkenergyequation of state w for future analysis; this can be achieved by a tomographic analysis of highprecision 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 nonlinear scales and baryonic analysis.
The KiDS data products are public and available through http://kids.strw.leidenuniv.nl/DR4
We tested that an arearescaled covariance matrix coming from over 600 fully independent Nbody simulations (see HarnoisDéraps et al. 2018, for a description of the SLICS simulation suite) results in similar constraints.
Currently available here: https://github.com/LSSTDESC/CCL
Also different values would not affect the posteriors as we discuss in Sect. 7.2.
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 KiDSinternal 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 NewtonKavliJunior Fellowship and by Churchill College Cambridge through a postdoctoral ByFellowship. 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/51) 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 ERC2019AdG 882679 and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC2094 – 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 PlanckHumboldt 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 CMSCSST2021A01 and CMSCSST2021B01, 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. ZDBSLY7013. 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
 Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526 [NASA ADS] [CrossRef] [Google Scholar]
 Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514 [NASA ADS] [CrossRef] [Google Scholar]
 Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540 [Google Scholar]
 Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asgari, M., Lin, C.A., Joachimi, B., et al. 2021, A&A, 645, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barthelemy, A., Codis, S., & Bernardeau, F. 2021, MNRAS, 503, 5204 [Google Scholar]
 Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992 [CrossRef] [Google Scholar]
 Bilicki, M., Jarrett, T. H., Peacock, J. A., Cluver, M. E., & Steward, L. 2014, ApJS, 210, 9 [Google Scholar]
 Bilicki, M., Dvornik, A., Hoekstra, H., et al. 2021, A&A, 653, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bridle, S., & King, L. 2007, New J. Phys., 9, 444 [Google Scholar]
 Brouwer, M. M., Demchenko, V., HarnoisDéraps, J., et al. 2018, MNRAS, 481, 5189 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. L., Taylor, A. N., Hambly, N. C., & Dye, S. 2002, MNRAS, 333, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Burger, P., Schneider, P., Demchenko, V., et al. 2020, A&A, 642, A161 [EDP Sciences] [Google Scholar]
 Burger, P., Friedrich, O., HarnoisDéraps, J., & Schneider, P. 2022, A&A, 661, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carron, J. 2013, A&A, 551, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Castro, T., Borgani, S., Dolag, K., et al. 2021, MNRAS, 500, 2316 [Google Scholar]
 Chisari, N., Codis, S., Laigle, C., et al. 2015, MNRAS, 454, 2736 [NASA ADS] [CrossRef] [Google Scholar]
 Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2 [Google Scholar]
 de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DES Collaboration (Abbott, T. M. C., et al.) 2022, Phys. Rev. D, 105, 023520 [CrossRef] [Google Scholar]
 Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021a, Astropart. Phys., 131, 102605 [NASA ADS] [CrossRef] [Google Scholar]
 Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021b, Astropart. Phys., 131, 102604 [NASA ADS] [CrossRef] [Google Scholar]
 Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971 [Google Scholar]
 Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32 [NASA ADS] [Google Scholar]
 Eifler, T., Schneider, P., & Hartlap, J. 2009, A&A, 502, 721 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fenech Conti, I., Herbonnet, R., Hoekstra, H., et al. 2017, MNRAS, 467, 1627 [NASA ADS] [Google Scholar]
 Friedrich, O., Gruen, D., DeRose, J., et al. 2018, Phys. Rev. D, 98, 023508 [Google Scholar]
 Friedrich, O., Uhlemann, C., VillaescusaNavarro, F., et al. 2020, MNRAS, 498, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Friedrich, O., Halder, A., Boyle, A., et al. 2022, MNRAS, 510, 5069 [NASA ADS] [CrossRef] [Google Scholar]
 Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725 [Google Scholar]
 Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Gruen, D., Friedrich, O., Amara, A., et al. 2016, MNRAS, 455, 3367 [Google Scholar]
 Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507 [Google Scholar]
 Halder, A., & Barreira, A. 2022, MNRAS, 515, 4639 [NASA ADS] [CrossRef] [Google Scholar]
 Halder, A., Friedrich, O., Seitz, S., & Varga, T. N. 2021, MNRAS, 506, 2780 [CrossRef] [Google Scholar]
 Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16 [Google Scholar]
 HarnoisDéraps, J., Pen, U.L., Iliev, I. T., et al. 2013, MNRAS, 436, 540 [Google Scholar]
 HarnoisDéraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337 [Google Scholar]
 HarnoisDéraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 HarnoisDéraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623 [CrossRef] [Google Scholar]
 HarnoisDéraps, J., Martinet, N., & Reischke, R. 2022, MNRAS, 509, 3868 [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2014, ApJ, 780, 111 [Google Scholar]
 Heydenreich, S., Brück, B., Burger, P., et al. 2022, A&A, 667, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hilbert, S., Hartlap, J., & Schneider, P. 2011, A&A, 536, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454 [Google Scholar]
 Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124 [EDP Sciences] [Google Scholar]
 Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304 [Google Scholar]
 Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338 [Google Scholar]
 Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894 [NASA ADS] [CrossRef] [Google Scholar]
 Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaiser, N. 1992, ApJ, 388, 272 [Google Scholar]
 Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kodwani, D., Alonso, D., & Ferreira, P. 2019, Open J. Astrophys., 2, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500 [Google Scholar]
 Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martinet, N., Castro, T., HarnoisDéraps, J., et al. 2021, A&A, 648, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858 [Google Scholar]
 Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution, (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Peacock, J. A., & Bilicki, M. 2018, MNRAS, 481, 1133 [Google Scholar]
 Percival, W. J., Friedrich, O., Sellentin, E., & Heavens, A. 2022, MNRAS, 510, 3207 [NASA ADS] [CrossRef] [Google Scholar]
 Pires, S., Leonard, A., & Starck, J.L. 2012, MNRAS, 423, 983 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pyne, S., & Joachimi, B. 2021, MNRAS, 503, 2300 [NASA ADS] [CrossRef] [Google Scholar]
 Ragagnin, A., Dolag, K., Biffi, V., et al. 2017, Astron. Comput., 20, 52 [Google Scholar]
 Reimberg, P., & Bernardeau, F. 2018, Phys. Rev. D, 97, 023524 [NASA ADS] [CrossRef] [Google Scholar]
 Sadeh, I., Abdalla, F. B., & Lahav, O. 2016, PASP, 128, 104502 [NASA ADS] [CrossRef] [Google Scholar]
 Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640 [Google Scholar]
 Schneider, P. 1996, MNRAS, 283, 837 [Google Scholar]
 Schneider, P. 1998, ApJ, 498, 43 [Google Scholar]
 Seitz, C., & Schneider, P. 1997, A&A, 318, 687 [NASA ADS] [Google Scholar]
 Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132 [Google Scholar]
 Smith, A., Cole, S., Baugh, C., et al. 2017, MNRAS, 470, 4646 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24 [Google Scholar]
 van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649 [Google Scholar]
 van den Busch, J. L., Wright, A. H., Hildebrandt, H., et al. 2022, A&A, 664, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, A. H., Hildebrandt, H., van den Busch, J. L., et al. 2020, A&A, 640, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Xavier, H. S., Abdalla, F. B., & Joachimi, B. 2016, MNRAS, 459, 3693 [NASA ADS] [CrossRef] [Google Scholar]
 Zürcher, D., Fluri, J., Sgier, R., et al. 2022, MNRAS, 511, 2075 [CrossRef] [Google Scholar]
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 b_{2} 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 , 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 KiDS1000 bins and lenses that mimic the KiDSbright sample. For the HOD analysis, we measure the data vector, and the covariance from 614 SLICS realisations with noisy KiDS1000 sources and GAMA lens mocks (see HarnoisDéraps et al. 2018, for a detailed description), where the covariance is scaled to approximately match the KIDS1000 footprint. It is clearly seen that the posterior for b_{2} > 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_{2} ≈ 0.
Fig. A.3. Redshift distribution of the full KiDSbright 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 bestestimated n_{be}(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 cosmoSLICS 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 χ^{2} distribution of mock data vectors around the MAP for the adapted filter with the bestestimated 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 KiDSbright samples in red and blue, with the full KiDSbright 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 KiDS1000 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(N_{ap}) is smaller than the predicted Map, resulting in the measured ⟨N_{ap}⟩ being smaller. 
Appendix C: Tophat filter analysis
Fig. C.1. MCMC results for the tophat filter using the model and simulations where we infused the pure shear signal with different amplitudes of IA. 
Marginalised MAP values and their 68% confidence intervals for the different lens n(z) of the full sample.
This section shows the corresponding plots that belong to the analysis with the tophat filter. We start by showing in Fig. C.1 the validation for the tophat with the cosmoSLICS, 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 tophat filter in Table C.1, which shows with the given reduced χ^{2} and corresponding pvalue 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 tophat filter to the adapted filter as shown in Fig. 12 and Table 4. The pvalue 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 χ^{2} 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 subsamples. 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 tophat 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 KiDSbright sample shown in green, and the joint red+blue posteriors that represent the results of the individual colourselected samples. In the latter case, by construction, the red and blue samples share the same cosmology (the dark blue contours). 
Marginalised MAP values and their 68% confidence intervals for the different lens samples.
All Tables
Marginalised MAP values and their 68% confidence intervals for the different lens n(z) of the full sample.
Marginalised MAP values and their 68% confidence intervals for the different lens samples.
Marginalised MAP values and their 68% confidence intervals for the different lens n(z) of the full sample.
Marginalised MAP values and their 68% confidence intervals for the different lens samples.
All Figures
Fig. 1. Redshift distributions, n(z), of the galaxy samples. The lens sample is obtained from the KiDSbright galaxies described in Bi21. The blue line shows n_{sp}(z), i.e. the redshift distribution of KiDS galaxies for which we have spectra from the GAMA survey. The red line shows n_{ph}(z), the distribution of the full KiDSbright sample as estimated by ANNz2 with a photometric redshift cut of z_{ph} < 0.1. The black line shows our bestestimated n_{be}(z): a smoothed version of n_{ph}(z) that better accounts for photometric redshift errors. The cyan, orange, and brown lines show the third, fourth and fifth redshift bins of the KiDS1000 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 KiDSbrightlike lenses and sources from the cosmoSLICS+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 IAfree 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 ⟨N_{ap}⟩ measured in the cosmoSLICS+IA bright mocks for four quantiles, and the value of ⟨N_{ap}⟩ predicted by the model at the same cosmology. The blue error bars show the KiDS1000 statistical uncertainty measured from FLASK. 

In the text 
Fig. 4. Pipeline validation: cosmological inference with the adapted filter using the cosmoSLICS 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 matteronly or full hydrodynamical Magneticum simulations. The residuals are with respect to the dark matteronly 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 ⟨N_{ap}⟩ to compare measurements from the hydro and dark matteronly Magneticum simulations. The relative difference is always well inside the expected statistical uncertainty of KiDS1000. 

In the text 
Fig. 7. Shear profiles measured in the data, compared to the bestfit 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 ⟨N_{ap}⟩ to compare measurements from the KiDSbright sample and our model, evaluated at the MAP values shown in Table 3. The measured ⟨N_{ap}⟩ are all greater than the predicted ⟨N_{ap}⟩ at MAP just indicating that the measured p(N_{ap}) is broader than the predicted one. 

In the text 
Fig. 9. Cosmological posteriors for the adapted filter for the bestestimated n_{be}(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 cosmoSLICS for smaller Ω_{m} values. The shearonly 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(z_{phot}) 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 bestfit parameters (red) and to those obtained by varying the covariance with the parameters (blue). 

In the text 
Fig. 12. Posterior for the full KiDSbright sample shown in green, while the joint red+blue posteriors represent the results of the colourselected samples. In the latter case, the red and blue samples share the same cosmology (the dark blue contours) by construction. The resulting measured and bestfit 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 b_{2} 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 KiDSbright 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 bestestimated n_{be}(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 cosmoSLICS 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 KiDSbright samples in red and blue, with the full KiDSbright 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 KiDS1000 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(N_{ap}) is smaller than the predicted Map, resulting in the measured ⟨N_{ap}⟩ being smaller. 

In the text 
Fig. C.1. MCMC results for the tophat 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 KiDSbright sample shown in green, and the joint red+blue posteriors that represent the results of the individual colourselected samples. In the latter case, by construction, the red and blue samples share the same cosmology (the dark blue contours). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.