Introduction
Our current understanding of the large-scale structure (LSS) of the Universe, obtained over the last three decades thanks to increasingly sophisticated astronomical observations, can be regarded as one of the great successes of modern scientific inquiry. Detailed observations of the LSS probe fundamental physics, since the structures we observe were formed by an interplay of physical effects. These primarily include the self-gravity of matter, cosmic expansion, astrophysical feedback processes, and the repulsive effect of the cosmological constant or dark energy.
When it comes to understanding the growth and evolution of LSS, the cosmological stage is set by observations of the cosmic microwave background (CMB) that reveal the structure of dark and baryonic matter at early times. Together with other probes, including big-bang nucleosynthesis (BBN) and baryon-acoustic oscillations (BAOs), the CMB has been used to establish and tightly constrain the so-called standard model of cosmology, Λ \Lambda Λ CDM. It is a strikingly simple model based on only a handful of numerical parameters ([103 ] (Qu et al. 2024; Ge et al. 2025)) and characterised, in particular, by two components that are at the origin of its name: a cosmological constant (Λ \Lambda Λ , which dominates the energy landscape at late times) and a cold dark matter (CDM) component that is dominant over baryonic matter by a factor of roughly five. Observations of the CMB constrain the parameters of this model at high redshift (z ∼ 1100 z \sim1100 z ∼ 1100 ), when the Universe was approximately 370 thousand years old. However, the model can be subsequently used to predict the properties of the LSS at late times (z ≲ 2 z \lesssim2 z ≲ 2 ; an extrapolation over more than 10 billion years), which can then be tested with low-redshift probes of cosmic structures.
Of the low-redshift probes, the relative prevalence of CDM poses a problem for those that rely on luminous tracers, such as galaxy clustering measurements. Modelling and constraining the relationship between luminous and dark matter (DM), also known as galaxy bias, constitutes one of the major systematic challenges of such measurements. In contrast, observations of the weak gravitational lensing effect of the LSS are intrinsically sensitive to all matter (dark or otherwise), eliminating the systematic uncertainty of galaxy bias. This unique feature has turned such ‘cosmic shear’ experiments into a major cosmological probe that is just entering its golden age. Indeed, it is the combination of cosmic shear with galaxy clustering – and their cross-correlation called galaxy-galaxy lensing (GGL) – that will potentially yield the most precise test of the standard model.
Over the past decade, the comparison between observations of low-z z z LSS and their CMB-based predictions has revealed a tension in the S 8 = σ 8 Ω m / 0.3 S_8 = \sigma_8\sqrt{\Omega_{\mathrm{m}}/0.3} S 8 = σ 8 Ω m /0.3 parameter. Therefore, it has been the focus of extensive study. All contemporary cosmic shear measurements to date, except one, have yielded values of S 8 S_8 S 8 that are lower than CMB predictions, indicating a universe with less structure and/or matter. The significance of these measurements varies, but there is a clear tendency towards lower values in the cosmic shear results of the past decade. Only the Deep Lens Survey (DLS; Jee et al. [64 ] ) has reported a value that is higher than that of Planck , and even then only slightly. In contrast, lower values have been reported by the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS; Heymans et al. [50 ] ), as well as from previous results from the ongoing ‘stage-III’ surveys by the Kilo-Degree Survey (KiDS; [8 ] (Li et al. 2023a)), the Dark Energy Survey (DES; [5 , 121 ] ), the Hyper-Suprime Camera survey (HSC; [27 ] (Li et al. 2023c)), and a joint analysis by Dark Energy Survey and Kilo-Degree Survey Collaboration et al. (2023). Other low-z z z LSS probes typically also share this trend but again with a varying level of significance (see Abdalla et al. [1 ] , for an overview).
Cosmic shear probes relatively low redshifts, and relatively small scales compared to other probes of LSS (see e.g. Figure 1 of Broxterman & Kuijken [21 ] , for a graphical summary). This makes cosmic shear somewhat more sensitive than other low-z z z probes to non-linear structure formation and astrophysical effects: the latter mainly involve feedback processes driven by supernovae and active galactic nuclei, collectively referred to below as baryonic feedback. It has therefore been suggested that inaccuracies in the modelling of non-linear scales might be the reason for the tension in S 8 S_8 S 8 , rather than a more fundamental cosmological cause ([143 , 142 , 4 , 104 ] ).
A definitive answer to the reality of the S 8 S_8 S 8 tension is within reach. Currently, three stage-III surveys have finished their observations and are preparing their final analyses; one of which we present here. The next generation of surveys, called stage-IV, are just beginning, represented by the ESA Euclid mission [33 ] , the Legacy Survey of Space and Time (LSST) conducted with the Vera C. Rubin Observatory (Ivezić et al. [61 ] ), the NASA Roman Space Telescope (Spergel et al. [126 ] ), and the Chinese Space Station Optical Survey (Gong et al. [45 ] ). We note that Euclid is expected to yield first cosmologically relevant results in 2026. These stage-IV surveys will soon provide datasets capable of measuring the growth of structure at low z z z with unprecedented precision.
New facilities and new instruments invariably bring along new challenges owing to, for instance, the novelty of their hardware. Comparatively, the final (usually called ‘legacy’) releases of a survey typically represent the culmination of years or decades of analysis, understanding, and refinement regarding their survey and the associated data. As such, the final releases of stage-III surveys will provide a crucial insight and reference for subsequent lensing studies of LSS with stage-IV surveys.
Here, we present the cosmic shear analysis of the completed Kilo-Degree Survey (KiDS) based on the fifth and final data release (KiDS-DR5; Wright et al. [140 ] ), dubbed KiDS-Legacy. With extensive wavelength coverage enabled by its unique pairing with its near-infrared (NIR) sister survey VIKING (VISTA Kilo-Degree Infrared Galaxy Survey) and with its superb image quality in the r r r -band, KiDS is the stage-III survey that is arguably most suited to simultaneously leverage high-redshift galaxies and accurately calibrate their redshift distributions. Over the past decade, KiDS has been pushing to increasingly high redshifts, while continuously improving the redshift calibration, shear measurements, shear calibration through image simulations, and methodology to accurately extract cosmological constraints from the imaging data.
This paper represents the culmination of these efforts, providing the most robust and precise analysis of KiDS to date. Our dataset is the largest analysed by KiDS thus far (Sect. 3.1), utilising the most extensive spectroscopic redshift compilation prepared for any cosmic shear analysis (Sect. 3.2), and has been calibrated using the most sophisticated simulations produced for this purpose (Sect. 3.4). Improvements in methodology (Sect. 2) have allowed us to access the highest redshifts (with the most tomographic bins) studied for cosmic shear to date (Sect. 3.4). An efficient workflow for this study was enabled by the new CosmoPipe analysis framework, which we have also made publicly available. We present null tests of our dataset in Sect. 4. The cosmological results are presented in Sect. 5, discussed in Sect. 6, and summarised in Sect. 7. Additionally, we include supplementary information regarding our power spectrum emulation (Appendix A), on the definition of parameters for our fiducial intrinsic alignment model (Appendix B), construction of band power spectra (Appendix C), covariance validation (Appendix D), B B B -mode exploration (Appendix E), correlation function cosmological analyses (Appendix F), detailed posterior analyses (Appendix G), tables of cosmological constraints (Appendix H), and details of changes made to the analysis and manuscript post-unblinding (Appendix I). Finally, this paper has been released alongside two companion papers: the work of Stölzner et al. [128 ] presents an extensive analysis of internal and external consistency tests of the KiDS-Legacy dataset (including joint constraints with other datasets) and the work of Wright et al. [138 ] describes the various redshift distribution and bias estimation developments made for KiDS-Legacy.
Modelling and theoretical framework
The main cosmic shear observable is the ellipticity correlation between pairs of galaxies. The signal depends on the two-point statistics of the shear and galaxy intrinsic alignment fields, which both trace the matter distribution and therefore depend on cosmological information. For our predictions and sampling, we used CosmoSIS, a modular cosmology pipeline described in Zuntz et al. [144 ] . We included a number of new and/or modified modules in our analysis, including a CosmoPower emulator for non-linear power spectra (Appendix A) and new intrinsic alignment models (Sect. 2.4). In this section, we start with the derivation of cosmic shear power spectra from matter power spectra and then detail our choices of intrinsic alignment models. Next, we introduce our observable two-point statistics, briefly discussing their measurement and prediction methodology. We end this section with a brief review of our analytical covariance modelling.
Shear power spectrum
Our two-point statistic predictions are calculated from shear (cross-)power spectra of galaxy samples i i i and j j j , C G G ( i j ) ( ℓ ) C_{\mathrm{GG}}^{(ij)}(\ell) C GG ( ij ) ( ℓ ) , where G stands for ‘gravitational’, which distinguishes these power spectra from intrinsic alignment contributions to the two-point signal. Shear power spectra are given by Limber-approximated projections (Kaiser [70 ] ; LoVerde & Afshordi [89 ] ; Kilbinger et al. [72 ] ) taking the form
C G G ( i j ) ( ℓ ) = ∫ 0 χ h o r d χ W G ( i ) ( χ ) W G ( j ) ( χ ) f K 2 ( χ ) P m , n l ( ℓ + 1 / 2 f K ( χ ) , z ( χ ) ) , (1) C_{\mathrm{GG}}^{(ij)}(\ell) = \int_{0}^{\chi_{\mathrm{hor}}} \mathrm{d}\chi\frac{W_{\mathrm{G}}^{(i)}(\chi)W_{\mathrm{G}}^{(j)}(\chi)}{f_{K}^{2}(\chi)} P_{\mathrm{m,nl}}\left(\frac{\ell+1/2}{f_{K}(\chi)},z(\chi)\right), \tag*{(1)} C GG ( ij ) ( ℓ ) = ∫ 0 χ hor d χ f K 2 ( χ ) W G ( i ) ( χ ) W G ( j ) ( χ ) P m , nl ( f K ( χ ) ℓ + 1/2 , z ( χ ) ) , (1)
where P m , n l P_{\mathrm{m,nl}} P m , nl is the non-linear matter power spectrum and f K ( χ ) f_K(\chi) f K ( χ ) is the comoving angular diameter distance. Although the integral in Eq. (1) technically runs over the entire line of sight to the horizon, χ h o r \chi_{\mathrm{hor}} χ hor , in practice, the weak lensing kernel, W G ( i ) ( χ ) W_{\mathrm{G}}^{(i)}(\chi) W G ( i ) ( χ ) , limits the integral to the comoving distance of the furthest galaxy in our sample, expressed as
W G ( i ) ( χ ) = 3 H 0 2 Ω m 2 c 2 f K ( χ ) a ( χ ) ∫ χ χ h o r d χ ′ n S ( i ) ( χ ′ ) f K ( χ ′ − χ ) f K ( χ ′ ) , (2) W_{\mathrm{G}}^{(i)}(\chi) = \frac{3H_{0}^{2}\Omega_{\mathrm{m}}}{2c^{2}}\frac{f_{K}(\chi)}{a(\chi)}\int_{\chi}^{\chi_{\mathrm{hor}}} \mathrm{d}\chi' \, n_{\mathrm{S}}^{(i)}(\chi') \frac{f_{K}(\chi'-\chi)}{f_{K}(\chi')}, \tag*{(2)} W G ( i ) ( χ ) = 2 c 2 3 H 0 2 Ω m a ( χ ) f K ( χ ) ∫ χ χ hor d χ ′ n S ( i ) ( χ ′ ) f K ( χ ′ ) f K ( χ ′ − χ ) , (2)
where n S ( i ) ( χ ) n_{\mathrm{S}}^{(i)}(\chi) n S ( i ) ( χ ) is the distribution of comoving distances of galaxies in source sample i i i ; in practice, we express this in terms of redshift, namely, N ( z ) d z = n S ( χ ) d χ N(z)\,\mathrm{d}z=n_{\mathrm{S}}(\chi)\,\mathrm{d}\chi N ( z ) d z = n S ( χ ) d χ . Here, Ω m \Omega_{\mathrm{m}} Ω m is the total matter density parameter today, H 0 H_{0} H 0 is the Hubble constant, c c c is the speed of light, and a ( χ ) a(\chi) a ( χ ) is the scale factor normalised to unity today.
Matter power spectrum
In our cosmic shear analysis, we remained mostly sensitive to Fourier scales 0.1 h / M p c ≲ k ≲ 1 h / M p c 0.1\,h/\mathrm{Mpc}\lesssim k\lesssim1\,h/\mathrm{Mpc} 0.1 h / Mpc ≲ k ≲ 1 h / Mpc , with a weighting of scales that depends on the two-point statistic (see Sect. 2.5). These scales are well within the (quasi-) non-linear regime and therefore require an accurate non-linear power spectrum model for their predictions. Currently, there are multiple approaches that can produce non-linear matter power spectra, with varying levels of accuracy, for this range of scales. Here, we followed the most popular approach for cosmic shear analyses, using a method inspired by the halo model and calibrated against simulations. Specifically, we used the Boltzmann solver CAMB (Lewis et al. [85 ] ; Howlett et al. [60 ] ) to estimate the linear matter power spectrum and the augmented halo model ‘HM-Code2020’ (Mead et al. [95 ] ) to derive its non-linear evolution. HMCode2020 has a similar level of accuracy as popular emulators trained directly on N-body simulations (2.5%; see e.g. Knabenhans et al. [73 ] ; Aricò et al. [6 ] ) and is twice as accurate as HALOFIT from Takahashi et al. [129 ] , which has been used in a number of previous cosmic shear analyses (see e.g. Amon et al. [5 ] ; Hikage et al. [53 ] ). Dark Energy Survey and Kilo-Degree Survey Collaboration et al. (2023) presented a comparison of HMCode2020 and Halofit with EuclidEmulatorv2, showing that HMCode2020 produces unbiased results when analysing mock data created with EuclidEmulatorv2; meanwhile Halofit tends to overestimate power for lower values of S 8 S_8 S 8 , resulting in an underestimated S 8 S_8 S 8 value. Similarly, an earlier version of HMCode, HMCode2015 [95 ] , which was used in previous KiDS analyses, is ∼ 30 % \sim30\% ∼ 30% less accurate than HMCode2020. This can subsequently produce small differences between cosmological constraints estimated using the two codes: Tröster et al. [132 ] demonstrated that the most important difference between these versions of HMCode is in their treatment of baryon feedback, which results in a slightly higher S 8 S_8 S 8 using HMCode2020. Lastly, we note that the parameter values at which HMCode was fitted to simulations are somewhat more constrained than our chosen priors for h 2 Ω m h^2\Omega_{\mathrm{m}} h 2 Ω m (see Mead et al. 2021, and Table 3 ). However, given that the fitting functions of HMCode are based on physically motivated parameters, we did not anticipate any significant issues just beyond this range. Furthermore, we note that our analysis remained unaffected by this choice.
Table 3. Fiducial prior ranges used in the analysis of KiDS-Legacy.
Baryonic feedback
The primary features of the matter power spectrum can be determined by assuming all matter is in the form of (pressureless) CDM, with modifications due to processes involving baryonic matter. At large scales (0.05 h / M p c ≲ k ≲ 0.5 h / M p c 0.05\ h/\mathrm{Mpc} \lesssim k \lesssim0.5\ h/\mathrm{Mpc} 0.05 h / Mpc ≲ k ≲ 0.5 h / Mpc ), BAOs produce oscillations in the power spectrum. Over a wide range of intermediate scales (0.1 h / M p c ≲ k ≲ 10 h / M p c 0.1\ h/\mathrm{Mpc} \lesssim k \lesssim10\ h/\mathrm{Mpc} 0.1 h / Mpc ≲ k ≲ 10 h / Mpc ), power is suppressed through feedback processes from active galactic nuclei (AGN), reaching a maximum suppression at around k ∼ 5 h / M p c k \sim5\ h/\mathrm{Mpc} k ∼ 5 h / Mpc . At even smaller scales (k ≳ 10 h / M p c k \gtrsim10\ h/\mathrm{Mpc} k ≳ 10 h / Mpc ), star formation processes lead to a sharp increase in power (see e.g. Semboloni et al. 2011).
To model baryonic feedback processes, HMCode2020 performs a modification to the concentration-mass relation of the halo profiles as a function of their gas fraction; this, in turn, depends on the strength of the AGN feedback (parametrised by log ( T A G N ) ∈ [ 7.1 , 8.3 ] \log(T_{\mathrm{AGN}}) \in[7.1, 8.3] log ( T AGN ) ∈ [ 7.1 , 8.3 ] ). This feedback model is then calibrated against a suite of hydrodynamical simulations from BA-HAMAS [93 ] , which had been generated using a wide range of effective feedback amplitudes and cosmological models. Recently, Schaye et al. (2023, FLAMINGO) presented a new suite of hydrodynamical simulations, including emulators for generating power spectra [112 ] , as a mechanism for modelling baryonic effects. Moving beyond hydrodynamical simulations, Salcido et al. (2023, SP(k)) presented alternative phenomenological models for the same purpose, while Schneider & Teyssier [115 ] ; Aricò et al. [7 ] presented perturbative models that act directly on the outputs of dark-matter only simulations and can be empirically calibrated [46 , 114 ] . We did not utilise these models for our cosmological analyses here, as these models were published over the course of the analysis of KiDS-Legacyor were not (at time of writing) implemented within our Cosmosis sampling framework. However, we note that a reanalysis of our data with these models is a natural and straightforward project that will be undertaken in a future publication.
Intrinsic alignments
The intrinsic alignment (IA) of galaxies is correlated with the observed ellipticities of galaxy images, making them indistinguishable from the cosmic shear effect we wish to probe for cosmological parameter estimation. In the limit of small shears and IA, the two-point statistics of galaxy ellipticities are given by
C E E ( i j ) ( ℓ ) = C G G ( i j ) ( ℓ ) + C G I ( i j ) ( ℓ ) + C I G ( i j ) ( ℓ ) + C I I ( i j ) ( ℓ ) , (3) C_{\mathrm{EE}}^{(ij)}(\ell) = C_{\mathrm{GG}}^{(ij)}(\ell) + C_{\mathrm{GI}}^{(ij)}(\ell) + C_{\mathrm{IG}}^{(ij)}(\ell) + C_{\mathrm{II}}^{(ij)}(\ell), \tag*{(3)} C EE ( ij ) ( ℓ ) = C GG ( ij ) ( ℓ ) + C GI ( ij ) ( ℓ ) + C IG ( ij ) ( ℓ ) + C II ( ij ) ( ℓ ) , (3)
where C G G ( i j ) C_{\mathrm{GG}}^{(ij)} C GG ( ij ) is the cosmic shear signal given by Eq. (1) . Analogous Limber equations can be derived for IA contributions, which are typically grouped into two categories. First, intrinsic-intrinsic (II) correlations represent the correlated alignment of galaxies that are tidally aligned by the same ambient LSS:
C I I ( i j ) ( ℓ ) = ∫ 0 χ h o r d χ n S ( i ) ( χ ) n S ( j ) ( χ ) f K 2 ( χ ) P I I , X ( ℓ + 1 / 2 f K ( χ ) , z ( χ ) ) , (4) C_{\mathrm{II}}^{(ij)}(\ell) = \int_0^{\chi_{\mathrm{hor}}} \mathrm{d}\chi\frac{n_{\mathrm{S}}^{(i)}(\chi)n_{\mathrm{S}}^{(j)}(\chi)}{f_{\mathrm{K}}^2(\chi)}P_{\mathrm{II},\mathrm{X}}\left(\frac{\ell+1/2}{f_{\mathrm{K}}(\chi)},z(\chi)\right), \tag*{(4)} C II ( ij ) ( ℓ ) = ∫ 0 χ hor d χ f K 2 ( χ ) n S ( i ) ( χ ) n S ( j ) ( χ ) P II , X ( f K ( χ ) ℓ + 1/2 , z ( χ ) ) , (4)
Second, shear-intrinsic (GI) correlations are driven by the mutual gravitational effects of LSS on (i) galaxies within the structures, and (ii) the lensed images of background sources that these structures produce [57 ] :
C G I ( i j ) ( ℓ ) = ∫ 0 χ h o r d χ W G ( i ) ( χ ) n S ( j ) ( χ ) f K 2 ( χ ) P δ I , X ( ℓ + 1 / 2 f K ( χ ) , z ( χ ) ) , (5) C_{\mathrm{GI}}^{(ij)}(\ell) = \int_0^{\chi_{\mathrm{hor}}} \mathrm{d}\chi\frac{W_{\mathrm{G}}^{(i)}(\chi)n_{\mathrm{S}}^{(j)}(\chi)}{f_{\mathrm{K}}^2(\chi)}P_{\mathrm{\delta I},\mathrm{X}}\left(\frac{\ell+1/2}{f_{\mathrm{K}}(\chi)},z(\chi)\right), \tag*{(5)} C GI ( ij ) ( ℓ ) = ∫ 0 χ hor d χ f K 2 ( χ ) W G ( i ) ( χ ) n S ( j ) ( χ ) P δ I , X ( f K ( χ ) ℓ + 1/2 , z ( χ ) ) , (5)
with C I G ( i j ) = C G I ( j i ) C_{\mathrm{IG}}^{(ij)}=C_{\mathrm{GI}}^{(ji)} C IG ( ij ) = C GI ( j i ) . In Eqs. (4) and (5) the power spectra P I I , X P_{\mathrm{II},\mathrm{X}} P II , X and P δ I , X P_{\mathrm{\delta I},\mathrm{X}} P δ I , X are given by the chosen IA model. We note that the IA contributions feature a more compact kernel along the line of sight than the lensing kernel in Eq. (1) , and hence the Limber approximation is poorer for these signals. However, this effect is not significant for our analysis, as the overall IA contributions to the total signal is small, making up only around 10 percent of the signal in our analysis. Additionally, we are only probing scales ℓ ≳ 50 \ell\gtrsim50 ℓ ≳ 50 , where the Limber approximation is sufficient even for more compact kernels [84 ] .
Since IA are caused by the local matter environment, they have a different redshift scaling than lensing, which however is insufficient to isolate the cosmological signal without substantial information loss [67 ] . Therefore, cosmic shear analyses opt for physically motivated IA models that are inferred jointly with the cosmology (see Lamman et al. 2024, for a review of the field). Guided by this, we implement four IA models in KiDS-Legacy: a widely used baseline with a single overall amplitude parameter (NLA, Sect. 2.4.1), two minimal extensions of NLA with either an additional redshift (NLA-z, Sect. 2.4.2) or scale dependence (via our restricted TATT model, NLA-k; Sect. 2.4.3), and a version of NLA that explicitly incorporates the well-established trends of IA with galaxy type and host halo mass (NLA-M, Sect. 2.4.4). The latter model is the fiducial model chosen for our cosmological analyses.
The non-linear linear alignment (NLA) model
The so-called NLA model [19 ] assumes that the intrinsic ellipticity of a galaxy depends linearly on its surrounding tidal quadrupole [24 ] . Moreover, it is commonly assumed that it is the tidal field at some early time during the formation of the galaxy that sets the amplitude of alignment [57 ] . This leads to the following expressions for the IA power spectra,
P δ I , N L A ( k , z ) = − A I A C 1 ρ c r Ω m D ( z ) P m , n l ( k , z ) , P I I , N L A ( k , z ) = A I A 2 ( C 1 ρ c r Ω m D ( z ) ) 2 P m , n l ( k , z ) , (6) \begin{aligned}
P_{\delta I,\mathrm{NLA}}(k,z) &= -A_{\mathrm{IA}}\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}P_{\mathrm{m,nl}}(k,z),\\
P_{II,\mathrm{NLA}}(k,z) &= A_{\mathrm{IA}}^2\left(\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\right)^2P_{\mathrm{m,nl}}(k,z), \tag*{(6)}
\end{aligned} P δ I , NLA ( k , z ) P I I , NLA ( k , z ) = − A IA D ( z ) C 1 ρ cr Ω m P m , nl ( k , z ) , = A IA 2 ( D ( z ) C 1 ρ cr Ω m ) 2 P m , nl ( k , z ) , (6)
where C 1 ρ c r ≈ 0.0134 C_1\rho_{\mathrm{cr}} \approx0.0134 C 1 ρ cr ≈ 0.0134 is a constant and D ( z ) D(z) D ( z ) is the linear growth factor normalised to unity at z = 0 z = 0 z = 0 . The model is empirically extended to non-linear scales by employing the full, non-linear matter power spectrum P m , n l P_{\mathrm{m,nl}} P m , nl , encapsulating a mixture of non-linear gravitational effects as well as alignment processes on the scales of dark matter halos (cf. Schneider & Bridle [116 ] ). The model has a single free parameter, namely, the global, dimensionless amplitude A I A A_{\mathrm{IA}} A IA , which necessitates the assumption that there are no coherent trends in how galaxy ellipticities respond to their local tidal field within the survey’s galaxy population. Since A I A A_{\mathrm{IA}} A IA is well-constrained by tomographic cosmic shear data, we chose a wide top-hat prior for A I A A_{\mathrm{IA}} A IA to ensure its posterior is data-driven. All previous KiDS cosmic shear analyses have used the NLA model as their default and, to date, there has been no evidence that this simple prescription is insufficient in representing current-generation datasets (e.g. Dark Energy Survey and Kilo-Degree Survey Collaboration et al. [28 ] ).
Additional redshift dependence (NLA-z)
There is currently no observational evidence of redshift evolution in the amplitude of IA (Joachimi et al. [66 ] ; Fortuna et al. [38 ] ), but the statistical uncertainties on measurements underpinning this statement are still large, especially when extrapolated to the redshift range of cosmic shear samples. Moreover, A I A A_{\mathrm{IA}} A IA could also acquire an effective redshift scaling; for example, due to a combination of a luminosity dependence and Malmquist bias. Hence, we tested whether the data prefer A I A A_{\mathrm{IA}} A IA to evolve over the tomographic bins. Since individual IA amplitudes per tomographic bin are poorly constrained (see the detailed analysis in Samuroff et al. [110 ] ), we chose an extended model that is linear in terms of the scale factor,
P δ I , z d e p ( i ) ( k , z ) = − ( A I A + B I A [ ⟨ a ⟩ ( i ) a p i v − 1 ] ) C 1 ρ c r Ω m D ( z ) P m , n l ( k , z ) , P I I , z d e p ( i j ) ( k , z ) = ( A I A + B I A [ ⟨ a ⟩ ( i ) a p i v − 1 ] ) ( A I A + B I A [ ⟨ a ⟩ ( j ) a p i v − 1 ] ) ( C 1 ρ c r Ω m D ( z ) ) 2 P m , n l ( k , z ) , (7) \begin{aligned}
P_{\delta I,\mathrm{zdep}}^{(i)}(k,z) &= -\left(A_{\mathrm{IA}}+B_{\mathrm{IA}}\left[\frac{\langle a\rangle^{(i)}}{a_{\mathrm{piv}}}-1\right]\right)\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}P_{\mathrm{m,nl}}(k,z),\\
P_{II,\mathrm{zdep}}^{(ij)}(k,z) &= \left(A_{\mathrm{IA}}+B_{\mathrm{IA}}\left[\frac{\langle a\rangle^{(i)}}{a_{\mathrm{piv}}}-1\right]\right)\left(A_{\mathrm{IA}}+B_{\mathrm{IA}}\left[\frac{\langle a\rangle^{(j)}}{a_{\mathrm{piv}}}-1\right]\right)\left(\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\right)^2P_{\mathrm{m,nl}}(k,z), \tag*{(7)}
\end{aligned} P δ I , zdep ( i ) ( k , z ) P I I , zdep ( ij ) ( k , z ) = − ( A IA + B IA [ a piv ⟨ a ⟩ ( i ) − 1 ] ) D ( z ) C 1 ρ cr Ω m P m , nl ( k , z ) , = ( A IA + B IA [ a piv ⟨ a ⟩ ( i ) − 1 ] ) ( A IA + B IA [ a piv ⟨ a ⟩ ( j ) − 1 ] ) ( D ( z ) C 1 ρ cr Ω m ) 2 P m , nl ( k , z ) , (7)
where ⟨ a ⟩ ( i ) \langle a\rangle^{(i)} ⟨ a ⟩ ( i ) is the average scale factor in tomographic bin i i i , evaluated using the estimated N ( z ) N(z) N ( z ) . The parameter A I A A_{\mathrm{IA}} A IA thus is the IA amplitude at a p i v ≈ 0.769 a_{\mathrm{piv}} \approx0.769 a piv ≈ 0.769 , which corresponds to the commonly chosen pivot of z = 0.3 z = 0.3 z = 0.3 for IA studies (Joachimi et al. [66 ] ). We set the prior range of the second parameter, B I A B_{\mathrm{IA}} B IA , by fitting our linear amplitude dependence to the compilation of KiDS and SDSS measurements studied in Fortuna et al. [38 ] , see their Figure 10. Increasing the resulting constraint by 50%, to account for the limited coverage of these calibration samples and for potential additional implicit redshift dependencies, we obtain a Gaussian prior on B I A B_{\mathrm{IA}} B IA with mean − 3.7 -3.7 − 3.7 and standard deviation 4.3. For direct comparisons between amplitudes estimated between this model and others, we define
A N L A z , t o t a l ( i ) = A ^ I A + B ^ I A ( ⟨ a ⟩ ( i ) a p i v − 1 ) , (8) A_{\mathrm{NLAz},\mathrm{total}}^{(i)}=\hat{A}_{\mathrm{IA}}+\hat{B}_{\mathrm{IA}}\left(\frac{\langle a\rangle^{(i)}}{a_{\mathrm{piv}}}-1\right), \tag*{(8)} A NLAz , total ( i ) = A ^ IA + B ^ IA ( a piv ⟨ a ⟩ ( i ) − 1 ) , (8)
where ^ \hat{\ } ^ denotes the best-fitted (i.e. maximum a posteriori) value of each parameter.
Restricted TATT model (NLA-k)
By design, the angular scale dependence of the NLA model is identical to that of the matter power spectrum projected over the redshift distribution of the sample. A natural way to extend this is to consider next-to-linear alignment effects, specifically: (i) a density weighting (Hirata & Seljak [57 ] ), accounting for the fact that alignments are only measured at local overdensities of the matter distribution where galaxies reside; and (ii) quadratic alignments as expected for tidal torquing effects (Crittenden et al. [25 ] ; Catelan et al. [24 ] ).
These contributions are subsumed into the tidal alignment, tidal torquing (TATT) model (Blazek et al. [18 ] ), each with a free amplitude parameter. As can be seen from Figure 1 of Blazek et al. [18 ] , the density-weighting and quadratic terms have very similar dependence on physical scale. Owing to the scale-mixing of angular statistics, cosmic shear measurements poorly distinguish small and localised scale differences, which leads to pronounced degeneracies among the four TATT parameters (e.g. Samuroff et al. [110 ] ). To avoid these, we implement the full TATT model within the analysis pipeline, but fix all but one of the free parameters that govern the non-linear terms. We treat this model as a more general modification to NLA, to capture potential deviations in scale from the model. We chose to keep the density weighting amplitude term free, thereby adding a single additional free parameter, but note that this choice is fairly inconsequential: given the aforementioned similarity between the scale dependence of the terms, implementing only the density weighting or only the quadratic term would produce a similar end product. The NLA-k k k model is then expressed as
Figure 1. Integrands of the transformation between the angular power spectrum and E n / B n E_n/B_n E n / B n (Eq. (17) , panel ‘a’), C E , B C_{\mathrm{E,B}} C E , B (Eq. C.5, panel ‘b’), and ξ ± \xi_{\pm} ξ ± (Eq. (15) , panels ‘c’ and ‘d’). All solid integrands are normalised by their maximum absolute value. For both COSEBIs and 2PCFs, we show the limiting ranges of the statistics (n ∈ { 1 , 6 } n \in\{1,6\} n ∈ { 1 , 6 } for E n E_n E n and θ ∈ { 2.0 , 300.0 } \theta\in\{2.0,300.0\} θ ∈ { 2.0 , 300.0 } arcmin for ξ ± \xi_{\pm} ξ ± ). For C E , B C_{\mathrm{E,B}} C E , B we show all eight bins used in KiDS-Legacy, defined between ℓ = 100 \ell=100 ℓ = 100 and ℓ = 1500 \ell=1500 ℓ = 1500 . Note that, unlike previous KiDS analyses, our new treatment of T ( θ ) T(\theta) T ( θ ) means that the C E , B C_{\mathrm{E,B}} C E , B initially go to zero within their desired ℓ \ell ℓ ranges.
P δ I , s c a l e d e p ( k , z ) = P δ I , N L A ( k , z ) + A I A b s r c C 1 ρ c r Ω m D ( z ) [ A 0 I E ( k , z ) + C 0 I E ( k , z ) ] , P_{\delta I,\mathrm{scaledep}}(k,z)=P_{\delta I,\mathrm{NLA}}(k,z)+A_{\mathrm{IA}}b_{\mathrm{src}}\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\left[A_{0\mathrm{IE}}(k,z)+C_{0\mathrm{IE}}(k,z)\right], P δ I , scaledep ( k , z ) = P δ I , NLA ( k , z ) + A IA b src D ( z ) C 1 ρ cr Ω m [ A 0 IE ( k , z ) + C 0 IE ( k , z ) ] ,
P I I , s c a l e d e p ( k , z ) = P I I , N L A ( k , z ) + A I A 2 b s r c 2 ( C 1 ρ c r Ω m D ( z ) ) 2 A 0 E 0 E ( k , z ) + 2 A I A 2 b s r c ( C 1 ρ c r Ω m D ( z ) ) 2 [ A 0 I E ( k , z ) + C 0 I E ( k , z ) ] , (9) \begin{aligned}
P_{II,\mathrm{scaledep}}(k,z)&=P_{II,\mathrm{NLA}}(k,z)\\
&\quad+A_{\mathrm{IA}}^2b_{\mathrm{src}}^2\left(\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\right)^2A_{0\mathrm{E}0\mathrm{E}}(k,z)\\
&\quad+2A_{\mathrm{IA}}^2b_{\mathrm{src}}\left(\frac{C_1\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\right)^2\left[A_{0\mathrm{IE}}(k,z)+C_{0\mathrm{IE}}(k,z)\right], \tag*{(9)}
\end{aligned} P I I , scaledep ( k , z ) = P I I , NLA ( k , z ) + A IA 2 b src 2 ( D ( z ) C 1 ρ cr Ω m ) 2 A 0 E 0 E ( k , z ) + 2 A IA 2 b src ( D ( z ) C 1 ρ cr Ω m ) 2 [ A 0 IE ( k , z ) + C 0 IE ( k , z ) ] , (9)
where the expressions for A 0 I E A_{0\mathrm{IE}} A 0 IE , C 0 I E C_{0\mathrm{IE}} C 0 IE , and A 0 E 0 E A_{0\mathrm{E}0\mathrm{E}} A 0 E 0 E are given in Blazek et al. [18 ] . We consider these terms in the Limber approximation where only transverse modes contribute, and assume that their redshift evolution is inherited from the linear matter power spectrum. In the original density-weighting term, the free parameter b s r c b_{\mathrm{src}} b src can be interpreted as an effective galaxy bias of the weak lensing source sample, albeit on fairly non-linear scales. For our model, we chose a top-hat prior in the range [ − 0.5 , 1.5 ] [-0.5, 1.5] [ − 0.5 , 1.5 ] to ensure it encompasses NLA (b s r c = 0 b_{\mathrm{src}} = 0 b src = 0 ) and reasonable values expected for the (unknown) source galaxy bias (b s r c ∼ 1 b_{\mathrm{src}} \sim1 b src ∼ 1 ). For comparisons between the alignment amplitude estimated with this model and others, we report the alignment amplitude of the NLA component of the model from Eq. (9) , namely,
A I A , t o t a l s c a l e d e e p = A ^ I A . (10) A_{\mathrm{IA,total}}^{\mathrm{scaledeep}}=\widehat{A}_{\mathrm{IA}}. \tag*{(10)} A IA , total scaledeep = A IA . (10)
Galaxy-type and mass-dependent model (NLA-M)
Two trends in IA strength are firmly established observationally: (i) a dichotomy between late- and early-type galaxies where only the latter have yielded significant detections on cosmological scales ([111 , 68 , 50 , 43 ] ; see also [91 , 131 , 94 ] for additional examples of null or marginal detections of late-type alignments); and (ii) a scaling of early-type alignment amplitude with host halo mass that is well described by a power-law ([66 , 36 ] ). Both have solid theoretical underpinning. First, rotationally supported late-type galaxies are expected to align through tidal torquing, which depends quadratically on the tidal quadrupole, whereas spheroid shapes should respond linearly and thus be well described by NLA-like models. Second, Piras et al. [102 ] motivated the halo mass dependence through an analytical argument and also reproduced it for simulated dark matter haloes.
We incorporated these IA trends explicitly into the amplitude of the NLA model. We assumed that, in each tomographic bin, only the fraction f r f_{\mathrm{r}} f r of ‘red’ (i.e. early-type) galaxies intrinsically align, whereas the remaining galaxies have zero alignment. The latter assumption is consistent with observations, but within a large statistical uncertainty and involving a substantial extrapolation in sample properties, primarily due to the imperfect relation between galaxy colour and rotational/pressure support (which is the property most-likely to cleanly stratify intrinsic alignment amplitudes). Moreover, we assumed that the red galaxy population IA amplitude has a power-law dependence on the average halo mass ⟨ M h ⟩ \langle M_{\mathrm{h}}\rangle ⟨ M h ⟩ within a tomographic bin. This leads to
P δ I , m a s s d e p ( i ) ( k , z ) = − A I A f r ( i ) ( ⟨ M h ⟩ ( i ) M h , p i v ) β C 1 ρ c r Ω m D ( z ) P m , n l ( k , z ) , (11) P_{\delta\mathrm{I},\mathrm{massdep}}^{(i)}(k,z)=-A_{\mathrm{IA}}f_{\mathrm{r}}^{(i)}\left(\frac{\langle M_{\mathrm{h}}\rangle^{(i)}}{M_{\mathrm{h,piv}}}\right)^{\beta}\frac{C_{1}\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}P_{\mathrm{m,nl}}(k,z), \tag*{(11)} P δ I , massdep ( i ) ( k , z ) = − A IA f r ( i ) ( M h , piv ⟨ M h ⟩ ( i ) ) β D ( z ) C 1 ρ cr Ω m P m , nl ( k , z ) , (11)
P I I , m a s s d e p ( i j ) ( k , z ) = A I A 2 f r ( i ) f r ( j ) ( ⟨ M h ⟩ ( i ) M h , p i v ) β ( ⟨ M h ⟩ ( j ) M h , p i v ) β ( C 1 ρ c r Ω m D ( z ) ) 2 P m , n l ( k , z ) , P_{\mathrm{II},\mathrm{massdep}}^{(ij)}(k,z)=A_{\mathrm{IA}}^{2}f_{\mathrm{r}}^{(i)}f_{\mathrm{r}}^{(j)}\left(\frac{\langle M_{\mathrm{h}}\rangle^{(i)}}{M_{\mathrm{h,piv}}}\right)^{\beta}\left(\frac{\langle M_{\mathrm{h}}\rangle^{(j)}}{M_{\mathrm{h,piv}}}\right)^{\beta}\left(\frac{C_{1}\rho_{\mathrm{cr}}\Omega_{\mathrm{m}}}{D(z)}\right)^{2}P_{\mathrm{m,nl}}(k,z), P II , massdep ( ij ) ( k , z ) = A IA 2 f r ( i ) f r ( j ) ( M h , piv ⟨ M h ⟩ ( i ) ) β ( M h , piv ⟨ M h ⟩ ( j ) ) β ( D ( z ) C 1 ρ cr Ω m ) 2 P m , nl ( k , z ) ,
where we adopt the pivot halo mass M h , p i v = 10 13.5 h − 1 M ⊙ M_{\mathrm{h,piv}}=10^{13.5}\,h^{-1}M_{\odot} M h , piv = 1 0 13.5 h − 1 M ⊙ from Fortuna et al. [36 ] . The same authors also fit the power-law dependence to a compilation of early-type IA measurements (see their Figure 4), and we adopt their joint posterior on A I A A_{\mathrm{IA}} A IA and β \beta β as our prior for cosmic shear inference. Their posterior is well approximated by a bivariate Gaussian with A I A = 5.74 ± 0.29 A_{\mathrm{IA}}=5.74\pm0.29 A IA = 5.74 ± 0.29 , β = 0.44 ± 0.03 \beta=0.44\pm0.03 β = 0.44 ± 0.03 , and a correlation coefficient of − 0.59 -0.59 − 0.59 . Details on how we determine f r f_{\mathrm{r}} f r and halo masses for our sample of cosmic shear sources, and how we treat their associated uncertainties, are provided in Appendix B. For comparisons between the alignment amplitude estimated with this model and others, we define
A I A , t o t a l m a s s d e p , ( i ) = A ^ I A f r ( i ) ( ⟨ M h ⟩ ( i ) M h , p i v ) β . (12) A_{\mathrm{IA,total}}^{\mathrm{massdep},(i)}=\widehat{A}_{\mathrm{IA}}f_{\mathrm{r}}^{(i)}\left(\frac{\langle M_{\mathrm{h}}\rangle^{(i)}}{M_{\mathrm{h,piv}}}\right)^{\beta}. \tag*{(12)} A IA , total massdep , ( i ) = A IA f r ( i ) ( M h , piv ⟨ M h ⟩ ( i ) ) β . (12)
Cosmic shear summary statistics
In cosmic shear analyses, it is most convenient to measure the signal using shear two-point correlation functions (2PCFs) while the modelling is done through shear power spectra (see Eq. (1) ). Alternatively, one could take the Fourier transform of the shear field and work in harmonic space, which brings the measurements closer to the theoretical predictions. We chose to use 2PCFs for the measurement, as they shield us from a particular form of E E E -/B B B -mode mixing that occurs in harmonic space in the presence of the (unavoidable) image masks. In addition, the shot-noise term does not impact 2PCFs (present only in their covariance), as it is only relevant at zero lag. In harmonic space, both shot noise and masking must be modelled and corrected, prior to beginning cosmological analyses. For more details see e.g. Alonso et al. [3 ] .
Our cosmological analysis here leverages two-point statistics: the complete orthogonal sets of E E E /B B B -integrals (COSEBIs, Sect. 2.5.2) and band power spectra (Sect. 2.5.3), constructed using linear combinations of finely-binned 2PCFs (Sect. 2.5.1).
2PCFs ()
The shear two-point correlation functions (ξ ± \xi_{\pm} ξ ± , Kaiser 1992) are functions of the tangential and cross components of shear (γ t \gamma_{\mathrm{t}} γ t and γ × \gamma_{\times} γ × , respectively) defined with respect to the line connecting pairs of galaxies with angular separation θ \theta θ (see e.g. Bartelmann & Schneider 2001),
ξ ± ( θ ) = ⟨ γ t γ t ⟩ ( θ ) ± ⟨ γ × γ × ⟩ ( θ ) . (13) \xi_{\pm}(\theta)=\langle\gamma_{\mathrm{t}}\gamma_{\mathrm{t}\rangle}(\theta)\pm\langle\gamma_{\times}\gamma_{\times}\rangle(\theta). \tag*{(13)} ξ ± ( θ ) = ⟨ γ t γ t ⟩ ( θ ) ± ⟨ γ × γ × ⟩ ( θ ) . (13)
The measured 2PCFs are binned in θ \theta θ to reduce measurement noise and compress the data. In practice, we measure ξ ± ( θ ) \xi_{\pm}(\theta) ξ ± ( θ ) for 1000 logarithmically spaced θ \theta θ -bins between θ min = 2 ′ \theta_{\min}=2' θ m i n = 2 ′ and θ max = 300 ′ \theta_{\max}=300' θ m a x = 30 0 ′ using TreeCorr1 ^{1} 1 ([62 ] ). The inputs to TreeCorr are galaxy catalogues, which include galaxy ellipticities (ϵ 1 \epsilon_{1} ϵ 1 and ϵ 2 \epsilon_{2} ϵ 2 ), galaxy positions (right ascension and declination), and galaxy weights (w w w ). The weights used by TreeCorr are equal to the recalibrated shape weights (Sect. 3.5) multiplied by the redshift distribution gold-weights (Sect. 3.4). Our ξ ± ( θ ) \xi_{\pm}(\theta) ξ ± ( θ ) are then estimated using
ξ ^ ± ( i j ) ( θ ˉ ) = ∑ a b w a w b [ ϵ t , a o b s ϵ t , b o b s ± ϵ × , a o b s ϵ × , b o b s ] Δ a b ( i j ) ( θ ˉ ) ∑ a b w a w b ( 1 + m a ) ( 1 + m b ) Δ a b ( i j ) ( θ ˉ ) , (14) \hat{\xi}_{\pm}^{(ij)}(\bar{\theta})=\frac{\sum_{ab}w_{a}w_{b}\left[\epsilon_{\mathrm{t},a}^{\mathrm{obs}}\epsilon_{\mathrm{t},b}^{\mathrm{obs}}\pm\epsilon_{\times,a}^{\mathrm{obs}}\epsilon_{\times,b}^{\mathrm{obs}}\right]\Delta_{ab}^{(ij)}(\bar{\theta})}{\sum_{ab}w_{a}w_{b}(1+m_{a})(1+m_{b})\Delta_{ab}^{(ij)}(\bar{\theta})}, \tag*{(14)} ξ ^ ± ( ij ) ( θ ˉ ) = ∑ ab w a w b ( 1 + m a ) ( 1 + m b ) Δ ab ( ij ) ( θ ˉ ) ∑ ab w a w b [ ϵ t , a obs ϵ t , b obs ± ϵ × , a obs ϵ × , b obs ] Δ ab ( ij ) ( θ ˉ ) , (14)
where Δ a b ( i j ) ( θ ˉ ) \Delta_{ab}^{(ij)}(\bar{\theta}) Δ ab ( ij ) ( θ ˉ ) is the bin selection function (which limits the sums to galaxy pairs of separation within the angular bin labelled by θ ˉ \bar{\theta} θ ˉ and the tomographic bins i i i and j j j ), and ϵ t o b s \epsilon_{\mathrm{t}}^{\mathrm{obs}} ϵ t obs and ϵ × o b s \epsilon_{\times}^{\mathrm{obs}} ϵ × obs are the projected tangential and radial components of observed ellipticities between pairs of galaxies. The estimator includes the galaxies’ multiplicative shear measurement biases m a m_{a} m a estimated via image simulations (Sect. 3.4).
Using the flat-sky approximation, we can relate the 2PCFs to the angular (shear) power spectra by Hankel transformations:
ξ ± ( i j ) ( θ ) = ∫ 0 ∞ ℓ d ℓ 2 π [ C ϵ ϵ , E ( i j ) ( ℓ ) ± C ϵ ϵ , B ( i j ) ( ℓ ) ] J 0 / 4 ( ℓ θ ) , (15) \xi_{\pm}^{(ij)}(\theta)=\int_{0}^{\infty}\frac{\ell\,\mathrm{d}\ell}{2\pi}\left[C_{\epsilon\epsilon,E}^{(ij)}(\ell)\pm C_{\epsilon\epsilon,B}^{(ij)}(\ell)\right]J_{0/4}(\ell\theta), \tag*{(15)} ξ ± ( ij ) ( θ ) = ∫ 0 ∞ 2 π ℓ d ℓ [ C ϵϵ , E ( ij ) ( ℓ ) ± C ϵϵ , B ( ij ) ( ℓ ) ] J 0/4 ( ℓ θ ) , (15)
where C ϵ ϵ , E / B ( i j ) ( ℓ ) C_{\epsilon\epsilon,E/B}^{(ij)}(\ell) C ϵϵ , E / B ( ij ) ( ℓ ) are the E E E /B B B -mode angular power spectra for redshift bins i i i and j j j , and J 0 / 4 ( x ) J_{0/4}(x) J 0/4 ( x ) are the zeroth- and fourth-order Bessel functions of the first kind. We expect cosmic shear to produce only E E E -mode signals up to first order (Schneider et al. 1998; Hilbert et al. [54 ] , while data systematic errors can in principle create both E E E and B B B modes. As a result, B-modes have long been used as a diagnostic tool in cosmic shear analyses (see e.g. Pen et al. [101 ] ; Hoekstra et al. [58 ] ; Heymans et al. [49 ] ; Kilbinger et al. [71 ] ).
As ξ ± \xi_{\pm} ξ ± mix E E E and B B B modes, they are not ideal for cosmic shear analysis. In addition, they variably weight a large range of physical scales in each angular bin, increasing sensitivity to scale-dependent effects such as baryonic feedback (which must be marginalised over using imperfect models). Figure 1 shows the range of ℓ \ell ℓ probed by the ξ ± \xi_{\pm} ξ ± over a range of angular separations. The J 0 J_0 J 0 filter in Eq. (15) means that ξ + \xi_+ ξ + mixes a very wide range of scales (1 ≲ ℓ ≲ 10 4 1 \lesssim\ell\lesssim10^4 1 ≲ ℓ ≲ 1 0 4 ) and J 4 J_4 J 4 makes ξ − \xi_- ξ − particularly sensitive to smaller scales (10 ≲ ℓ ≲ 10 5 10 \lesssim\ell\lesssim10^5 10 ≲ ℓ ≲ 1 0 5 ). This breadth can create a problem when using ξ ± \xi_{\pm} ξ ± for inference, as the power spectra and intrinsic alignment models become increasingly uncertain at smaller scales. Consequently, we expect the cosmological constraints to be possibly affected by a range of systematic effects that may introduce biases and, as such, we do not include ξ ± \xi_{\pm} ξ ± in our primary cosmological analyses. Nevertheless, results for ξ ± \xi_{\pm} ξ ± are provided for completeness and comparison with previous work, in Appendix F.
COSEBIs (En/Bn)
The COSEBIs (Schneider et al. [117 ] ) were originally designed to cleanly separate E E E - and B B B -modes over a finite θ \theta θ -range. This was achieved through the use of filter functions T ± n ( θ ) T_{\pm n}(\theta) T ± n ( θ ) , which form complete bases, creating well-defined COSEBIs E E E - and B B B -modes (E n / B n E_n/B_n E n / B n ). COSEBIs can be calculated or measured via a linear transformation of ξ ± \xi_{\pm} ξ ± ,
E ^ n = 1 2 ∫ θ min θ max d θ θ [ T + n ( θ ) ξ + ( θ ) + T − n ( θ ) ξ − ( θ ) ] , \hat{E}_n = \frac{1}{2} \int_{\theta_{\min}}^{\theta_{\max}} \mathrm{d}\theta\,\theta\left[T_{+n}(\theta)\xi_+(\theta) + T_{-n}(\theta)\xi_-(\theta)\right], E ^ n = 2 1 ∫ θ m i n θ m a x d θ θ [ T + n ( θ ) ξ + ( θ ) + T − n ( θ ) ξ − ( θ ) ] ,
B ^ n = 1 2 ∫ θ min θ max d θ θ [ T + n ( θ ) ξ + ( θ ) − T − n ( θ ) ξ − ( θ ) ] . (16) \hat{B}_n = \frac{1}{2} \int_{\theta_{\min}}^{\theta_{\max}} \mathrm{d}\theta\,\theta\left[T_{+n}(\theta)\xi_+(\theta) - T_{-n}(\theta)\xi_-(\theta)\right]. \tag*{(16)} B ^ n = 2 1 ∫ θ m i n θ m a x d θ θ [ T + n ( θ ) ξ + ( θ ) − T − n ( θ ) ξ − ( θ ) ] . (16)
The COSEBIs modes n n n start from 1 and go to infinity. However, the first few COSEBIs are sufficient to capture essentially all the cosmological information in each pair of tomographic bins (Asgari et al. [9 ] ). This means that COSEBIs naturally compress the data into well-defined modes, once θ min \theta_{\min} θ m i n and θ max \theta_{\max} θ m a x are defined, without any need for subsequent data compression. In previous KiDS analyses (e.g. Asgari et al. [8 ] ), the first five COSEBIs modes were used for cosmological analyses. In this work, however, we add the sixth mode, as we found that (with our increased constraining power, and larger redshift baseline) the sixth mode contains relevant cosmological information and thus improves our constraining power on S 8 S_8 S 8 . Conversely, the addition of modes 7–20 did not result in any improvement of constraining power, returning constraints identical to those measured in the 6 mode case.
We use the relation between COSEBIs and shear power spectra to make our predictions,
E n = ∫ 0 ∞ d ℓ ℓ 2 π C E E ( ℓ ) W n ( ℓ ) , E_n = \int_0^\infty\frac{\mathrm{d}\ell\,\ell}{2\pi} C_{\mathrm{EE}}(\ell)W_n(\ell), E n = ∫ 0 ∞ 2 π d ℓ ℓ C EE ( ℓ ) W n ( ℓ ) ,
B n = ∫ 0 ∞ d ℓ ℓ 2 π C B B ( ℓ ) W n ( ℓ ) , (17) B_n = \int_0^\infty\frac{\mathrm{d}\ell\,\ell}{2\pi} C_{\mathrm{BB}}(\ell)W_n(\ell), \tag*{(17)} B n = ∫ 0 ∞ 2 π d ℓ ℓ C BB ( ℓ ) W n ( ℓ ) , (17)
where the weight functions W n ( ℓ ) W_n(\ell) W n ( ℓ ) are
W n ( ℓ ) = ∫ θ min θ max d θ θ T + n ( θ ) J 0 ( ℓ θ ) = ∫ θ min θ max d θ θ T − n ( θ ) J 4 ( ℓ θ ) . (18) \begin{aligned}
W_n(\ell) = \int_{\theta_{\min}}^{\theta_{\max}} \mathrm{d}\theta\,\theta\,T_{+n}(\theta)J_0(\ell\theta) \\
&= \int_{\theta_{\min}}^{\theta_{\max}} \mathrm{d}\theta\,\theta\,T_{-n}(\theta)J_4(\ell\theta). \tag*{(18)}
\end{aligned} W n ( ℓ ) = ∫ θ m i n θ m a x d θ θ T + n ( θ ) J 0 ( ℓ θ ) = ∫ θ m i n θ m a x d θ θ T − n ( θ ) J 4 ( ℓ θ ) . (18)
Figure 1 shows two example integrands from Eq. (17) . W n ( ℓ ) W_n(\ell) W n ( ℓ ) are highly oscillatory weight functions with a limited range of support over ℓ \ell ℓ . This allows us to truncate our integrals in Eq. (17) to a finite range in ℓ \ell ℓ , as the E n E_n E n integrand in Eq. (17) is practically zero outside 10 ≲ ℓ ≲ 10 3 10 \lesssim\ell\lesssim10^3 10 ≲ ℓ ≲ 1 0 3 .
Band powers (CE,B)
Band power spectra (hereafter ‘band-powers’) refer to binned angular power spectra measured by integrating over ξ ± ( θ ) \xi_{\pm}(\theta) ξ ± ( θ ) (Schneider et al. [119 ] ; van Uitert et al. [134 ] ; Joachimi et al. 2021). Ideally, we want a finite number of bins, labelled as l l l , where each bin is only sensitive to a narrow range of ℓ \ell ℓ -values between ℓ l o , l \ell_{\mathrm{lo},l} ℓ lo , l and ℓ u p , l \ell_{\mathrm{up},l} ℓ up , l . We also want the E E E - and B B B -mode bandpowers to be pure, namely, no leakage from B B B modes into E E E modes and vice versa. For this ideal case, we can express the following:
C E / B , l ( i j ) : = 1 N l ∫ 0 ∞ d ℓ ℓ S l ( ℓ ) C ε ε , E / B ( i j ) ( ℓ ) , (19) C^{(ij)}_{E/B,l} := \frac{1}{N_l}\int_{0}^{\infty}\mathrm{d}\ell\,\ell\,S_l(\ell)\,C^{(ij)}_{\varepsilon\varepsilon,E/B}(\ell), \tag*{(19)} C E / B , l ( ij ) := N l 1 ∫ 0 ∞ d ℓ ℓ S l ( ℓ ) C εε , E / B ( ij ) ( ℓ ) , (19)
for redshift bins i i i and j j j , where S l ( ℓ ) S_l(\ell) S l ( ℓ ) is the response function for bin l l l . Here we have normalised the integral with
N l = ∫ 0 ∞ d ℓ ℓ − 1 S l ( ℓ ) . (20) N_l = \int_{0}^{\infty}\mathrm{d}\ell\,\ell^{-1}S_l(\ell). \tag*{(20)} N l = ∫ 0 ∞ d ℓ ℓ − 1 S l ( ℓ ) . (20)
In this work, we consider the ideal case of a top-hat response function between ℓ l o , l \ell_{\mathrm{lo},l} ℓ lo , l and ℓ u p , l \ell_{\mathrm{up},l} ℓ up , l which simplifies the normalisation to
N l = log ( ℓ u p , l ) − log ( ℓ l o , l ) , (21) N_l = \log(\ell_{\mathrm{up},l})-\log(\ell_{\mathrm{lo},l}), \tag*{(21)} N l = log ( ℓ up , l ) − log ( ℓ lo , l ) , (21)
and produces a band-power that traces ℓ 2 C ( ℓ ) \ell^2 C(\ell) ℓ 2 C ( ℓ ) at the logarithmic centre of the bin. We define eight band-powers logarithmically binned between ℓ min = 100 \ell_{\min}=100 ℓ m i n = 100 and ℓ max = 1500 \ell_{\max}=1500 ℓ m a x = 1500 . In practice, measuring band-powers with perfect top-hat functions is infeasible, as ξ ± ( θ ) \xi_{\pm}(\theta) ξ ± ( θ ) is only measured for a finite θ \theta θ -range. This has the adverse effect that our band-powers will inevitably mix E E E and B B B modes4 . The mixing increases for narrower ℓ \ell ℓ -bins and when the available range of θ \theta θ -scales is smaller ((Asgari & Schneider 2015)). With this definition for band-powers, we forgo perfect E E E -B B B -mode separation in favour of more control over the ℓ \ell ℓ -values that contribute to our analyses. We present the detailed derivation of our band-power windows C E , l C_{E,l} C E , l in Appendix § 10 , and show the integrands of C E , l C_{E,l} C E , l in Figure 1 . Note, however, that our integrands behave differently to those in Asgari et al. [8 ] , due to our different treatment of the edges of T ( θ ) T(\theta) T ( θ ) (Appendix § 10 ) and due to our different scale-cuts.
Analytic covariance modelling
We presented the covariance modelling for KiDS-Legacy and the OneCovariance code5 in Reischke et al. [107 ] . In general, the methodology for the cosmic shear part still follows closely the covariance modelling of KiDS-1000 (Joachimi et al. [65 ] ). Major changes to the previous analysis include: (i) the possibility of including non-uniform source distributions in the mixed term of the covariance from the triplet counts of the survey; (ii) including non-binary masks for the area calculation and survey response; and (iii) offering a more generalised framework for a consistent treatment of all summary statistics used in KiDS-Legacy.
The cosmic shear covariance, C \mathbf{C} C , generally consists of four contributions:
C = C G + C N G + C S S C + C m u l t , (22) \mathbf{C}=\mathbf{C}_{\mathrm{G}}+\mathbf{C}_{\mathrm{NG}}+\mathbf{C}_{\mathrm{SSC}}+\mathbf{C}_{\mathrm{mult}}, \tag*{(22)} C = C G + C NG + C SSC + C mult , (22)
with the Gaussian contribution G, the non-Gaussian term NG due to non-linear in-survey modes, the super-sample covariance SSC from modes larger than the survey, and lastly the contribution due to the uncertainty in the multiplicative shear bias calibration, C m u l t C_{\mathrm{mult}} C mult . Asgari et al. [8 ] ; Secco et al. [121 ] included the multiplicative shear bias in the theory prediction instead and marginalised over in the sampling process. If the residual uncertainties on m a m_a m a are small, these two approaches are equivalent. The Gaussian contribution can further be split into three parts:
C G = C s v a + C s n + C m i x , (23) \mathbf{C}_{\mathrm{G}}=\mathbf{C}_{\mathrm{sv}a}+\mathbf{C}_{\mathrm{sn}}+\mathbf{C}_{\mathrm{mix}}, \tag*{(23)} C G = C sv a + C sn + C mix , (23)
describing the sample variance term ‘sva’ due to the finite number of available modes in the survey volume, the pure shot noise term ‘sn’ due to the finite number of sources and their composition, and the mix term ‘mix’.
We calculated the covariances for angular power spectra and transform those to the three summary statistics considered. The only exception was the pure noise term, which was directly computed from the weighted number of pairs in the catalogue. SSC and NG terms were modelled with a halo model approach and the survey variance was directly calculated from the survey mask. Our analytical prescription was tested against GLASS (Tessore et al. [130 ] ) mocks to test the impact of variable depth and the overall agreement. These mocks are discussed in more detail in Appendix § 11 . In Reischke et al. [107 ] , we showed that the agreement between our analytic covariance and that constructed from GLASS is within 10 percent at large angular separations and below 5 percent for scales below 0.5 degrees. Furthermore, we found that our idealised mix term (which does not account for non-uniformity of the source distribution) did not affect our inference results significantly, and as a result was neglected in our analyses.
Data and analysis
In this section, we summarise the KiDS-Legacy lensing data and its calibration, the associated image simulations, and our analysis blinding strategy. Our analysis was performed using a new version of the CosmoPipe infrastructure presented in Wright et al. [139 ] . Our fiducial pipeline includes all steps required for a catalogues-to-cosmology analysis, leveraging three primary inputs: a (blinded) shear catalogue (Sect. 3.3), a spectroscopic compilation (Sect. 3.2), and a simulated galaxy lightcone (Sect. 3.4). With these data products, the pipeline performs nine primary functions: pre-processing of the wide-field catalogue, simulated catalogue construction, redshift calibration, N ( z ) N(z) N ( z ) construction, shear calibration, data vector construction, covariance construction, cosmological inference, and post-processing. Details of each processing step are provided alongside the public code base6 .
Figure 2 shows the footprint of the KiDS-Legacy sample, split between the two survey patches of KiDS (one in each of the northern and southern Galactic caps). The figure is coloured by the mean effective weight (i.e. including both shape measurement and redshift distribution estimation weights) of sources on sky, which is an approximate tracer of survey properties such as the size of the point-spread function and variable depth (see Appendix § 11 and Yan et al. [141 ] ). The weight is defined as described in Section 7 of Wright et al. [140 ] , and has the range w ∈ [ 0 , σ p o p − 2 ] w \in[0,\sigma_{\mathrm{pop}}^{-2}] w ∈ [ 0 , σ pop − 2 ] , where σ p o p = 0.255 \sigma_{\mathrm{pop}}=0.255 σ pop = 0.255 is an approximate intrinsic ellipticity dispersion of all sources in KiDS (without e.g. tomographic or brightness selections) defined for Legacy in Li et al. (2023b).
Figure 2. Distribution of the mean weight of the KiDS-Legacy source galaxies in on-sky bins (0.1 × 0.1 degrees). The weights here include both shape measurement weights (Sect. 3.5) and redshift distribution gold weighting (Sect. 3.4).
Notes. All values are computed for the final gold-selected sample of sources, and with the appropriate shape-measurement and gold weights applied. Where relevant values are computed using an effective area of 967.4 deg2 ^2 2 . Values of σ ϵ \sigma_{\epsilon} σ ϵ and n e f f n_{\mathrm{eff}} n eff are computed using Eqs. C.9 and C.12 of Joachimi et al. [65 ] , respectively. The σ ϵ \sigma_{\epsilon} σ ϵ values are the ellipticity dispersion averaged over both components. m g o l d m_{\mathrm{gold}} m gold is the multiplicative shear measurement bias, averaged over both ellipticity components, estimated from our SKiLLS simulations. σ ϵ \sigma_{\epsilon} σ ϵ values include the appropriate correction for multiplicative shear biases (i.e. the m g o l d m_{\mathrm{gold}} m gold column).
KiDS and VIKING
The KiDS-Legacy lensing sample was drawn from the fifth and final data release of the KiDS survey7 , KiDS-dr5, described in Wright et al. [140 ] . It includes 1347 square-degree tiles observed in the u u u , g g g , r r r , and i i i filters with the VST (Capaccioli et al. 2005) and OmegaCAM (Kuijken [78 ] ) between 2011 and 2019, as well as the 2009–2018 VIKING survey VISTA/VIRCAM observations (Edge et al. [32 ] ) of the same area in Z Z Z , Y Y Y , J J J , H H H , and K s K_{\mathrm{s}} K s . Additionally, the dataset contains spectroscopic calibration fields (the ‘KiDZ’ fields, 27 square-degree tiles in total, including four that are part of the main survey).
KiDS-dr5 significantly extends the previous (fourth) data release (Kuijken et al. [79 ] ), which was the basis for the KiDS-1000 analyses ([8 , 51 , 132 ] ). It represents a substantial increase in survey area (by 34%) and in i i i -band depth (doubling the exposure time with a second pass). The volume of photometric redshift calibration data has also grown, including photometry for 5 times more spectroscopic sources than in the original KiDS-1000 analysis (Hildebrandt et al. [56 ] ). The deeper i i i -band depth and more extensive calibration data has allowed us to push the photometric redshift limit for galaxies in the lensing sample from 1.2 to 2. The combined effect of the increase in area and depth is an increase in the survey volume of a factor of 3.5.
Several aspects of the data analysis were changed between DR4 and DR5, as detailed in Wright et al. [140 ] . The astrometry reference catalogue for the THELI reduction moved from the Sloan Digital Sky Survey (SDSS Ahumada et al. [2 ] ) and the 2-micron All Sky Survey (2MASS Skrutskie et al. [125 ] ) to Gaia DR2 (Gaia Collaboration: Brown et al. [41 ] ), and the astrometric solutions and photometric zero-point determinations in the Astro-WISE u u u -band reductions were made more robust. Also, the manual and automated masking of artefacts in the images, mostly due to ghost reflections, artificial satellite tracks, and detector instabilities, was revisited and improved for the final release. Most importantly, the weak lensing shape measurement method uses an updated version of lensfit (v321), with improvements in the sampling of the ellipticity likelihood surface and of the treatment of residual PSF leakage in the shapes and associated weights (Li et al. 2023b). Details and tests of these new implementations are given in Wright et al. [140 ] , and in the sections below. A summary of the properties of the KiDS-Legacy sample is given in Table 1 , including details of tomographic bin limits, final effective number densities of galaxies (after redshift distribution calibration), and redshift and shape measurement bias parameters.
Table 1. Properties of the six KiDS-Legacy tomographic bins.
Lastly, we note that a pernicious astrometric issue was identified during the investigation of our B-mode null tests (Sect. 4.3), which resulted in additional masking of approximately 4% of the survey sources to be used in the KiDS-Legacy analysis. More details are given in Sect. 4.3 and Appendix E. The total area after masking is 967.4 deg2 ^2 2 , with n e f f n_{\mathrm{eff}} n eff = 8.79 per arcminute within the photometric redshift limits of KiDS-Legacy.
KiDZ
The KiDZ data consist of VST and VISTA images that target well-studied spectroscopic survey fields, including the Cosmic Evolutions Survey field (COSMOS; Scoville et al. [120 ] ), the DEEP2 Galaxy Redshift Survey (DEEP2; Newman et al. [98 ] ) 02h and 23h fields, the Chandra Deep Field South (CDFS), the Visible Multi-object Spectrograph (VIMOS) Very Large Telescope (VLT) Deep Survey (VVDS; Le Fèvre et al. [82 ] ) 14h field, and the VIMOS Public Extragalactic Redshift Survey (VIPERS; Guzzo et al. [47 ] ) W1 and W4 fields (a full list is given in Wright et al. [140 ] ). They were purposely taken under very similar observing conditions, and with the same observing setup, as the main KiDS/VIKING survey and serve as the calibration sample for photometric redshifts of the faint KiDS sources. A total of 126 085 sources with spectroscopic redshifts have nine-band KiDZ photometry. Additional photometric redshift calibration via cross-correlations uses the overlap between KiDS and shallower, wide-angle spectroscopic redshift surveys, principally: SDSS, the Galaxy And Mass Assembly (Driver et al. [30 ] ) survey, 2-degree Field Lensing Survey (Blake et al. [16 ] ), WiggleZ survey (Blake et al. [17 ] ), and the Dark Energy Spectroscopic Instrument (DESI Collaboration et al. [29 ] ).
Blinding
We used the same blinding methodology as the one implemented in Asgari et al. [8 ] , first presented in Kuijken et al. [77 ] . This process involves the double-blind generation of three realisations of our shape measurement catalogues. One of these is the true catalogue, while in the other two there are induced systematic differences in measured shapes that result in (up to) ± 2 σ \pm2\sigma ± 2 σ change in the inferred S 8 S_8 S 8 . The true catalogue is equally likely to be the one that yields the highest, middle, or lowest S 8 S_8 S 8 . These blinded data vectors were analysed throughout the project until the analyses were finalised, cosmological chains were run, the results written up in this publication, and the manuscripts had undergone internal review within the KiDS collaboration. As such, this manuscript (for example) was written containing results for all three blinds, with only the discussion and conclusion sections written after unblinding.
The authors note (for the reference of future surveys) that this blinding scheme effectively hampered the diagnosis of subtle systematic effects in the dataset (see Appendix E), as we were limited in our ability to directly compare measured shapes between previously constructed catalogues and our KiDS-Legacy catalogues. Such comparisons are extremely useful for the validation of spatially localised systematic patterns in the shape catalogue, which are obfuscated when translated into two-point statistics. This blinding strategy, however, makes such direct comparisons between shape catalogues problematic, as it can also lead to unblinding.
Redshift distributions and calibration
Redshift distributions (hereafter ‘N ( z ) N(z) N ( z ) ’) used in KiDS-Legacy were estimated and calibrated leveraging both colour-based direct calibration using self-organising maps (SOMs Kohonen [76 ] ) and spatial cross-correlations (CC) between tracer spectroscopic samples and the KiDS-Legacy sources. For a detailed description of the N ( z ) N(z) N ( z ) construction and validation process used in this work, we refer to Wright et al. [138 ] . In brief, the N ( z ) N(z) N ( z ) here differ in their construction (compared to previous work from KiDS) in three main regards: a considerably larger redshift calibration sample from KiDZ (see Sect. 3.2); a new weighting scheme that replaces the binary ‘gold-class’ with a continuous ‘gold-weight’; and through the optional inclusion of weights for spectroscopic objects, such as our so-called ‘prior-volume weighting’. Gold-weights primarily account for noise in the colour-based association between wide-field sources and their calibrating spectra, and are computed by repeating the gold classification N N N times (where we arbitrarily chose N = 10 N = 10 N = 10 ). Similarly, our prior-volume weighting applies an a priori weighting to the distribution of calibrating spectra as a function of redshift, to reduce bias caused by colour-redshift degeneracies and the complex redshift selection function of our spectroscopic compilation.
The redshift distribution calibration for KiDS-Legacy made use of two suites of simulations. For the calibration of the SOM redshift distributions (and calibration of source shape measurements, Sect. 3.5), we used the multi-colour SKiLLS simulation of Li et al. (2023b). For additional SOM redshift distribution calibration, and cross-correlation redshift distribution calibration, we used the MICE2 simulated lightcone (Fosalba et al. [39 ] ; Crocce et al. [26 ] ; Fosalba et al. [40 ] ; Carretero et al. [23 ] ; Hoffmann et al. [59 ] ) adapted to KiDS noise levels and sample selections (van den Busch et al. [133 ] ).
Figure 3 shows the fiducial redshift distributions used in this cosmological analysis. These N ( z ) N(z) N ( z ) include the appropriate δ z \delta z δ z shifts in each of their mean, as inferred using our simulations and described in Wright et al. [138 ] . As such, the N ( z ) N(z) N ( z ) presented are shown as-used in the generation of the theory/models during parameter inference. In Wright et al. [138 ] we quantified the impact of a series of different analysis choices for N ( z ) N(z) N ( z ) estimation and calibration, demonstrating that they have no significant impact on the final N ( z ) N(z) N ( z ) that we used for cosmological inference. The figure also presents the off-diagonal elements of the correlation matrices for δ z \delta z δ z , used to construct our correlated δ z \delta z δ z priors. The lower triangle shows the correlation matrix estimated from our SOM redshift distribution calibration, while the upper triangle shows the correlation matrix for the CC estimated redshift distributions.
Figure 3. Left: N ( z ) N(z) N ( z ) used in the fiducial cosmic shear analyses of KiDS-Legacy. Each panel is one tomographic bin whose definition is determined by the grey band in photo-z z z . The N ( z ) N(z) N ( z ) have been shifted according to the mean bias estimated using our SKiLLS simulations, making each N ( z ) N(z) N ( z ) a correct representation of the effective N ( z ) N(z) N ( z ) that is used in the cosmological analyses. Fiducial biases are provided in Table 1 . Right: δ z \delta z δ z correlation matrix for our fiducial N ( z ) N(z) N ( z ) (δ z S O M \delta z^{\mathrm{SOM}} δ z SOM , lower triangle ) and the correlation matrix inferred from our cross correlation N ( z ) N(z) N ( z ) on the data (D z C C D_z^{\mathrm{CC}} D z CC , upper triangle ).
Shape measurement and calibration
The measurement and calibration of source shapes in KiDS-Legacy were performed with the latest KiDS methods, which have already been established in the literature. We refer to Li et al. (2023b) for an extensive description of the state-of-the-art for shape measurement and calibration in KiDS. Topics discussed therein include: the lensfit algorithm (v321; see also Miller et al. [97 , 96 ] ; Fenech Conti et al. [35 ] ), used for measuring shapes of galaxies (including modelling of the PSF); shape calibration using SKiLLS image simulations; and PSF leakage corrections (with minor updates to the PSF leakage estimation as described in Wright et al. [140 ] ). The SKiLLS image simulations are multi-band image simulations, constructed specifically to match the observational properties of KiDS-Legacy in all available bands, and featuring realistic clustering and correlations between galaxy properties and environment, allowing for realistic treatment of blending effects. The simulations also feature redshift-dependent shear (MacCrann et al. [90 ] ), although this is determined to be of negligible importance to the computation of multiplicative bias ((Li et al. 2023b)) and redshift distribution bias ([138 ] ) at the sensitivity of KiDS-Legacy.
The multiplicative bias estimation in KiDS-Legacy was performed using the same suite of multi-colour image simulations from SKiLLS used for redshift distribution calibration. Redshift distribution calibration utilises many realisations of the sample to estimate uncertainties on tomographic bias parameters (see [138 ] ). Afterwards, these samples were also used to estimate the multiplicative shear measurement bias, m m m , per tomographic bin. This means that, for the first time, we have been able to directly estimate the covariance between our redshift distribution bias calibration and our multiplicative shear bias estimation. The resulting covariance matrix (and the corresponding correlation matrix) are functionally block-diagonal: there is no cross covariance between m m m and δ z \delta z δ z estimates.
We checked for the degree of PSF leakage measured in the lensing sample through a direct empirical calculation of the correlation between PSF and galaxy ellipticities. Here, we document the measured PSF contamination in our KiDS-Legacy sample, and the measured residual one- and two-dimensional constant ellipticity terms (‘c c c -terms’) in the dataset. We adopted the first-order systematics model from Heymans et al. [52 ] :
ε k o b s = ( 1 + m k ) ( ε k i n t + γ k ) + α k ε k P S F + c k , k ∈ { 1 , 2 } , (24) \varepsilon_k^{\mathrm{obs}} = (1 + m_k)(\varepsilon_k^{\mathrm{int}} + \gamma_k) + \alpha_k\varepsilon_k^{\mathrm{PSF}} + c_k,\quad k \in\{1,2\}, \tag*{(24)} ε k obs = ( 1 + m k ) ( ε k int + γ k ) + α k ε k PSF + c k , k ∈ { 1 , 2 } , (24)
where k k k denotes the individual ellipticity components, γ \gamma γ is the (true) shear, and the superscripts ‘obs’, ‘int’, and ‘PSF’ refer to a source’s observed, intrinsic, and modelled PSF ellipticity, respectively. In the limit of large numbers of sources averaged over a large sky area, we can assume ⟨ ( 1 + m k ) ( ε k i n t + γ k ) ⟩ = 0 \langle(1+m_k)(\varepsilon_k^{\mathrm{int}}+\gamma_k)\rangle=0 ⟨( 1 + m k ) ( ε k int + γ k )⟩ = 0 , and therefore model the residual PSF contamination using a simple linear regression between observed source ellipticities (i.e. after recalibration to remove PSF leakage [140 ] (Li et al. 2023b)) and modelled PSF ellipticities.
Figure 4. PSF leakage measured in the KiDS-Legacy lensing sample after recalibration of shapes, computed per tomographic bin and ellipticity component k k k , and following the first-order systematics model in Eq. (24) .
Figure 4 shows the α k \alpha_k α k and c k c_k c k PSF leakage terms in each ellipticity component of the KiDS-Legacy sample, measured in each tomographic bin. Note that the shape recalibration process used in KiDS-Legacy necessitates that the measured PSF leakage factors be approximately zero, as the recalibration specifically addresses any residual in this parameter, in bins of source signal-to-noise ratio (SN) and resolution (Li et al. 2023b). The residual α k \alpha_k α k are all negligible, with the α 2 \alpha_2 α 2 having a slightly larger residual amplitude overall (α 2 ≈ 0.01 \alpha_2 \approx0.01 α 2 ≈ 0.01 , roughly three orders of magnitude smaller than that measured before any recalibration effort). Furthermore, we note that the constant term shown in the figure (c k c_k c k ) is not the same as the additive bias which is removed from our galaxy sample. In practice, the additive component of our ellipticities was computed directly (using a weighted mean) from the sample of galaxies per tomographic bin and hemisphere: c ^ k = ⟨ ε k ⟩ \hat{c}_k = \langle\varepsilon_k\rangle c ^ k = ⟨ ε k ⟩ . The c c c -terms that were subtracted from our recalibrated shape catalogues prior to computation of correlation functions are provided in Table 2 . In contrast, the figure presents the constant component of our linear regression (Eq (24) ) to all data per tomographic bin after removal of the additive bias. As such it is expected that the c k c_k c k -terms are small, but not necessarily zero.
Table 2. Subtracted additive shear per ellipticity component, tomographic bin, and hemisphere.
As a further quality test, we computed the non-tomographic 2D c c c -terms in camera coordinates per ellipticity component, which is shown in Figure 5 . The 2D c c c -terms were calculated as a weighted mean per bin (using both shape- and gold-weights), with the resulting map smoothed with a Gaussian with one-bin standard deviation (approximately 1′.2), and the computed c c c -term scaled by the (also smoothed) measured ellipticity dispersion in each bin. The residual c c c -term can be seen to have amplitudes ∣ c i / σ ϵ i ∣ ≤ 1 % |c_i/\sigma_{\epsilon_i}| \le1\% ∣ c i / σ ϵ i ∣ ≤ 1% , indicating that these 2D patterns are negligible in the context of our cosmic shear measurements. The primary structure that is visible in c 1 c_1 c 1 was identified in Hildebrandt et al. [55 ] as being due to an electronic effect that introduces an anomalous signal in the read-out direction of one detector on OmegaCAM. However, as shown in Asgari et al. [8 ] , the amplitude of this signal (≲ 1 % \lesssim1\% ≲ 1% of the ellipticity dispersion) is negligible in the context of our cosmic shear measurements.
Figure 5. Average of each ellipticity component in bins of the focal plane, for sources in all tomographic bins, relative to the (weighted) dispersion of ellipticities in the same bin. The weighted ellipticity dispersion per component over the bins in the focal plane are 0.287 ± 0.001 0.287 \pm0.001 0.287 ± 0.001 and 0.288 ± 0.001 0.288 \pm0.001 0.288 ± 0.001 , indicating that the observed pattern is not driven by variable ellipticity dispersion across the focal plane. The figures still indicate the presence of the two-dimensional pattern reported in Figure 2 of Hildebrandt et al. [55 ] .
Null tests
In this section, we outline the various null tests that were completed prior to the unblinding of the cosmological analysis of KiDS-Legacy; the reader interested in cosmological parameter estimates only may wish to skip to § 5 . We performed three primary tests: an analysis of the Paulin-Henriksson et al. [100 ] systematics model, which quantifies the impact of PSF modelling residuals on correlation functions; an analysis of Bacon et al. [11 ] statistics, which test specifically for contamination of the galaxy shape measurements by the PSF; and an analysis of B B B modes, which tests for a systematic signal that cannot be generated by cosmic shear.
paulin-henriksson/etal:2008 model
To validate the fidelity of our shape measurement catalogue, we first compute the Paulin-Henriksson et al. [100 ] systematics model for KiDS-Legacy, which aims to quantify the impact of PSF modelling residuals on cosmic shear correlation functions. This systematics model, in the context of the shear two-point
correlation function, reduces to the following:
⟨ ε o b s ε o b s ⟩ ≃ ( 1 + 2 [ δ T P S F T g a l ] ) ⟨ ε o b s p e r f e c t ε o b s p e r f e c t ⟩ + [ 1 T g a l ] 2 ⟨ ( ε P S F δ T P S F ) ( ε P S F δ T P S F ) ⟩ + 2 [ 1 T g a l ] 2 ⟨ ( ε P S F δ T P S F ) ( δ ε P S F T P S F ) ⟩ + [ 1 T g a l ] 2 ⟨ ( δ ε P S F T P S F ) ( δ ε P S F T P S F ) ⟩ , (25) \begin{aligned}
\langle\varepsilon_{\mathrm{obs}}\varepsilon_{\mathrm{obs}}\rangle&\simeq\left(1+2\left[\frac{\delta T_{\mathrm{PSF}}}{T_{\mathrm{gal}}}\right]\right)\langle\varepsilon_{\mathrm{obs}}^{\mathrm{perfect}}\varepsilon_{\mathrm{obs}}^{\mathrm{perfect}}\rangle\\
&\quad+ \left[\frac{1}{T_{\mathrm{gal}}}\right]^2 \langle(\varepsilon_{\mathrm{PSF}}\,\delta T_{\mathrm{PSF}})(\varepsilon_{\mathrm{PSF}}\,\delta T_{\mathrm{PSF}})\rangle\\
&\quad+ 2\left[\frac{1}{T_{\mathrm{gal}}}\right]^2 \langle(\varepsilon_{\mathrm{PSF}}\,\delta T_{\mathrm{PSF}})(\delta\varepsilon_{\mathrm{PSF}}\,T_{\mathrm{PSF}})\rangle\\
&\quad+ \left[\frac{1}{T_{\mathrm{gal}}}\right]^2 \langle(\delta\varepsilon_{\mathrm{PSF}}\,T_{\mathrm{PSF}})(\delta\varepsilon_{\mathrm{PSF}}\,T_{\mathrm{PSF}})\rangle, \tag*{(25)}
\end{aligned} ⟨ ε obs ε obs ⟩ ≃ ( 1 + 2 [ T gal δ T PSF ] ) ⟨ ε obs perfect ε obs perfect ⟩ + [ T gal 1 ] 2 ⟨( ε PSF δ T PSF ) ( ε PSF δ T PSF )⟩ + 2 [ T gal 1 ] 2 ⟨( ε PSF δ T PSF ) ( δ ε PSF T PSF )⟩ + [ T gal 1 ] 2 ⟨( δ ε PSF T PSF ) ( δ ε PSF T PSF )⟩ , (25)
where we only keep terms that correlate to first order. It should be noted that we use a shorthand notation in Eqs. (25) , (26) , and (27) as in § 3.3 of Giblin et al. [44 ] , so that ⟨ a b ⟩ \langle ab\rangle ⟨ ab ⟩ denotes ξ ± \xi_{\pm} ξ ± . Here, ϵ o b s \epsilon_{\rm obs} ϵ obs is the observed ellipticity of a source, ϵ o b s p e r f e c t \epsilon_{\rm obs}^{\rm perfect} ϵ obs perfect is the noiseless and unbiased measurement of source ellipticity (which we approximate as being equal to the modelled shear: ϵ o b s p e r f e c t = γ t h e o r y \epsilon_{\rm obs}^{\rm perfect}=\gamma_{\rm theory} ϵ obs perfect = γ theory ), and T T T is the R 2 R^2 R 2 object size estimate.
Figure 6 shows the various contributions to the additive systematic, δ ξ + s y s ( θ ) \delta\xi_{+}^{\rm sys}(\theta) δ ξ + sys ( θ ) , from the Paulin-Henriksson et al. [100 ] systematics model. The four terms from Eq. (25) , shown in various colours, cause ⟨ ϵ o b s ϵ o b s ⟩ \langle\epsilon_{\rm obs}\epsilon_{\rm obs}\rangle ⟨ ϵ obs ϵ obs ⟩ to deviate from ⟨ ϵ o b s p e r f e c t ϵ o b s p e r f e c t ⟩ \langle\epsilon_{\rm obs}^{\rm perfect}\epsilon_{\rm obs}^{\rm perfect}\rangle ⟨ ϵ obs perfect ϵ obs perfect ⟩ . The total systematic δ ξ + s y s \delta\xi_{+}^{\rm sys} δ ξ + sys :
Figure 6. Contributions to the additive systematic δ ξ + s y s ( θ ) \delta\xi^{\rm sys}_{+}(\theta) δ ξ + sys ( θ ) from the Paulin-Henriksson et al. [100 ] systematics model. The various lines show the four components, with their summation in thick blue-green. The yellow band depicts 10% of the cosmic shear uncertainty on ξ + \xi_{+} ξ + . The figure is linear between the black dashed lines, and logarithmic outside. This demonstration is for tomographic bin six, but all bins show quantitatively similar results.
δ ξ + s y s ( θ ) ≔ ⟨ ϵ o b s ϵ o b s ⟩ − ⟨ ϵ o b s p e r f e c t ϵ o b s p e r f e c t ⟩ (26) \delta\xi_{+}^{\rm sys}(\theta) \coloneqq\langle\epsilon_{\rm obs}\epsilon_{\rm obs}\rangle-\left\langle\epsilon_{\rm obs}^{\rm perfect}\epsilon_{\rm obs}^{\rm perfect}\right\rangle \tag*{(26)} δ ξ + sys ( θ ) : = ⟨ ϵ obs ϵ obs ⟩ − ⟨ ϵ obs perfect ϵ obs perfect ⟩ (26)
is given by the summation of these four terms8 and can be compared to the yellow band which encloses 10% of the uncertainty on the cosmic shear correlation function ξ + \xi_{+} ξ + , from the OneCovariance (Sect. 2.6). The figure presents the analysis of our sixth tomographic bin (quantitatively similar results are obtained for all other bins), and demonstrates that errors in our PSF modelling are sufficiently small (consistently less than 10% of our correlation function uncertainty) that there is no significant contamination of our two-point statistics from this source of bias.
As such, we conclude that PSF modelling errors are negligible for the analysis of cosmic shear in KiDS-Legacy.
Bacon et al. (2003) statistics
We also used the method from Bacon et al. [11 ] to infer the degree of PSF leakage into the two-point correlation function measurement using
ξ ± s y s = ⟨ ϵ o b s ϵ P S F ⟩ 2 ⟨ ϵ P S F ϵ P S F ⟩ . (27) \xi_{\pm}^{\rm sys}=\frac{\langle\epsilon_{\rm obs}\epsilon_{\rm PSF}\rangle^2}{\langle\epsilon_{\rm PSF}\epsilon_{\rm PSF}\rangle}. \tag*{(27)} ξ ± sys = ⟨ ϵ PSF ϵ PSF ⟩ ⟨ ϵ obs ϵ PSF ⟩ 2 . (27)
In Figure 7 , we show the value of ξ ± s y s \xi_{\pm}^{\rm sys} ξ ± sys in each tomographic bin autocorrelation, as a fraction of the theory prediction of the cosmic shear signal from our best-fit ξ ± \xi_{\pm} ξ ± cosmology (ξ Λ C D M \xi^{\Lambda{\rm CDM}} ξ Λ CDM ). Results from the cross-correlation bins are quantitatively similar to the autocorrelation. As a reference, we also show 10% of the 1 σ 1\sigma 1 σ uncertainty on the cosmic shear signal as reported by the O
Figure 7. Ratio of PSF contamination in the cosmic shear correlation function, ξ + s y s \xi_{+}^{\mathrm{sys}} ξ + sys , to the predicted best-fit cosmic shear correlation function, ξ + Λ C D M \xi_{+}^{\Lambda\mathrm{CDM}} ξ + Λ CDM . The statistic is shown relative to 10% of the 1 σ 1\sigma 1 σ uncertainty on our fiducial ξ + \xi_{+} ξ + correlation functions (blue). We show only tomographic bin auto-correlations here, as the cross-correlation signals are extremely similar to the auto-correlations (i.e. combinations i − j i-j i − j ; ∀ j > i \forall j>i ∀ j > i are quantitatively similar to the i − i i-i i − i signal).
neCovariance (Sect. 2.6). The measured fractional PSF contamination is less than the benchmark 10% of the cosmic shear correlation function uncertainty for all but five data points (out of 189 across all tomographic bin combinations). As such, we conclude that the PSF contamination of cosmic shear correlation functions is negligible in KiDS-Legacy, as expected from our PSF leakage analysis (Sect. 3.5).
B modes
Our third null test examined the significance of B B B -mode signals in the data vector used for cosmological analyses. This test was based on the understanding that, at the sensitivity available to stage-III cosmological imaging surveys, B B B -mode signals of cosmological and astrophysical origin ought to be negligible. As such, any significant detection of B B B modes implies a contamination of the data vector by systematic effects, which may similarly (but undiagnosably) also impact the E E E -modes used for cosmological inference.
We used a tomographic analysis, equivalent to that of the fiducial cosmological analysis, to derive the significance of B B B modes (estimated via the COSEBIs B n B_n B n ) in our data vector. In the early stages of our blinded cosmological analysis, it became clear that our initial KiDS-Legacy sample failed this null test, at high significance. This led to a re-examination of the cosmological sample, which in turn led to a redefinition of the survey footprint to conservatively exclude 46.6 square-degrees with increased astrometric scatter between exposures in our lensing imaging (see Appendix E).
The B B B -mode signal is presented in Figure 8 . The p p p -values for the measured B B B modes across all tomographic bin combinations are 0.04 in the case of 6 B n B_n B n modes, and 0.09 in the case of 20 modes, and thus both pass our required threshold of p > 0.01 p>0.01 p > 0.01 (adopted from previous work within both the KiDS and DES teams). Looking at the individual tomographic bins, only the sixth bin autocorrelation, with 6 B n B_n B n modes, has a p p p -value below 0.01. However, with 21 tests (one per bin combination) it is expected that one will occasionally find p < 0.01 p<0.01 p < 0.01 even with correlated data drawn from the null hypothesis; indeed, this is essentially the calculation that is performed in the computation of the ‘all bins’ p p p -value, which passes.
Figure 8. COSEBIs B B B modes measured in KiDS-Legacy. Each panel is annotated with the p p p -value of the B B B -mode signals for twenty modes (black) and our fiducial 6 modes (blue). The total p p p -value for the full data vector, also for six and twenty modes, is annotated in the upper right corner of the figure. The B B B -mode signals are highly correlated across modes within a tomographic bin, so readers are cautioned against so-called ‘χ \chi χ -by-eye’.
Results
In this section, we present the cosmological results from KiDS-Legacy split between our fiducial results (Sect. 5.1), and variations of analysis and modelling choices (Sect. 5.2). Furthermore, we refer the reader to our companion paper (Stölzner et al. [128 ] ) for detailed analyses of the internal consistency of KiDS-Legacy, and for analysis of the external consistency between KiDS-Legacy and datasets from DES, DESI, Pantheon+ (Brout et al. [20 ] ), and Planck.
Fiducial constraints
The fiducial analysis of KiDS-Legacy used two statistics (E n E_n E n , C E C_E C E ) measured over 2 ′ .0 ≤ θ ≤ 300 ′ 2'.0 \leq\theta\leq300' 2 ′ .0 ≤ θ ≤ 30 0 ′ , an NLA-M M M intrinsic alignment model, correlated δ z \delta z δ z prior marginalisation, emulated CosmoPower P ( k ) P(k) P ( k ) trained on camb and HMCode2020, covariances constructed with the OneCovariance and augmented with the uncertainty in the shear calibration, and sampling performed by nautilus (Lange [81 ] ) within CosmoSIS, using a Gaussian likelihood. We chose the nautilus sampler due to its superior efficiency, robustness of uncertainties, and accurate evidence estimation (see e.g. Lemos et al. [83 ] ; Dark Energy Survey and Kilo-Degree Survey Collaboration et al. Notes. Parameters in the first block are inputs to the P ( k ) P(k) P ( k ) emulator, and all use flat priors denoted by their lower and upper limits [lower, upper]. We sample over five cosmological parameters: S 8 S_8 S 8 , Ω c \Omega_{\rm c} Ω c and Ω b \Omega_{\rm b} Ω b (the density parameters for cold dark matter and baryonic matter, respectively); h h h (the dimensionless Hubble parameter); and n s n_s n s (the spectral index of the primordial power spectrum). Additional inputs to the power spectrum emulator include four fixed cosmological parameters: Ω K \Omega_K Ω K (the energy density parameter associated with spatial curvature); w 0 w_0 w 0 and w a w_a w a (describing the equation of state of dark energy and its redshift evolution); and ∑ m ν \sum m_\nu ∑ m ν (the total mass of all families of neutrinos). We also sample over the baryon feedback parameter from HMCode2020 log 10 ( T A G N / K ) \log_{10}(T_{\rm AGN}/K) log 10 ( T AGN / K ) . Parameters in the second block are implemented within additional Cosmosis modules, and are used in the intrinsic alignment model (A I A A_{\rm IA} A IA , b I A b_{\rm IA} b IA , and log 10 ( ⟨ M ⟩ ( i ) ) \log_{10}(\langle M\rangle^{(i)}) log 10 (⟨ M ⟩ ( i ) ) ) or our redshift distribution bias mitigation (δ z ( i ) \delta z^{(i)} δ z ( i ) ). These sets of parameters are specified with multivariate Gaussian priors denoted by N ( μ x , C x ) \mathcal{N}(\mu_x,C_x) N ( μ x , C x ) , where μ \mu μ is the vector of prior means and C C C is the related covariance matrix. Values of μ I A , x \mu_{\rm IA,x} μ IA , x are given in Appendix B. Values of μ z \mu_z μ z , diag ( C δ z ) \operatorname{diag}(C_{\delta z}) diag ( C δ z ) are given in Table 1 , and the correlation matrix of C δ z C_{\delta z} C δ z is shown in Figure 3 .
2023). Priors for the various cosmological and nuisance parameters are given in Table 3 . Correlated priors used for our redshift distribution bias parameters are given in Sect. 3.4, and for the IA model are given in Appendix B.
Figure 9 presents the fiducial marginal constraints in σ 8 \sigma_8 σ 8 , S 8 S_8 S 8 , and Ω m \Omega_{\rm m} Ω m . Marginal constraints over both S 8 S_8 S 8 and Σ 8 = σ 8 ( Ω m / 0.3 ) α \Sigma_8 = \sigma_8(\Omega_{\rm m}/0.3)^\alpha Σ 8 = σ 8 ( Ω m /0.3 ) α (with the exponent fitted to the Ω m \Omega_{\rm m} Ω m vs. σ 8 \sigma_8 σ 8 contour) are presented in Table 4, for three estimation types: marginal mean and standard deviation, marginal mode and 1D 68% highest posterior density interval (HPDI), and maximum a posteriori (MAP) and 68% projected joint-distribution highest posterior density interval (PJ-HPD). Dark Energy Survey and Kilo-Degree Survey Collaboration et al. (2023, Sect. 2.6) present a discussion of the benefits and detriments of each metric. We also report the χ 2 \chi^2 χ 2 of our cosmological model evaluated at the MAP, and the model probability to exceed (PTE). The MAP was evaluated using an adaptive Nelder-Mead (downhill simplex) algorithm, as described in Joachimi et al. [65 ] .
Figure 9. Fiducial Ω m , S 8 \Omega_{\mathrm{m}}, S_8 Ω m , S 8 (left), and Ω m , σ 8 \Omega_{\mathrm{m}}, \sigma_8 Ω m , σ 8 (right) constraints for our fiducial two-point statistics (E n E_n E n and C E C_E C E , orange and purple respectively), compared to the CMB results from Planck–Legacy (red). The grey contours outline the marginal distribution of our a priori volume, with (inner) 1 σ 1\sigma 1 σ and (outer) 2 σ 2\sigma 2 σ boundary levels.
Data vectors for the fiducial E n E_n E n and C E C_E C E are given in Figure 10 , alongside translated posterior distributions [75 for the first application of these consistency metrics to cosmic shear datasets] drawn from the fiducial posterior computed for each statistic.
Figure 10. Data vector and TPDs for the fiducial analysis of E n E_n E n (upper panel) and C E C_E C E (lower panel). TPDs are shown as polygons spanning the 68th (red), and 95th (orange) percentiles of the posterior models. The E n E_n E n signals are highly correlated across modes within a tomographic bin, so readers are cautioned against so-called ‘χ-by-eye’. The bin 2 autocorrelation signal, for example, is consistent with the best fit model (PTE = 0.09), despite the apparent divergence of the data from the model. The overall PTEs for the full tomographic dataset in each statistic are provided in Table 4 .
Our fiducial cosmological results, from COSEBIs, estimate a value of S 8 = 0.813 ± 0.018 S_8 = 0.813 \pm0.018 S 8 = 0.813 ± 0.018 , evaluated using the marginal mean and standard deviation. This is the tightest constraint made with our COSEBIs measurements, although all constraints are relatively consistent in their constraining power: the marginal mode and HPDI estimate is S 8 = 0.81 5 − 0.021 + 0.016 S_8 = 0.815^{+0.016}_{-0.021} S 8 = 0.81 5 − 0.021 + 0.016 , and the maximum a posteriori and PJ-HPD estimate is S 8 = 0.81 1 − 0.015 + 0.022 S_8 = 0.811^{+0.022}_{-0.015} S 8 = 0.81 1 − 0.015 + 0.022 . The consistency between the measurements is a reflection of the Gaussianity of the marginal posterior, and of the N N N -dimensional posterior (in the case of the PJ-HPD). This despite the necessary smoothing of the posterior distribution that is required of the HPDI computation, which is performed using a Gaussian kernel with standard deviation defined using the method of Sheather & Jones [123 ] . Testing of the HPDI method indicates that, for our posterior distributions in S 8 S_8 S 8 , the chosen smoothing kernel is unlikely to introduce additional spurious uncertainty to our cosmological parameter constraints. As such, we expect that the constraints returned from our various summarisation methods to be robust. In the following, all values and uncertainties correspond to the marginal mode and HPDI estimate, unless otherwise specified.
The above discussion about the stability of our summarisation statistics is contrasted by the observed variance between our S 8 S_8 S 8 constraints from different statistics. The variance in S 8 S_8 S 8 between our statistics is seemingly large, with our C E C_E C E and E n E_n E n differing by 0.55 σ 0.55\sigma 0.55 σ . However, examination of the 2D posterior distribution in Ω m \Omega_{\rm m} Ω m and S 8 S_8 S 8 (Fig. 9) indicates that the optimal degeneracy direction for KiDS-Legacy is somewhat different than the assumed α = 0.5 \alpha= 0.5 α = 0.5 parametrised by S 8 S_8 S 8 . This is expected, given the different scales probed by each of our statistics, and has the effect of exacerbating apparent differences between our statistics when marginalising on the S 8 S_8 S 8 parameter (as translation along the degeneracy can give the false impression of internal inconsistency between final parameter constraints).
When fitting for the optimal α \alpha α under the more general Σ 8 = σ 8 ( Ω m / 0.3 ) α \Sigma_8 = \sigma_8(\Omega_{\rm m}/0.3)^\alpha Σ 8 = σ 8 ( Ω m /0.3 ) α , we found the best fitting α \alpha α is 0.58 for E n E_n E n and 0.60 for C E C_E C E . These significant deviations from α = 0.5 \alpha= 0.5 α = 0.5 can cause offsets between derived values of S 8 S_8 S 8 , as is seen between our C E C_E C E and E n E_n E n constraints. Additionally, the change in degeneracy direction has a significant impact on the apparent constraining power of each of our statistics. For C E C_E C E in particular, we find a significant ∼ 40 % \sim40\% ∼ 40% reduction in the relative marginal uncertainty between our S 8 S_8 S 8 and Σ 8 \Sigma_8 Σ 8 parameters. As such, we opt to preferentially present results in the Σ 8 \Sigma_8 Σ 8 parametrisation in this manuscript, especially when quantifying the impact of analysis choices on recovered constraints. We note, however, that the different α \alpha α values between the different statistics mean that Σ 8 \Sigma_8 Σ 8 values recovered using each statistic are not directly comparable to one another.
Analysis variations
In addition to the fiducial constraints, we provide a range of additional cosmological constraints computed with various permutations of our modelling and optimisation choices. Specifically, we present: the impact of the emulation of the matter power spectrum (Sect. 5.2.1), modelling of the signal with only the Gaussian component of the covariance (Sect. 5.2.2), iterative computation of covariances (which require a pre-defined cosmology; Sect. 5.2.3), the impact of different choices of intrinsic alignment models (Sect. 5.2.4), the role of observational systematics (Sect. 5.2.5), the impact of different redshift calibrations (Sect. 5.2.6), the influence of different scale cuts (Sect. 5.2.7), and the impact of individual tomographic bins on our estimated cosmological parameters (Sect. 5.2.8).
Figure 11 presents a whisker diagram of all the various Σ 8 \Sigma_8 Σ 8 (that is, for the optimal Ω m , σ 8 \Omega_{\rm m}, \sigma_8 Ω m , σ 8 -degeneracy) results presented in this work. The figure includes our fiducial results in each statistic, with constraints computed using the three methods previNotes. α \alpha α in this table refers to the exponent in the equation Σ 8 = σ 8 ( Ω m / 0.3 ) α \Sigma_8 = \sigma_8 (\Omega_{\mathrm{m}}/0.3)^\alpha Σ 8 = σ 8 ( Ω m /0.3 ) α . When α = 0.5 \alpha= 0.5 α = 0.5 (upper section) reported constraints are for S 8 S_8 S 8 . Otherwise, constraints are for the optimal degeneracy direction for each statistic, thereby showing the maximal constraining power. We report the χ 2 \chi^2 χ 2 of the best-fit (i.e. maximum a posteriori) model for each statistic, and the number of effective degrees of freedom (dof, computed as N d a t a − N p a r a m e f f N_{\mathrm{data}} - N_{\mathrm{param}}^{\mathrm{eff}} N data − N param eff where N p a r a m e f f N_{\mathrm{param}}^{\mathrm{eff}} N param eff is estimated as described in Appendix B of Stölzner et al. 2025). We compute the best-fit probability-to-exceed (PTE), which indicates the fraction of possible data vectors, drawn from a generative model described by our best fit, that would be at least as extreme as our observed data vector. We provide three constraints over our primary cosmological parameters per setup: the marginal mode and marginal highest posterior density interval (HPDI), the marginal mean and standard deviation, and the maximum a posteriori point and projected joint highest posterior density interval (PJ-HPD). The latter statistic can be noisy in the limit of few samples N s a m p P J − H P D N_{\mathrm{samp}}^{\mathrm{PJ-HPD}} N samp PJ − HPD residing within the PJ-HPD, and so we also provide the number of samples in the PJ-HPD region as a guide.
!Whisker diagram with two panels, COSEBIs (E n E_n E n ) and Bandpowers (C E C_E C E ), showing Σ 8 = σ 8 ( Ω m / 0.3 ) α \Sigma_8=\sigma_8(\Omega_m/0.3)^\alpha Σ 8 = σ 8 ( Ω m /0.3 ) α cosmological results and analysis variations. Left-panel labels: KiDS–Legacy Fiducial; (Marginal Max + HPDI); (Marginal Mean CI); (Marg. Post. + PJ–HPD); Modelling variations; CAMB + HMCode(2020); Gaussian Covariance; IA model variations; NLA; NLA–z z z ; NLA–k k k ; Alternative redshift calibration; MICE2 δ z \delta_z δ z ; Calib. weight: PV only; Calib weight: shape only; CC δ ( z ) \delta(z) δ ( z ) ; Systematics variation; σ z = 0 \sigma_z=0 σ z = 0 ; σ z = 0 \sigma_z=0 σ z = 0 & σ m = 0 \sigma_m=0 σ m = 0 & σ I A ≃ 0 \sigma_{\mathrm{IA}}\simeq0 σ IA ≃ 0 ; Alternate scale cuts; θ ∈ [ 0.5 , 300 ] \theta\in[0.5,300] θ ∈ [ 0.5 , 300 ] arcmin; Bin combination variations; w/o z z z –bin 1; w/o z z z –bin 2; w/o z z z –bin 3; w/o z z z –bin 4; w/o z z z –bin 5; w/o z z z –bin 6; Planck Legacy. Right-panel labels: Bandpowers (C E C_E C E ). Annotated values: 0.73 σ 0.73\sigma 0.73 σ , − 0.12 σ -0.12\sigma − 0.12 σ , 0.02 σ 0.02\sigma 0.02 σ , − 0.54 σ -0.54\sigma − 0.54 σ , − 0.15 σ -0.15\sigma − 0.15 σ , − 0.54 σ -0.54\sigma − 0.54 σ , 1.25 σ 1.25\sigma 1.25 σ , 0.11 σ 0.11\sigma 0.11 σ , 0.17 σ 0.17\sigma 0.17 σ , − 0.01 σ -0.01\sigma − 0.01 σ , − 0.07 σ -0.07\sigma − 0.07 σ , 0.42 σ 0.42\sigma 0.42 σ , − 0.25 σ -0.25\sigma − 0.25 σ , 0.05 σ 0.05\sigma 0.05 σ , − 0.01 σ -0.01\sigma − 0.01 σ , 0.13 σ 0.13\sigma 0.13 σ , 0.17 σ 0.17\sigma 0.17 σ , 0.22 σ 0.22\sigma 0.22 σ , − 0.32 σ -0.32\sigma − 0.32 σ ; right-panel values 1.01 σ 1.01\sigma 1.01 σ , − 0.08 σ -0.08\sigma − 0.08 σ , 0.03 σ 0.03\sigma 0.03 σ , − 0.14 σ -0.14\sigma − 0.14 σ , 0.13 σ 0.13\sigma 0.13 σ , − 0.17 σ -0.17\sigma − 0.17 σ , 0.99 σ 0.99\sigma 0.99 σ , 0.11 σ 0.11\sigma 0.11 σ , 0.22 σ 0.22\sigma 0.22 σ , − 0.01 σ -0.01\sigma − 0.01 σ , − 0.02 σ -0.02\sigma − 0.02 σ , 0.03 σ 0.03\sigma 0.03 σ , 0.08 σ 0.08\sigma 0.08 σ , 0.00 σ 0.00\sigma 0.00 σ , 0.18 σ 0.18\sigma 0.18 σ , − 0.08 σ -0.08\sigma − 0.08 σ , − 0.22 σ -0.22\sigma − 0.22 σ , 0.11 σ 0.11\sigma 0.11 σ . bbox=0.18,0.08,0.83,0.83
Figure 11. Whisker diagram of Σ 8 \Sigma_8 Σ 8 cosmological results presented in this work. The best-fitting α \alpha α used for each statistic is annotated on the x-axis. Each panel shows the compilation of results for one of our two-point statistics, as annotated. In each panel, we show the fiducial constraints in the relevant statistic with coloured points (colour-coded by the statistic), and highlight the fiducial marginal HPDI with a band of the equivalent colour. We also highlight the marginal HPDI in Σ 8 \Sigma_8 Σ 8 from Planck -Legacy (shaded red). All variations to the analysis are annotated, with the difference relative to the fiducial result quantified using the metric ( Σ 8 f i d − Σ 8 v a r ) / σ Σ 8 v a r (\Sigma_8^{\mathrm{fid}}-\Sigma_8^{\mathrm{var}})/\sigma_{\Sigma_8}^{\mathrm{var}} ( Σ 8 fid − Σ 8 var ) / σ Σ 8 var (listed in blue). The Hellinger tension between the fiducial analyses and Planck is annotated in red.
ously described in Table 4 . All discussions below refer to our fiducial metric, the marginal mode and HPDI summary, unless stated otherwise. Additionally, we have annotated the figure with a measure of the difference between our various results, computed using the metric ( Σ 8 f i d − Σ 8 v a r ) / σ Σ 8 v a r (\Sigma_8^{\mathrm{fid}} - \Sigma_8^{\mathrm{var}})/\sigma_{\Sigma_8}^{\mathrm{var}} ( Σ 8 fid − Σ 8 var ) / σ Σ 8 var . Differences are computed between the Σ 8 \Sigma_8 Σ 8 values of each setup and the fiducial. We also annotate the fiducial constraints with their Euclidean distance (see e.g. Appendix G.1 of Heymans et al. 2021) to the constraints from the public, fiducial Planck-Legacy analysis [103 ] . Finally, we present a similar figure for S 8 S_8 S 8 in Appendix G, and tables of all constraints are presented in Appendix H.
Table 4. Constraints on S 8 S_8 S 8 (α = 0.50 \alpha= 0.50 α = 0.50 ) and Σ 8 \Sigma_8 Σ 8 (α \alpha α free).
Modelling P ( k ) P(k) P ( k ) with camb
We tested the impact of the CosmoPower emulator by re-running our fiducial analysis directly with camb and comparing the cosmological constraints between these runs in Figure 11 . We found that, for all statistics, the camb posteriors were consistent with those estimated with the emulated power spectrum, showing maximal changes in the marginal value of Σ 8 \Sigma_8 Σ 8 of 0.03 σ 0.03\sigma 0.03 σ . We consider this level of change to be consistent with chain-to-chain variations, and therefore negligible.
Gaussian covariance
We tested the impact of our covariance modelling on the recovered value of Σ 8 \Sigma_8 Σ 8 in Figure 11 , by performing the fiducial measurement with only the Gaussian component of the covariance activated in the OneCovariance code. We found a reduction in the uncertainty on Σ 8 \Sigma_8 Σ 8 of between 12% (C E C_E C E ) and 19% (E n E_n E n ), and no change in the value of Σ 8 \Sigma_8 Σ 8 with respect to the fiducial (Δ Σ 8 = 0.02 σ \Delta\Sigma_8 = 0.02\sigma Δ Σ 8 = 0.02 σ ). This is an indication that, while shape noise is still sufficient to shroud information on non-linear scales, KiDS-Legacy has a non-negligible contribution to its constraining power from these non-Gaussian scales.
Iterative covariance
The covariance used in our cosmological inference modelling assumes a fixed cosmology. If this fixed cosmology differs from the best-fit cosmology preferred by the dataset, this can add a source of inconsistency in our fiducial results. Therefore, as with previous KiDS analyses (see e.g. van Uitert et al. [136 ] ), we verified the fidelity of our cosmological inference by recomputing the fiducial sample covariances at the maximum a posteriori cosmology for each chain, and rerun our cosmological inference. This iterative evaluation of the covariance and cosmological parameters do not need to converge after a single iteration, and so we ran this process for many (≲ 10 \lesssim10 ≲ 10 ) iterations. We found that the cosmological constraints (and constraining power) converged after two iterations, even for our most complex model choices. This is in good agreement with previous investigations into the impact of varying cosmological parameters on covariances in the context of stage-III surveys ([106 , 74 ] ).
Intrinsic alignment modelling
We tested the cosmological parameter constraints recovered when modelling intrinsic alignments with our four different models. The resulting cosmological parameters are presented in Figure 11 .
First, we see that the impact of changing our IA models on the marginal constraints of S 8 S_8 S 8 and Σ 8 \Sigma_8 Σ 8 are mild. For our C E C_E C E statistic, we found maximal differences of 0.17 σ 0.17\sigma 0.17 σ across all model variations. Also, E n E_n E n values are more sensitive to our chosen IA model, showing maximal differences of 0.54 σ 0.54\sigma 0.54 σ (for the NLA and NLA-k k k models and Σ 8 \Sigma_8 Σ 8 ). Second, we see that changing the IA model has a minor impact on the constraining power of our survey, with systematic changes in constraining power being most evident for the NLA-z z z model: we see a consistent reduction in
constraining power (compared to the fiducial NLA-M M M model) of 12 − 13% for Σ 8 \Sigma_8 Σ 8 across all statistics. For our NLA and NLA-k k k models, the reduction in constraining power compared to the fiducial is slightly smaller, at ∼ 5 % \sim5\% ∼ 5% . This reduced constraining power (despite the fact that e.g. the NLA model has fewer parameters) is attributed to the informative prior applied in our NLA-M M M model (see Sect. 2.4.4 and Appendix B). Finally, we also note mild changes to our recovered goodness-of-fit under the various IA models. In all cases, the χ 2 \chi^2 χ 2 of the modified IA model is consistent with the fiducial (changes in PTE at the order of 1 − 5%).
The stability of our cosmological constraints under variation of the IA modelling provides confidence that any possible misspecification of the IA model is unlikely to introduce bias in our cosmological parameter constraints. Furthermore, we can compare the implied amplitude of IA under each of our models by examining the A I A , t o t a l A_{\mathrm{IA,total}} A IA , total values (see Eqs. (8) , (10) , (12) ). Figure 12 shows the value of A I A , t o t a l A_{\mathrm{IA,total}} A IA , total inferred from our E n E_n E n analysis, as a function of redshift. All statistics show quantitatively similar behaviour. In Appendix G, we present the posterior distributions in all parameters (Figs. G.2 and G.3) for each IA model.
Figure 12. Constraints on the total intrinsic alignment amplitude from the various intrinsic alignment models, per tomographic bin (where appropriate), computed using the E n E_n E n posteriors. Bands and error bars represent the 1 σ 1\sigma 1 σ posterior constraint on each model. Our NLA-M M M model (red) is shown per tomographic bin, with the point positioned at the mean photo-z z z of the bin.
For our NLA and NLA-k k k models, there is no explicit dependence on tomographic bin number (or, more accurately, redshift), and thus the IA amplitude is necessarily constant across the redshift axis for these models. Both of these models are consistent with IA amplitudes of zero. A change in IA amplitude with redshift is nominally expected, though, as the combination of Malmquist bias and the mass-scaling of IA amplitudes combine to increase the effective measured IA with increasing redshift, thereby motivating our NLA-z z z and NLA-M M M models.
For our fiducial NLA-M M M model, the IA do change with tomographic bin, due to their varying fraction of passive/red galaxies and their typical mass (blue galaxies have zero IA in this model). In our analysis, the alignment amplitudes initially climb through bins one to five, but rapidly drop in bin six where there are very few galaxies with an early-type SED (see Appendix B and Table B.1).9 As such, the NLA-M M M model represents a functional form that is not accessible to any of our other models. It is therefore interesting to note the similarity between the inferred values of A I A , t o t a l A_{\mathrm{IA,total}} A IA , total in our NLA-z z z and NLA-M M M models. Both show the expected increase in alignment amplitudes with increasing redshift, although the NLA-z z z model is at slightly lower amplitudes overall. We attribute this difference to the sixth tomographic bin: as the NLA-z z z model requires that all bins follow the same (monotonic) trend, any reduction in A I A , t o t a l A_{\mathrm{IA,total}} A IA , total between the fifth and sixth tomographic bins will require a compensatory reduction in alignment amplitude across all bins.
Impact of observational and systematic nuisance parameters
To demonstrate the impact of the various nuisance parameters and the relative importance of marginalisation over various priors, we ran a posterior analysis with our intrinsic alignment model uncertainty (σ I A \sigma_{\mathrm{IA}} σ IA ) and redshift and shear calibration uncertainty (σ z , σ m \sigma_z,\sigma_m σ z , σ m ) nuisance parameters set to zero. The associated parameters to these uncertainties were fixed to the centre of their fiducial prior. We repeated the posterior calculation using only the remaining (cosmological and physical) parameters. Figure 11 shows these results. We found that the posterior constraints on Σ 8 \Sigma_8 Σ 8 in the absence of marginalisation over these observational and systematic nuisance parameters are tighter by 10% (20%) for our E n E_n E n (C E C_E C E ) measurements, and are shifted with respect to the fiducial by 0.28 σ 0.28\sigma 0.28 σ (0.14 σ 0.14\sigma 0.14 σ ). This suggests that KiDS-Legacy is statistics-limited.
Alternative redshift distribution calibration
We analysed the changes in parameter constraints found under different redshift calibration scenarios. In particular, we demonstrate the constraints found when we switch off shape-measurement weighting of the calibration sample, and when we switch off prior volume weighting of the calibration sample. In both cases, we found that the cosmological constraints are consistent with our fiducial constraints, showing a mild increase of Δ Σ 8 ≤ 0.17 σ \Delta\Sigma_8 \le0.17\sigma Δ Σ 8 ≤ 0.17 σ . This is, to some degree, expected: if our redshift distribution estimation changed appreciably when modifying our calibration method, this would also be reflected in our estimated bias parameters δ z \delta_z δ z . As such, this test is particularly relevant for demonstrating possible sensitivity to source galaxies that reside outside the redshift baseline of our SKiLLS simulations. The lack of significant change in our cosmological parameters under changing the redshift calibration methodology suggests that our analysis is not strongly influenced by sources outside the redshift limits of SKiLLS.
Conversely, we also show in Figure 11 the result of calibrating our redshift distributions using simulations from MICE2. These are the simulations that were used for KiDS-1000 N ( z ) N(z) N ( z ) calibration, but differ here in that we applied our KiDS-Legacy calibration methodology to these simulations. Wright et al. [138 ] discuss the accuracy of the redshift distribution bias parameters estimated from MICE2, in particular possible biases due to the MICE2 redshift limits of 0.1 ≲ z ≲ 1.4 0.1 \lesssim z \lesssim1.4 0.1 ≲ z ≲ 1.4 . As our bin six N ( z ) N(z) N ( z ) (in particular) extends beyond this upper limit, it is expected that our redshift distribution bias parameters there will be inaccurate. Therefore, we include them here predominantly as a demonstration of the effect of switching from MICE2 to SKiLLS. The maximum difference in the mean redshifts between our fiducial (effective) N ( z ) N(z) N ( z ) and that from MICE2 is Δ μ z , e f f = 0.016 \Delta\mu_{z,\mathrm{eff}} = 0.016 Δ μ z , eff = 0.016 , in bin 3. However, the difference found is coherent between bins 3–6, which is therefore expected to cause a shift in S 8 S_8 S 8 . As seen in Figure 11 , our analysis using N ( z ) N(z) N ( z ) calibrated with the MICE2 simulations increases the central value of our Σ 8 \Sigma_8 Σ 8 constraints, by up to 1.25 σ 1.25\sigma 1.25 σ (for E n E_n E n ). We consider these results inferior to the fiducial redshift calibration, though, given the considerably less realistic simulations (compared to SKiLLS), and the smaller redshift baseline that truncates the sixth bin to half of its supposed width.
In addition to our SOM calibrated redshift distributions, we also show results using our redshift distributions calibrated to cross-correlation measurements (Wright et al. [138 ] ). Our cosmological constraints from this analysis are essentially unchanged in their central value for both E n E_n E n and C E C_E C E , showing a 0.01 σ 0.01\sigma 0.01 σ shift in the central value of Σ 8 \Sigma_8 Σ 8 (computed with the fiducial α \alpha α for each statistic). However, the results have lower constraining power, given the larger uncertainty in (particularly) the sixth tomographic bin: uncertainties in the optimal projection of Σ 8 \Sigma_8 Σ 8 are 23 % 23\% 23% (29 % 29\% 29% ) larger for our CC N ( z ) / δ ( z ) N(z)/\delta(z) N ( z ) / δ ( z ) marginalisation than in the fiducial case for E n E_n E n (C E C_E C E ). Nonetheless, the strong agreement between the estimates made with our CC and SOM N ( z ) N(z) N ( z ) provides confidence that redshift distribution biases in KiDS-Legacy are under control.
Scale cuts
In KiDS-1000, Asgari et al. [8 ] presented cosmological constraints using a range of angular scales from 0.5 ′ ≤ θ ≤ 300 ′ 0.5' \leq\theta\leq300' 0. 5 ′ ≤ θ ≤ 30 0 ′ . Subsequent investigation by Dark Energy Survey and Kilo-Degree Survey Collaboration et al. [28 ] , however, found that an optimal treatment of baryonic systematic effect mitigation (for E n / B n E_n/B_n E n / B n ) was achieved using 2.0 ′ ≤ θ ≤ 300 ′ 2.0' \leq\theta\leq300' 2. 0 ′ ≤ θ ≤ 30 0 ′ . For KiDS-1000 (Asgari et al. [8 ] ), this resulted in a shift in the marginal parameter constraints of S 8 S_8 S 8 at the level of 0.8 σ 0.8\sigma 0.8 σ . Reanalysis of KiDS-1000 by Li et al. (2023a), with improved shape measurement and calibration from Li et al. (2023b), showed a reduced sensitivity to the choice of scale cut, seeing a shift of 0.3 σ 0.3\sigma 0.3 σ . We also tested the influence of changing our scale cuts, shown in the marginal projection of Ω m , S 8 \Omega_m,S_8 Ω m , S 8 in Figure 13 , finding that including smaller scales results in a slight (0.1 σ 0.1\sigma 0.1 σ ) decrease in the marginal value of Σ 8 \Sigma_8 Σ 8 and S 8 S_8 S 8 . We also saw a slight change in the preferred value of our baryonic feedback parameter (see Figure 13 ) where there is a reduction in posterior probability mass at negligible feedback amplitudes (i.e. log 10 ( T A G N / K ) ≈ 7.3 \log_{10}(T_{\mathrm{AGN}}/\mathrm{K}) \approx7.3 log 10 ( T AGN / K ) ≈ 7.3 , consistent with dark-matter only) when adding in smaller scale information. However, this must be interpreted with caution, as both posteriors are formally prior dominated in this parameter. This extra data is expected to be sensitive to baryon feedback, but it is difficult to distinguish from a simple change in the noise realisation. In either case, though, it is clear that the additional scales have at most a minor impact on our cosmological constraints.
Figure 13. Comparison of the marginal constraints on Ω m \Omega_m Ω m and S 8 S_8 S 8 (top) and log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}}/\mathrm{K}) log 10 ( T AGN / K ) (bottom) when using a wider range of scales (0.5 ′ ≤ θ ≤ 300 ′ 0.5' \leq\theta\leq300' 0. 5 ′ ≤ θ ≤ 30 0 ′ , orange), compared to the fiducial (blue-green). We compare the constraints to the full fiducial Planck -Legacy posterior (solid red, left) and the marginal constraint on log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}/\mathrm{K})} log 10 ( T AGN / K ) when analysing our fiducial data vector with cosmology fixed to the Planck-Legacy best-fit and marginalisation over nuisance parameters only (dashed red, right).
We note that the expanded scale cut has a non-negligible effect on the significance of our B B B mode null test, causing it to once again fail our p > 0.01 p > 0.01 p > 0.01 requirement. A full reanalysis of the additional astrometric selections required to reach non-significance in our null test with this expanded scale cut has not been explored here, and so the accuracy of the constraints with these expanded scales is not guaranteed. Nonetheless, the stability of
the cosmological constraints in the presence of both baryonic feedback contamination and possible systematic contamination is encouraging.
Removing tomographic redshift bins
In previous analyses of KiDS, the removal of individual tomographic bins (including all cross-correlations) introduced changes in marginal cosmological constraints of up to 1.8 σ 1.8\sigma 1.8 σ in Σ 8 \Sigma_8 Σ 8 (see e.g. Asgari et al. [8 ] ). We performed the same test here, and show the results in Figure 11 . We found that our cosmological results are more robust to the removal of individual tomographic bins than previous work from KiDS: we observe a maximal difference in Σ 8 \Sigma_8 Σ 8 when removing individual bins of 0.24 σ 0.24\sigma 0.24 σ for E n E_n E n , 0.17 σ 0.17\sigma 0.17 σ for C E C_E C E . We attribute this decreased sensitivity at least partly to the reduction in relative importance of each tomographic bin (as signal is shared more evenly between the bins using our approximately equal-number tomography, and sixth bin). Furthermore, in our companion paper [128 ] , we performed a series of internal consistency tests aimed at specifically estimating the constraints from different combinations of tomographic bins, finding that they were all consistent with each other, unlike KiDS-1000. As such, this represents a considerable improvement in internal consistency for KiDS that we attribute to the enhanced redshift calibration samples and methodology, as well as with the general improvements in data quality seen in Legacy compared to previous KiDS data releases.
Discussion
The results from KiDS-Legacy are consistent with the results of Planck Legacy, with a marginal Σ 8 \Sigma_8 Σ 8 consistency of 0.73 σ 0.73\sigma 0.73 σ for E n E_n E n and 1.01 σ 1.01\sigma 1.01 σ for C E C_E C E . The constraints from KiDS-Legacy are also of comparable constraining power to those from Planck when projected into Σ 8 \Sigma_8 Σ 8 , with Planck having a marginal uncertainty of ± 0.015 \pm0.015 ± 0.015 compared to Legacy’s 0.016 0.016 0.016 .
Compared to previous cosmological constraints from KiDS, Legacy represents a decrease in the marginal tension with Planck (computed using the Hellinger distance) of 2.2 σ 2.2\sigma 2.2 σ , compared to previous work from Li et al. (2023a), who utilised methods and modelling choices that most resemble those used in Legacy. Such a difference in marginal constraints from individual KiDS analyses was deemed acceptable prior to unblinding (indeed, such a judgement call was required prior to unblinding), however was nonetheless surprising to the authors. This prompted an investigation into the underlying drivers of the change in cosmological constraints between the two samples.
A complete description of the investigation into the consistency between KiDS-1000 and KiDS-Legacy is provided in Appendix I. In summary, we found that improvements to our redshift distribution estimation and calibration methodology are responsible for approximately two thirds (1.32 σ 1.32\sigma 1.32 σ ) of the reduction in marginal tension with Planck , as reported by KiDS-1000 and Legacy. An additional reduction is attributable to a reduction in statistical noise from the sources outside the KiDS-1000 footprint. These results are shown graphically in Figure 1 .
We also explored the consistency between our new Legacy cosmological parameters with recent results from the literature. Figure 14 presents the E n E_n E n cosmological constraints from Legacy in the context of two recent stage-III cosmic shear analyses: DES-Y3 hybrid-pipeline results from Dark Energy Survey and Kilo-Degree Survey Collaboration et al. (2023), and 2PCF measurements from HSC-Y3 (Li et al. 2023c). DES-Y3 is a largely independent analysis (there is a small on-sky overlap), which is also most consistent with Legacy: differences between marginal cosmological parameter constraints are 0.75 σ 0.75\sigma 0.75 σ in Σ 8 \Sigma_8 Σ 8 and 0.41 σ 0.41\sigma 0.41 σ in S 8 S_8 S 8 . Comparisons to the joint analysis of DESY3 and KiDS-1000, or to HSC, are both more covariate (due to area overlap) and show larger differences (1.5 σ 1.5\sigma 1.5 σ and 1.8 σ 1.8\sigma 1.8 σ in Σ 8 \Sigma_8 Σ 8 respectively), but are otherwise consistent. Overall, we found good consistency between all samples presented here, and with Planck . To this end, there is no evidence of significant tension in S 8 S_8 S 8 (or in the more constraining Σ 8 \Sigma_8 Σ 8 ) from KiDS-Legacy. Furthermore, the consistency that we found between KiDS-Legacy and other low-redshift cosmological probes motivates the joint analysis of cosmological parameters with these surveys, which can be found in our companion paper [128 ] . Lastly, as mentioned before, we do not include ξ ± \xi_{\pm} ξ ± into our fiducial analysis. However, the constraints and resulting consistency metrics are very similar and is discussed in Appendix F.
Figure 14. Comparison between KiDS-Legacy E n E_n E n cosmological constraints and recent stage-III cosmic shear survey results from the literature. We show the 2D posteriors for each analysis Ω m \Omega_{\mathrm{m}} Ω m , S 8 S_8 S 8 (left) as well as the marginal projections in Σ 8 \Sigma_8 Σ 8 (α = 0.58 \alpha= 0.58 α = 0.58 , middle), and S 8 S_8 S 8 (right). The 1D marginal ‘whisker’ plots are annotated as previously: circular points represent marginal mode and HPD intervals, while diamonds represent maximum a posteriori points and their associated PJ-HPD. Values are all calculated directly from the chains released by each analysis/survey, and in all cases except Legacy we take the MAP directly from the chain. We see that all surveys presented here are consistent with both KiDS-Legacy (black annotations) and Planck (red annotations) using our fiducial consistency metric of ∣ Δ X ∣ ≤ 2.6 σ |\Delta X| \leq2.6\sigma ∣Δ X ∣ ≤ 2.6 σ (where X X X is either Σ 8 \Sigma_8 Σ 8 or S 8 S_8 S 8 ), as estimated using the Hellinger tension between marginal mode summaries.
We performed an analysis of the effective constraining power provided by KiDS-Legacy over the matter power spectrum P ( k ) P(k) P ( k ) and the amplitude of baryon feedback. Figure 15 shows our constraint over P ( k ) P(k) P ( k ) from the fiducial E n E_n E n , estimated by sampling the full posterior of KiDS-Legacy and generating posterior predictive distributions for P ( k ) P(k) P ( k ) using CAMB and HMcode2020 . We show the P ( k ) P(k) P ( k ) constraint normalised to the mean of all posterior samples, with a 0.1 σ 0.1\sigma 0.1 σ region shaded around this central reference. This allows us to directly compare our P ( k ) P(k) P ( k ) constraint to the distribution of P ( k ) P(k) P ( k ) returned by two power spectrum emulators, BACCO [6 ] and EuclidEmulator2 [34 ] , also evaluated at the same posterior samples. We note, though, that these emulators have a reduced range of validity compared to our fiducial priors, particularly in n s n_s n s . Nonetheless, we see that the differences between the emulated power spectra and those returned by CAMB are within the 0.1 σ 0.1\sigma 0.1 σ credible interval of our posterior constraint on P ( k ) P(k) P ( k ) , demonstrating that differences between the various P ( k ) P(k) P ( k ) models (evaluated at the same cosmological and nuisance parameters) are consistently less than 10% of our constraining power, and are thus unlikely to be a significant source of bias in our analysis. With increased constraining power (from e.g. stage-IV surveys), however, the uncertainty in our modelling of P ( k ) P(k) P ( k ) will become a significant fraction of the final uncertainty.
Figure 15. Upper panel: Estimate of the systematic uncertainty in our modelling of the matter power spectrum using camb and HMCode2020, compared to two P ( k ) P(k) P ( k ) emulators BACCO (green) and EuclidEmulator2 (orange). The figure shows the ratio between the power spectrum from camb and HMCode2020 to those returned by each emulator, evaluated at 10 000 representative samples from the fiducial E n / B n E_n/B_n E n / B n posterior with a tenth of a standard-deviation. Lower panel: Constraints on baryonic feedback suppression from KiDS-Legacy. The marginal HPDI of feedback amplitudes from our fiducial E n / B n E_n/B_n E n / B n is shown in red, relative to a range of hydrodynamical models of baryonic feedback suppression from Flamingo Schaye et al. [113 ] .
We also explored the importance of baryon feedback modelling on our cosmological constraints, by examining the posterior on our feedback amplitude parameter log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}}/K) log 10 ( T AGN / K ) in the context of our P ( k ) P(k) P ( k ) constraints. We found that the contribution to the power spectrum from our inferred feedback amplitude is negligible for scales k ≲ 0.2 h / M p c k \lesssim0.2\,h/\mathrm{Mpc} k ≲ 0.2 h / Mpc . Above this scale our posterior constraint on the feedback model introduces a 1 − 10 % 1-10\% 1 − 10% suppression in power (compared to the dark-matter only case), and is consistent with the fiducial baryon feedback model from FLAMINGO at all scales. Above approximately k = 1 h / M p c k = 1\,h/\mathrm{Mpc} k = 1 h / Mpc , our marginal posterior on P ( k ) P(k) P ( k ) excludes the more extreme FLAMINGO feedback models. Our MAP estimate of log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}/K}) log 10 ( T AGN / K ) is at the lower boundary limit of the prior, indicating that our best-fitting model prefers negligible baryon feedback amplitudes. Our marginal constraint over log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}/K}) log 10 ( T AGN / K ) is unconstrained (see Figure 13 ), but is nonetheless skewed towards the MAP value.
Finally, given the agreement between our cosmological constraints and those from Planck , we investigated the inferred amplitude of baryon feedback with cosmological parameters fixed to those reported by Planck . This has the effect of placing a strong informative prior that the cosmological model is as described by Planck , and used cosmic shear to marginalise effectively over only observational systematic effects (IA, δ z \delta z δ z , and baryon feedback). Our posterior constraints over all systematic parameters are essentially unchanged with respect to the fiducial, with the only small change of note being that our baryon feedback posterior skewed slightly more strongly towards the lower prior boundary of log 10 ( T A G N / K ) = 7.3 \log_{10}(T_{\mathrm{AGN}}/K) = 7.3 log 10 ( T AGN / K ) = 7.3 , while the MAP remains at the boundary (see Figure 13 ). In the fiducial case our marginal posterior on log 10 ( T A G N / K ) \log_{10}(T_{\mathrm{AGN}}/K) log 10 ( T AGN / K ) was unconstrained (using the metric of Asgari et al. 2021, see their Appendix A), however with the informative Planck prior we were able to place an upper limit constraint of log 10 ( T A G N / K ) ≤ 7.71 \log_{10}(T_{\mathrm{AGN}}/K) \leq7.71 log 10 ( T AGN / K ) ≤ 7.71 (8.01) at one-sided 68% (95%) confidence. As such, we argue that there is little evidence in KiDS-Legacy for significant amplitudes of baryon feedback, even when placing strong informative priors on cosmological parameters.
Summary
We present a cosmological analysis of the completed Kilo-Degree Survey (KiDS). Utilising data from KiDS Data Release 5 [140 ] , we constructed our cosmic shear source sample from 1347 square degrees of optical and near-infrared imaging, containing deeper imaging, new astrometric and photometric calibrations, and cleaner masking of spurious sources and artefacts. Our improved data was then coupled with a greatly expanded spectroscopic calibration dataset named ‘KiDZ’ (consisting of 2.4 – 5.0 times the spectra as used in the analysis of KiDS-1000), along with an updated redshift distribution calibration methodology (‘gold-weight’; see Wright et al. [138 ] ). This allowed us to calibrate sources at higher redshifts than those in any previous KiDS analysis. We defined six tomographic bins, spanning photometric redshifts between 0.1 < z B ≤ 2.0 0.1 < z_{\mathrm{B}} \leq2.0 0.1 < z B ≤ 2.0 , and calibrated the source redshift distributions with a precision up to six times higher than previously possible. We calibrated the source shapes using multi-band image simulations (Li et al. 2023b), finding that the KiDS-Legacy sample has the lowest systematic contamination (by e.g. PSF leakage) of any KiDS dataset to date. As per standard practices in the field, our cosmological analysis was performed fully blind.
Our fiducial cosmological parameter constraints were estimated with COSEBIs (E n E_n E n ), utilising conservative θ ∈ [ 2.0 , 300.0 ] \theta\in[2.0, 300.0] θ ∈ [ 2.0 , 300.0 ] arcmin scales and a new physically-motivated intrinsic alignment (IA) model. Our fiducial E n E_n E n measurement constrains S 8 = σ 8 Ω m / 0.3 = 0.81 5 − 0.021 + 0.016 S_8 = \sigma_8 \sqrt{\Omega_{\mathrm{m}}/0.3} = 0.815^{+0.016}_{-0.021} S 8 = σ 8 Ω m /0.3 = 0.81 5 − 0.021 + 0.016 , which is in full agreement (Hellinger distance of 0.76 σ 0.76\sigma 0.76 σ ) with the marginal prediction S 8 = 0.83 4 − 0.016 + 0.016 S_8 = 0.834^{+0.016}_{-0.016} S 8 = 0.83 4 − 0.016 + 0.016 from the Planck Legacy cosmic microwave background experiment. Our fiducial measurement represents a ∼ 15 % \sim15\% ∼ 15% improvement in constraining power over the previous results from KiDS. However, we also found a significant change in optimal degeneracy between σ 8 \sigma_8 σ 8 and Ω m \Omega_{\mathrm{m}} Ω m , with Σ 8 ≡ σ 8 ( Ω m / 0.3 ) α \Sigma_8 \equiv\sigma_8(\Omega_{\mathrm{m}}/0.3)^\alpha Σ 8 ≡ σ 8 ( Ω m /0.3 ) α and α = 0.58 \alpha= 0.58 α = 0.58 being preferred by our E n E_n E n measurement. We found Σ 8 = 0.82 1 − 0.016 + 0.014 \Sigma_8 = 0.821^{+0.014}_{-0.016} Σ 8 = 0.82 1 − 0.016 + 0.014 , which is a ∼ 32 % \sim32\% ∼ 32% improvement in optimal constraining power compared to previous KiDS analyses. In this projection, the corresponding Hellinger distance to Planck is 0.73 σ 0.73\sigma 0.73 σ . We further computed the cosmological parameters using band power spectra (C ℓ C_\ell C ℓ ) and found full agreement with our fiducial COSEBIs. We performed a number of analysis variations, with our cosmological parameter constraints shown to be consistent across all modelling choices, statistics, tomographic bins, and summary statistics. Our companion paper [128 ] presents a full suite of internal and external consistency tests, finding the KiDS-Legacy dataset to be the most internally robust sample produced by KiDS to date.
We tested for systematic bias with respect to previous results from KiDS, finding that most of the differences in parameter constraints can be attributed to improved redshift distribution estimation and calibration, as well as new survey area and improved image reduction. We leveraged our consistency with Planck to estimate the preferred baryon feedback amplitude of KiDS-Legacy, allowing us to place an upper limit on log 10 T A G N ≤ 7.71 \log_{10} T_{\mathrm{AGN}} \leq7.71 log 10 T AGN ≤ 7.71 (8.01) at 68% (95%) confidence. As such, the results from KiDS-Legacy paint a harmonious picture of cosmology at the end of stage-III, with full agreement between the cosmological parameters estimated by both low- and high-redshift probes of cosmic structure.
Acknowledgements. The authors would like to sincerely thank Matthias Bartelmann, who has acted as the external blinding coordinator for every blinded analysis performed by KiDS in the last 10+ years. AHW and HHi are supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR), made possible by the Bundesministerium für Wirtschaft und Klimaschutz, under project 50QE2305, and acknowledge funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 “Cosmic Interacting Matters - From Source to Signal”. BS, CH, & ZY acknowledge support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. MA is supported by the UK Science and Technology Facilities Council (STFC) under grant number ST/Y002652/1 and the Royal Society under grant numbers RGS/R222226 and ICAR1231094. MB & PJ are supported by the Polish National Science Center through grant no. 2020/38/E/ST9/00395. MB is also supported by grant no. 2020/39/B/ST9/03494. BG acknowledges support from the UKRI Stephen Hawking Fellowship (grant reference EP/Y017137/1). CH also acknowledges support from the UK Science and Technology Facilities Council (STFC) under grant ST/V000594/1. HHi, CM, RR, & AD are supported by an ERC Consolidator Grant (No. 770935). HHi is also supported by a DFG Heisenberg grant (Hi 1495/5-1). HHo & MY acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program with Grant agreement No. 101053992. BJ acknowledges support by the ERC-selected UKRI Frontier Research Grant EP/Y03015X/1 and by STFC Consolidated Grant ST/V000780/1. KK acknowledges support from the Royal Society and Imperial College. SSL is receiving funding from the programme “Netzwerke 2021”, an initiative of the Ministry of Culture and Science of the State of North-Rhine Westphalia. MvWK acknowledges the support by the UKSA and STFC (grant no. ST/X001075/1). PB acknowledges financial support from the Canadian Space Agency (Grant No. 23EXPROSS1) and the Waterloo Centre for Astrophysics. CG is funded by the MICINN project PID2022-141079NB-C32. JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1) SJ acknowledges the Ramón y Cajal Fellowship (RYC2022-036431-I) from the Spanish Ministry of Science and the Dennis Sciama Fellowship at the University of Portsmouth. LL is supported by the Austrian Science Fund (FWF) [ESP 357-N]. AL acknowledges support from the research project grant ‘Understanding the Dynamic Universe’ funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067. CM acknowledges support from the Beecroft Trust, and the Spanish Ministry of Science under the grant number PID2021-128338NB-I00. LM acknowledges support from STFC grant ST/N000919/1. LMo acknowledges the financial contribution from the grant PRIN-MUR 2022 2022PRIN:LY3 “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 ‘Attività scientifiche per la missione Euclid – fase E’. NRN acknowledges financial support from the National Science Foundation of China, Research Fund for Excellent International Scholars (grant n. 12150710511), and from the research grant from China Manned Space Project n. CMS-CSST-2021-A01. LP acknowledges
support from the DLR grant 50QE2002. MR acknowledges financial support from the INAF grant 2022. TT acknowledges funding from the Swiss National Science Foundation under the Ambizione project PZ00P2_193352. YZ acknowledges the studentship from the UK Science and Technology Facilities Council (STFC). Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 179.A-2004, 177.A-3016, 177.A-3017, 177.A-3018, 298.A-5015. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.
Author Contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (AHW, BS), followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products of this manuscript and release. The second group covers those who have either made a significant contribution to the preparation of data products or to the scientific analysis of KiDS-Legacy. {sup}12` Institute for Theoretical Physics, Utrecht University, Princetonplein 5, 3584CC Utrecht, The Netherlands.
{sup}13` Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
{sup}14` Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
{sup}15` School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, NE1 7RU, Newcastle-upon-Tyne, UK
{sup}16` Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Av. Complutense 40, E-28040 Madrid, Spain
{sup}17` Institute of Cosmology & Gravitation, Dennis Sciama Building, University of Portsmouth, Portsmouth, PO1 3FX, United Kingdom
{sup}18` Dipartimento di Fisica e Astronomia “Augusto Righi” - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy
{sup}19` INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy
{sup}20` Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria
{sup}21` The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova University Centre, SE-106 91 Stockholm, Sweden
{sup}22` Imperial Centre for Inference and Cosmology (ICIC), Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2AZ, UK
{sup}23` Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, United Kingdom.
{sup}24` Donostia International Physics Center, Manuel Lardizabal Ibilbidea, 4, 20018 Donostia, Gipuzkoa, Spain.
{sup}25` Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, D-69120 Heidelberg, Germany; Institute for Theoretical Physics, Philosophenweg 16, D-69120 Heidelberg, Germany
{sup}26` Istituto Nazionale di Fisica Nucleare (INFN) - Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy
{sup}27` Department of Physics “E. Pancini” University of Naples Federico II C.U. di Monte Sant’Angelo Via Cintia, 21 ed. 6, 80126 Naples, Italy
{sup}28` INAF - Osservatorio Astronomico di Padova, via dell’Osservatorio 5, 35122 Padova, Italy
{sup}29` Institute for Particle Physics and Astrophysics, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland
Emulating matter power spectra with CosmoPower
The most computationally expensive section of our modelling pipeline is the prediction of the matter power spectra with camb. Fortunately, we can circumvent this time-consuming computation through the use of artificial neural network emulators, which act to provide a direct non-linear mapping between input cosmological parameters and the corresponding output power spectrum. Such an emulator can be created with relative ease using the CosmoPower framework of Spurio Mancini et al. [127 ] . The original CosmoPower emulator was trained to reproduce the linear matter power spectrum, P m , l i n ( k ) P_{\mathrm{m,lin}}(k) P m , lin ( k ) , and the perturbations caused by non-linear evolution, P m , n l ( k ) / P m , l i n ( k ) P_{\mathrm{m,nl}}(k)/P_{\mathrm{m,lin}}(k) P m , nl ( k ) / P m , lin ( k ) , given a set of relevant input parameters. We updated this process to emulate residuals log 10 [ P ( k ) / P r e f ( k ) ] \log_{10}[P(k)/P_{\mathrm{ref}}(k)] log 10 [ P ( k ) / P ref ( k )] between logarithmic power spectra (for both linear and non-linear spectra) and a single reference power spectrum instead. The reference spectrum was chosen to be near the centre of our sampling/prior volume, which is defined by the ranges of the various parameters in Table 3 . These ranges subsequently also define the range of validity of our emulator.
Construction of the emulator requires a large sample of training and testing power-spectra, which were generated by camb, over the seven-dimensional10 volume that we wish to sample. For this purpose, we generated 450 000 training and 50 000 testing spectra, produced by a Latin hypercube sampling over our volume.
The consistency between our emulated spectra and the target spectra returned by camb and HMcode2020 is validated in the publicly available emulator repository on GitHub11 , through computation of the distribution of absolute fractional residuals between the target and estimated spectra for all 50 000 testing samples. We found that our updated method of emulation produces a more uniform error distribution than the original emulation method from Spurio Mancini et al. [127 ] , across the wide range of scales that we emulate (k ∈ [ 10 − 5 , 20 ] k \in[10^{-5}, 20] k ∈ [ 1 0 − 5 , 20 ] h h h /Mpc), and generally reduces the maximum error over all scales. Our emulator exhibits a maximal error below 0.1% for 68% of our test samples, in both the linear and non-linear power spectra.
Galaxy sample properties for intrinsic alignment modelling
For our fiducial galaxy-type and mass-dependent intrinsic alignment (IA) model, we need to measure the average host halo mass M h M_{\mathrm{h}} M h and the red galaxy fraction f r f_{\mathrm{r}} f r in each tomographic bin. These, and relevant intermediate properties, are summarised in Table B.1. We note that, where appropriate, we computed averages as the means over logarithmic quantities as these propagate straightforwardly to the averages of derived quantities in power-law relations. Moreover, all averages were weighted by the shear measurement weight and redshift estimation gold weight.
To determine f r f_{\mathrm{r}} f r , we used the output of the photometric redshift estimation from the template-fitting code BPZ [14 ] . BPZ employs a set of six model templates, ordered in terms of star formation activity, and interpolates linearly between adjacent template to determine a best-fitting spectral energy distribution. The BPZ output T B T_{\mathrm{B}} T B encodes the best-fit templates and their admixture in steps of 0.1. We chose T B < 1.9 T_{\mathrm{B}} < 1.9 T B < 1.9 to define early-type galaxies, encompassing all galaxies with some contribution of an elliptical galaxy spectrum (template no. 1). We validated that, at least on a matched sample between the KiDS and GAMA (Driver et al. [31 ] ) surveys (i.e. for relatively bright and low-redshift galaxies), this T B T_{\mathrm{B}} T B cut is capable of isolating the red sequence.
We found that in the first five tomographic bins f r f_{\mathrm{r}} f r fluctuates around 0.2, broadly in line with expectations for a flux-limited sample. The sixth bin contains a negligible fraction of early-type galaxies, a feature that we also see in our simulated sixth bin from KiLLS. The statistical error on f r f_{\mathrm{r}} f r is very small. We repeated the measurement varying T B T_{\mathrm{B}} T B by ± 0.1 \pm0.1 ± 0.1 and obtained changes in f r f_{\mathrm{r}} f r of less than 0.02. Therefore, we decided to consider f r f_{\mathrm{r}} f r as fixed in the IA model.
To compute halo masses, we made use of the r r r -band luminosities derived as part of the KiDS Data Release 5 (Wright et al. [140 ] ). van Uitert et al. [135 ] measured effective halo masses for SDSS (LOWZ and CMASS) early-type galaxies via weak lensing, as a function of luminosity and redshift. We used their reported halo mass results and galaxy properties to construct the following linear relation
log 10 M h ( i ) [ h − 1 M ⊙ ] = γ ( 1 1 + ⟨ z ⟩ ( i ) − 1 1 + z p i v ) + δ log 10 ( ⟨ L r ⟩ ( i ) L 0 ) + Δ M (28) \log_{10} M_{\mathrm{h}}^{(i)} [h^{-1} M_{\odot}] = \gamma\left( \frac{1}{1+\langle z\rangle^{(i)}} - \frac{1}{1+z_{\mathrm{piv}}} \right) + \delta\log_{10}\left(\frac{\langle L_{\mathrm{r}}\rangle^{(i)}}{L_0}\right) + \Delta M \tag*{(28)} log 10 M h ( i ) [ h − 1 M ⊙ ] = γ ( 1 + ⟨ z ⟩ ( i ) 1 − 1 + z piv 1 ) + δ log 10 ( L 0 ⟨ L r ⟩ ( i ) ) + Δ M (28)
where we chose z p i v = 0.3417 z_{\mathrm{piv}} = 0.3417 z piv = 0.3417 , such that the parameters γ \gamma γ and Δ M \Delta M Δ M are uncorrelated to a good approximation. An additional parameter δ \delta δ describes the dependence on the average luminosity ⟨ L r ⟩ \langle L_{\mathrm{r}}\rangle ⟨ L r ⟩ , normalised by the reference luminosity L 0 = 10 11 ( h / 0.7 ) − 2 L ⊙ L_0 = 10^{11}(h/0.7)^{-2}L_{\odot} L 0 = 1 0 11 ( h /0.7 ) − 2 L ⊙ . van Uitert et al. [135 ] showed that the luminosity dependence evolves only mildly with redshift, so that we determine δ \delta δ as the inverse variance-weighted average scaling over the four samples they considered. Thus, we obtained γ = 1.25 ± 0.33 \gamma= 1.25 \pm0.33 γ = 1.25 ± 0.33 , δ = 1.42 ± 0.16 \delta= 1.42 \pm0.16 δ = 1.42 ± 0.16 , and Δ M = 13.56 ± 0.02 \Delta M = 13.56 \pm0.02 Δ M = 13.56 ± 0.02 . In line with expectations, both the average luminosity and halo mass increase monotonically by more than an order of magnitude from the first to the sixth tomographic bin (see Table B.1).
Statistical uncertainties on the halo mass are propagated via Monte Carlo sampling from the fit and measurement uncertainties. The standard deviation of log 10 M h \log_{10} M_{\mathrm{h}} log 10 M h is approximately 0.2 in the first tomographic bin, driven by the extrapolation of Eq. (B.1) to fainter samples, and of order 0.1 in the remaining bins, where
Table B.1. Early-type galaxy sample properties for the mass-dependent intrinsic alignment model per tomographic bin.
Notes. Columns, from left to right: tomographic bin number, mean redshift, logarithmic normalised galaxy luminosity in the r r r -band, logarithm of the effective halo mass, and the red galaxy fraction. a reduced error due to luminosities more compatible with the range fitted by van Uitert et al. [135 ] is gradually traded off with larger uncertainty because of extrapolation to higher redshifts. Halo mass estimates are strongly correlated (correlation coefficient greater than 0.9) between tomographic bins due to the extrapolation from the luminosity and redshift ranges considered by van Uitert et al. [135 ] . Therefore, we propagated the uncertainty in halo masses into our IA model by assuming a multivariate Gaussian prior with standard deviation of 0.2 in log 10 M h \log_{10} M_{\mathrm{h}} log 10 M h for all bins and the correlation matrix as obtained from the Monte Carlo error propagation.
Throughout the construction of our astrophysics-informed galaxy IA model, we have not distinguished between central and satellite galaxies, although their alignment mechanisms are believed to be different [116 ] . The NLA model assumes a fixed relation between linear and non-linear scales where centrals and satellites tend to dominate IA signals, respectively, although our scale-dependent IA model tests whether the data prefers a different ratio on average over the whole survey. However, the satellite fraction in our early-type galaxy sample is unknown, yet almost certainly evolves with halo mass and redshift (see e.g. Shuntov et al. [124 ] ). For a given mass and redshift, early-type satellite fractions are typically twice as large as those of late types [92 ] , so that the impact of potential IA signals of blue satellites is less of a concern.
Moreover, the satellite fractions in most galaxy samples used in direct IA constraints are also unknown. Fortuna et al. [36 ] estimated a satellite fraction of 30 % and possibly higher (see also van Uitert et al. [134 ] ) in a subset of the GAMA-KiDS sample and subsequently excluded it from their fit. Johnston et al. [68 ] studied alignments amongst and between centrals and satellites. They tentatively found a steeper slope of the red satellite IA signal on small transverse scales. The mass scaling of the IA amplitude in red satellites is qualitatively similar to that of centrals [37 ] .
In summary, neither cosmic shear surveys nor direct measurements to calibrate IA currently provide sufficient information to treat the IA signals of satellites and their associated trends separately. We are forced to treat satellite effects as part of an effective small-scale IA model, with any discrepancies in the combined dependencies on scale and galaxy sample properties expected to be comfortably absorbed in the still wide priors. Once surveys reach the statistical power to significantly constrain IA models beyond an overall amplitude, calibration measurements of IA will need to be (re-)analysed separately for central and satellite populations. Halo models [37 ] will be a natural choice to accurately represent different IA mechanisms and scalings on large and small scales.
Band-powers derivation
To arrive at the correct response functions for our band-powers, we first write C E / B , l C_{\mathrm{E/B},l} C E/B , l as a function of ξ ± ( θ ) \xi_{\pm}(\theta) ξ ± ( θ ) , in the same manner as for the COSBIs in Eq. (16) ,
C E / B , l = π N l ∫ 0 ∞ d θ θ T ( θ ) [ g + l ( θ ) ξ + ( θ ) ± g − l ( θ ) ξ − ( θ ) ] . (29) C_{\mathrm{E/B},l} = \frac{\pi}{N_l} \int_{0}^{\infty} \mathrm{d}\theta\,\theta\,T(\theta)\left[g_{+}^{l}(\theta)\,\xi_{+}(\theta) \pm g_{-}^{l}(\theta)\,\xi_{-}(\theta)\right]. \tag*{(29)} C E/B , l = N l π ∫ 0 ∞ d θ θ T ( θ ) [ g + l ( θ ) ξ + ( θ ) ± g − l ( θ ) ξ − ( θ ) ] . (29)
where g ± l ( θ ) g_{\pm}^{l}(\theta) g ± l ( θ ) are the band-power filter functions for bin l l l and T ( θ ) T(\theta) T ( θ ) is a function that limits the range of θ \theta θ -scales that contribute to the integral. If ξ ± \xi_{\pm} ξ ± is known for all scales then T ( θ ) = 1 T(\theta)=1 T ( θ ) = 1 everywhere. We now write g ± l ( θ ) g_{\pm}^{l}(\theta) g ± l ( θ ) in terms of the idealised response function,
g ± l ( θ ) = ∫ 0 ∞ d l l S l ( l ) J 0 / 4 ( l θ ) . (30) g_{\pm}^{l}(\theta) = \int_{0}^{\infty} \mathrm{d}l\,l\,S_l(l) J_{0/4}(l\theta). \tag*{(30)} g ± l ( θ ) = ∫ 0 ∞ d l l S l ( l ) J 0/4 ( l θ ) . (30)
If S l ( l ) S_l(l) S l ( l ) is a top-hat then,
g + l ( θ ) = 1 θ 2 [ θ l u p , l J 1 ( θ l u p , l ) − θ l l o , l J 1 ( θ l l o , l ) ] , (31) g_{+}^{l}(\theta) = \frac{1}{\theta^2}\left[\theta l_{\mathrm{up},l} J_1(\theta l_{\mathrm{up},l})-\theta l_{\mathrm{lo},l} J_1(\theta l_{\mathrm{lo},l})\right], \tag*{(31)} g + l ( θ ) = θ 2 1 [ θ l up , l J 1 ( θ l up , l ) − θ l lo , l J 1 ( θ l lo , l ) ] , (31)
g − l ( θ ) = 1 θ 2 [ G − ( θ l u p , l ) − G − ( θ l l o , l ) ] , g_{-}^{l}(\theta) = \frac{1}{\theta^2}\left[\mathcal{G}_{-}(\theta l_{\mathrm{up},l})-\mathcal{G}_{-}(\theta l_{\mathrm{lo},l})\right], g − l ( θ ) = θ 2 1 [ G − ( θ l up , l ) − G − ( θ l lo , l ) ] ,
with,
G − ( x ) = ( x − 8 x ) J 1 ( x ) − 8 J 2 ( x ) . (32) \mathcal{G}_{-}(x)=\left(x-\frac{8}{x}\right)J_1(x)-8J_2(x). \tag*{(32)} G − ( x ) = ( x − x 8 ) J 1 ( x ) − 8 J 2 ( x ) . (32)
Inserting ξ ± \xi_{\pm} ξ ± from Eq. (15) into Eq. (C.1) we arrive at,
C E , l ( i j ) = 1 2 N l ∫ 0 ∞ d l l [ W E E l ( l ) C E E , E ( i j ) ( l ) + W E B l ( l ) C E E , B ( i j ) ( l ) ] , (33) C_{\mathrm{E},l}^{(ij)} = \frac{1}{2N_l}\int_{0}^{\infty}\mathrm{d}l\,l\left[W_{\mathrm{EE}}^{l}(l)C_{\mathrm{EE},\mathrm{E}}^{(ij)}(l)+W_{\mathrm{EB}}^{l}(l)C_{\mathrm{EE},\mathrm{B}}^{(ij)}(l)\right], \tag*{(33)} C E , l ( ij ) = 2 N l 1 ∫ 0 ∞ d l l [ W EE l ( l ) C EE , E ( ij ) ( l ) + W EB l ( l ) C EE , B ( ij ) ( l ) ] , (33)
C B , l ( i j ) = 1 2 N l ∫ 0 ∞ d l l [ W B E l ( l ) C E E , E ( i j ) ( l ) + W B B l ( l ) C E E , B ( i j ) ( l ) ] , C_{\mathrm{B},l}^{(ij)} = \frac{1}{2N_l}\int_{0}^{\infty}\mathrm{d}l\,l\left[W_{\mathrm{BE}}^{l}(l)C_{\mathrm{EE},\mathrm{E}}^{(ij)}(l)+W_{\mathrm{BB}}^{l}(l)C_{\mathrm{EE},\mathrm{B}}^{(ij)}(l)\right], C B , l ( ij ) = 2 N l 1 ∫ 0 ∞ d l l [ W BE l ( l ) C EE , E ( ij ) ( l ) + W BB l ( l ) C EE , B ( ij ) ( l ) ] ,
where the weight functions are,
W E E l ( l ) = ∫ 0 ∞ d θ θ T ( θ ) [ J 0 ( l θ ) g + l ( θ ) + J 4 ( l θ ) g − l ( θ ) ] , W_{\mathrm{EE}}^{l}(l)=\int_{0}^{\infty}\mathrm{d}\theta\,\theta\,T(\theta)\left[J_0(l\theta)g_{+}^{l}(\theta)+J_4(l\theta)g_{-}^{l}(\theta)\right], W EE l ( l ) = ∫ 0 ∞ d θ θ T ( θ ) [ J 0 ( l θ ) g + l ( θ ) + J 4 ( l θ ) g − l ( θ ) ] ,
W E B l ( l ) = ∫ 0 ∞ d θ θ T ( θ ) [ J 0 ( l θ ) g + l ( θ ) − J 4 ( l θ ) g − l ( θ ) ] , (34) W_{\mathrm{EB}}^{l}(l)=\int_{0}^{\infty}\mathrm{d}\theta\,\theta\,T(\theta)\left[J_0(l\theta)g_{+}^{l}(\theta)-J_4(l\theta)g_{-}^{l}(\theta)\right], \tag*{(34)} W EB l ( l ) = ∫ 0 ∞ d θ θ T ( θ ) [ J 0 ( l θ ) g + l ( θ ) − J 4 ( l θ ) g − l ( θ ) ] , (34)
with W E E l ( l ) = W B B l ( l ) W_{\mathrm{EE}}^{l}(l)=W_{\mathrm{BB}}^{l}(l) W EE l ( l ) = W BB l ( l ) and W E B l ( l ) = W B E l ( l ) W_{\mathrm{EB}}^{l}(l)=W_{\mathrm{BE}}^{l}(l) W EB l ( l ) = W BE l ( l ) . Since ξ ± \xi_{\pm} ξ ± are measured between θ min \theta_{\min} θ m i n and θ max \theta_{\max} θ m a x , T ( θ ) T(\theta) T ( θ ) has to be defined such that no information beyond this range is included. To avoid ringing, we apodise T ( θ ) T(\theta) T ( θ ) with a Hahn window and smooth its transitions to zero:
T ( θ ) = { 0 ; x < x l o − Δ x 2 cos 2 [ π 2 x − ( x l o + Δ x / 2 ) Δ x ] ; x l o − Δ x 2 ≤ x < x l o + Δ x 2 1 ; x l o + Δ x 2 ≤ x < x u p − Δ x 2 cos 2 [ π 2 x − ( x u p − Δ x / 2 ) Δ x ] ; x u p − Δ x 2 ≤ x < x u p + Δ x 2 0 ; x ≥ x u p + Δ x 2 . (35) T(\theta)=
\begin{cases}
0;&x<x_{\mathrm{lo}}-\dfrac{\Delta_x}{2}\\
\cos^2\left[\dfrac{\pi}{2}\dfrac{x-(x_{\mathrm{lo}}+\Delta_x/2)}{\Delta_x}\right];&x_{\mathrm{lo}}-\dfrac{\Delta_x}{2}\le x<x_{\mathrm{lo}}+\dfrac{\Delta_x}{2}\\
1;&x_{\mathrm{lo}}+\dfrac{\Delta_x}{2}\le x<x_{\mathrm{up}}-\dfrac{\Delta_x}{2}\\
\cos^2\left[\dfrac{\pi}{2}\dfrac{x-(x_{\mathrm{up}}-\Delta_x/2)}{\Delta_x}\right];&x_{\mathrm{up}}-\dfrac{\Delta_x}{2}\le x<x_{\mathrm{up}}+\dfrac{\Delta_x}{2}\\
0;&x\ge x_{\mathrm{up}}+\dfrac{\Delta_x}{2}.
\end{cases} \tag*{(35)} T ( θ ) = ⎩ ⎨ ⎧ 0 ; cos 2 [ 2 π Δ x x − ( x lo + Δ x /2 ) ] ; 1 ; cos 2 [ 2 π Δ x x − ( x up − Δ x /2 ) ] ; 0 ; x < x lo − 2 Δ x x lo − 2 Δ x ≤ x < x lo + 2 Δ x x lo + 2 Δ x ≤ x < x up − 2 Δ x x up − 2 Δ x ≤ x < x up + 2 Δ x x ≥ x up + 2 Δ x . (35)
Here x = log ( θ ) x=\log(\theta) x = log ( θ ) , Δ x \Delta_x Δ x is the log-width of the apodisation, x l o − Δ x / 2 = log ( θ min ) x_{\mathrm{lo}}-\Delta_x/2=\log(\theta_{\min}) x lo − Δ x /2 = log ( θ m i n ) and x u p + Δ x / 2 = log ( θ max ) x_{\mathrm{up}}+\Delta_x/2=\log(\theta_{\max}) x up + Δ x /2 = log ( θ m a x ) . We require Δ x > log ( θ u p / θ l o ) \Delta_x>\log(\theta_{\mathrm{up}}/\theta_{\mathrm{lo}}) Δ x > log ( θ up / θ lo ) , and in practice have used Δ x = 0.5 \Delta_x=0.5 Δ x = 0.5 for all KiDS-Legacy C E , B C_{\mathrm{E,B}} C E , B analyses. We note, also, that this apodisation method is different to how we applied the apodisation for KiDS-1000 analysis, where we instead had x l o = log ( θ min ) x_{\mathrm{lo}}=\log(\theta_{\min}) x lo = log ( θ m i n ) and x u p = log ( θ max ) x_{\mathrm{up}}=\log(\theta_{\max}) x up = log ( θ m a x ) . This change was made as it was deemed more intuitive/appropriate that the band power windows be contained primarily within the requested l l l ranges after apodisation, rather than before.
KiDS-Legacy variable depth mocks and covariance validation
Our analysis used an analytical covariance matrix, as summarised in Sect. 2.6. Here, we describe the mocks which were used to validate the covariance matrix [107 ] in more detail.
Figure D.1. Cartesian spatial map (corresponding N s i d e = 1024 N_{\rm side}=1024 N side = 1024 , showing a local flat projection of the healpix map) of the TG weights w T G w_{\rm TG} w TG throughout a single field of our KiDS-Legacy-like mock catalogue. Larger TG weights correspond to a higher local source density.
Figure D.2. Dependence of the galaxy density, n e f f n_{\rm eff} n eff , on TG weights, w T G w_{\rm TG} w TG in the KiDS-Legacy-like mock catalogue. The data points represent the mean n e f f n_{\rm eff} n eff of ten equi-populated bins in w T G w_{\rm TG} w TG . The solid line shows the linear fit to the aforementioned data points of their respective tomographic bin, labelled S1-S6. The dashed horizontal lines show the mean values of n e f f n_{\rm eff} n eff calculated from the galaxy samples with variable depth per tomographic bin, while the dashed horizontal lines show the values of n e f f n_{\rm eff} n eff for the respective galaxy samples without any spatial variations in the observational depth.
The log-normal simulations are based on GLASS10 ^{10} 10 (Generator for Large-Scale Structure, Tessore et al. 2023) with a detailed description given in von Wietersheim-Kramsta et al. [137 ] as well as Section 8.1 in Reischke et al. [107 ] . The resulting density field realisations, in concentric shells, were subsequently integrated along the line-of-sight weighted by the weak lensing kernel, Eq. (2) , to create a total of 4224 realisations of the shear field.
For each realisation, galaxies were sampled, factoring in systematic effects such as the survey footprint and variable depth (as caused by seeing, sky transparency fluctuations, and differences in the number of overlapping pointings) using SALMO11 ^{11} 11 (Speedy Acquisition for Lensing and Matter Observables, Joachimi et al. 2021). The inclusion of a sixth tomographic bin in KiDS-Legacy incorporates more faint galaxies in the sample, enhancing the overall magnitude of variable depth and boosting the lensing signal relative to the noise. Although previous analyses suggested that depth variability effects on the signal in KiDS-1000 remain below the noise threshold [48 ] , this might not be the case for KiDS-Legacy.
To model variable depth, we rely on organised randoms (ORs, introduced in Johnston et al. 2021; Yan et al. 2025). ORs use SOMs to generate random samples that reflect the spatial variability caused by a high-dimensional systematic space. The ORs provide Tiaogeng (TG) weights, w T G w_{\rm TG} w TG , [see 141 for a detailed explanation] , quantifying the spatial variability of the selection function of the galaxy sample as a function of position on the sky. The ORs robustly trace the variations in depth for KiDS-Legacy. Yan et al. [141 ] demonstrated that the TG weights are uncorrelated with r r r -band magnitude, photometric redshift, and underlying cosmological density variations. In contrast, they are highly correlated with the magnitude limit of the survey in the r r r -band, which is used for both source extraction and shape measurement.
Figure D.1 shows the TG weights for a single field in our KiDS-Legacy-like mock, as a function of position on the simulated sky. Additionally, Fig. D.2 shows the corresponding local source density of sources as a function of TG weights, colour-scaled for the different source bins. The horizontal dashed lines depict the effective homogeneous number density. Solid lines are linear fits to the measured effective number density in ten equi-populated bins of the TG weights. This linear interpolation allows for an effective evaluation of n e f f n_{\rm eff} n eff at each pixel.
The local density variations seen here also affect the redshift distributions. This issue is addressed by reapplying the photometric redshift calibration using SOMs [see 138 for details] to each tomographic bin and each of the ten TG weight bins. The resulting redshift distribution then enters the shear maps via the weak lensing kernel.
On each simulation, we measured the shear two-point correlation functions as indicated in Sect. 2.5, from which the sample covariance could be estimated directly. These mocks were then compared to the ONECOVARIANCE code. The results are presented in Reischke et al. [107 ] with 10 percent agreement for the variances, very similar to what was found in previous analysis [65 ] . In particular, we found that the limiting factor in the covariance modelling is not the effect of variable depth but rather the survey geometry on large scales. All these effects, however, do not significantly affect the cosmological results presented in this paper.
B-mode analysis
During the initial phases of the blinded KiDS-Legacy analysis, our B B B -mode null-test analyses (see Sect. 4.3) failed our required significance threshold (p ≥ 0.01 p \geq0.01 p ≥ 0.01 ). This led to a dedicated effort to investigate the possible causes of the B B B -mode signal in the dataset, starting with the highest level data products and systematically working our way back through the data production and analysis pipelines, to identify the root cause of the B B B -mode signal.
As KiDS-1000 did not show significant B B B modes (specifically, as analysed by Li et al. 2023a using the same shape-measurement code and analysis choices, such as scale cuts), the initial test in our investigation was to determine whether the observed B B B -mode significance in KiDS-Legacy was owing to a decrease in covariance (specifically shot noise) caused by the increased area of KiDS-Legacy over KiDS-1000. We therefore started by analysing the KiDS-Legacy B B B -mode significance in the same footprint as KiDS-1000, finding that the KiDS-Legacy data in this footprint still exceeded the B B B -mode significance requirement.
Subsequent investigation of the B B B -mode origin focussed on PSF contamination (via Paulin-Henriksson statistics; see Sect. 4), one- and two-dimensional c c c -terms, covariance construction, and shape measurement biases (as a function of, for example, galaxy morphology and position on-sky). Each of these tests was performed locally (per square degree tile), in the detector plane (i.e. coadding all sources in the frame of the focal plane), with spatial jackknife subsets of the full survey, and globally (e.g. at the data vector level, and looking for population-level trends). Each of these investigations failed to illuminate the origin of the B B B -mode signal. Finally, we investigated the lowest-level possible cause of the anomalous B B B -mode signal: systematic failures in the astrometric solution used in the theli calibration (with Gaia) of r r r -band images. The global astrometric accuracy of KiDS-Legacy is presented in Wright et al. [140 ] , where we demonstrated the accuracy of the theli (constructed with a Gaia reference) and Astro-WISE (constructed with a SDSS and 2MASS reference) catalogues both globally (see their Figures 10, 11, C1, and C2) and as a function of right ascension and declination (RA and Dec; see their Figures 8 and 9), compared to SDSS, 2MASS, and Gaia. In our new investigation, we explored the astrometric accuracy in finer detail (i.e. higher on-sky resolution) and using the individual astrometric solutions of theli exposures and Astro-WISE co-adds.
We computed two test statistics for investigation of the B B B -mode signal. The first test statistic is sensitive to the astrometric consistency between theli co-added images, which are used for source detection, and Astro-WISE co-added images, which are used for photometric measurements (and therefore impacts photometric redshift estimation, which is used for tomographic binning). This metric (δ r i n t \delta r_{\rm int} δ r int ) measured the systematic difference between the centroid of all targets extracted from these images. This was achieved by matching all sources (within 2 ′ ′ 2^{\prime\prime} 2 ′′ ) extracted from Astro-WISE and theli co-added r r r -band images. With these matches, we computed the radial separation r r r between each source’s centroid in the two images and then calculated the median of these separations in 1 × 1 1 \times1 1 × 1 arcminute on-sky bins. The 2D map of median separations was then assigned back to each source.
Our second test statistic was computed between individual exposures of theli r r r -band images, which are used for shape measurement. This metric (δ r e x p \delta r_{\rm exp} δ r exp ) measures the systematic difference between the centroids of targets extracted from individual theli exposures, which are point-like and have high S/N. This was achieved by matching sources extracted from Astro-WISE (used here purely as a static reference) and the individual theli exposures (within 2 ′ ′ 2^{\prime\prime} 2 ′′ ). For the Astro-WISE static reference, we applied selections on Source Extractor [15 ] output variables
F L U X _ R A D I U S < 4.0 (36) \mathrm{FLUX\_RADIUS} < 4.0 \tag*{(36)} FLUX_RADIUS < 4.0 (36)
and
F L U X _ I S O / F L U X E R R _ I S O > 10.0 , (37) \mathrm{FLUX\_ISO}/\mathrm{FLUXERR\_ISO} > 10.0, \tag*{(37)} FLUX_ISO / FLUXERR_ISO > 10.0 , (37)
to select point-like high-S/N sources. Settings used in the extraction of sources are provided in Table 4 of Wright et al. [140 ] . For the sources extracted from the individual theli exposures, we performed the same point source selection (Eq. (E.1)) and a modified S/N criteria, to account for the reduced depth of the individual exposures compared to the co-add:
F L U X _ I S O / F L U X E R R _ I S O > 10.0 / ( 5 ) . (38) \mathrm{FLUX\_ISO}/\mathrm{FLUXERR\_ISO} > 10.0/\sqrt{(5)} . \tag*{(38)} FLUX_ISO / FLUXERR_ISO > 10.0/ ( 5 ) . (38)
With these matches, we then computed the radial separation r i j r_{ij} r ij between all exposures in { i , j ∈ [ 1 , 5 ] ∣ i ≠ j } \{i,j \in[1,5] \mid i \ne j\} { i , j ∈ [ 1 , 5 ] ∣ i = j } , computed the maximum of these separations in 2 × 2 2 \times2 2 × 2 arcminute on-sky bins, and smoothed these maps with a 1 1 1 arcminute Gaussian kernel. We then combined these individual 15 i ↔ j 15\ i \leftrightarrow j 15 i ↔ j maps into a final single map, by selecting again the maximum separation per on-sky bin. This resulting map was then once-again smoothed using a 1 1 1 arcminute Gaussian kernel.
These maps provided us with two test statistics that could be used as an additional mask for the survey, given specific thresholds in the statistic values. The first metric is relevant for possible failures in photometry: if the astrometric solution between Astro-WISE and theli differs, then the apertures extracted from theli imaging will be incorrectly placed on the Astro-WISE images. The second metric is relevant for possible failures in shape measurement: the shape measurement algorithm lensfit, which is used in KiDS, assumes perfect astrometric alignment of exposures within the theli r r r -band imaging. As such, any residual astrometric misalignment will manifest as an erroneously estimated shape, and coherent misalignments may manifest as an E E E - or B B B -mode cosmic-shear signal.
For our first statistic, the degree of inconsistency between the imaging that can lead to systematic errors in photometry is relatively benign: our photometric measurement code GAAP has a minimum aperture size of 0\farcs7 , driven largely by the fact that our imaging has seeing that is consistently larger than 0\farcs6 (see e.g. Figure 7 of Wright et al. [140 ] ). Systematic errors in aperture centroids at the level of the FWHM and larger will lead to significant flux to be lost from the target aperture (assuming a point-like source), and thus may pose a problem for our photometry. However, there are two additional relevant considerations. Firstly, in our cosmological analyses, all sources that are used must be resolved and thus (in some direction) larger than the PSF FWHM. Secondly, even if significant flux is missed, if the Astro-WISE imaging are all well aligned to each other, then the flux will be consistently missed in all bands, thereby keeping the source SED largely consistent. As such, we expect the effect to be rarely relevant, and (when it is relevant) only to have a minor impact on photometric redshifts. Nonetheless, we chose to mask all areas of the survey where our first statistic exceeds the minimum PSF FWHM: that is, where \delta r_{\rm int} > 0\farcs6 . This corresponds to approximately 0.7 d e g 2 0.7\ {\rm deg}^2 0.7 deg 2 of removed data over the full KiDS-Legacy footprint.
For our astrometric statistic, there is much less intuition available regarding the level of astrometric error that may be problematic for our analysis. As such, we chose to investigate (while remaining fully blinded; outcomes of the various tests were very similar for all blinds) the impact of various masks (defined by thresholds in δ r e x p \delta r_{\rm exp} δ r exp ) on our null test results, and on our cosmological inference. Figure E.1 shows the results of this investigation. On the x x x -axis of each panel in the figure, we show the threshold in our astrometric metric δ r e x p \delta r_{\rm exp} δ r exp that was used to define a mask. This mask was then applied to the full KiDS-Legacy cosmic shear source sample, and the end-to-end cosmological analysis pipeline was run (including N ( z ) N(z) N ( z ) estimates, shape recalibration, two-point statistics, covariances, and cosmological chains). For each threshold, we were therefore able to report the change in sample size (bottom panel of Fig. E.1; similar to the change in area), the change in the B B B -mode null test p p p -value (middle panel), and the change in the inferred marginal value of Σ 8 \Sigma_8 Σ 8 (top panel).
Firstly, we see that the null test p p p -value passes our predefined threshold of p = 0.01 p = 0.01 p = 0.01 with \delta r_{\rm exp} \lesssim0\farcs25 , corresponding to a masked fraction of sources of at least 1 % 1\% 1% . Secondly, we see that the value of Σ 8 \Sigma_8 Σ 8 is robust to the application of increasingly conservative thresholds, demonstrating a Δ Σ 8 ≈ 0.5 σ \Delta\Sigma_8 \approx0.5\sigma Δ Σ 8 ≈ 0.5 σ variation when discarding up to 100 d e g 2 \mathrm{deg}^2 deg 2 of possibly contaminated data and causing the B B B -mode significance to decrease to p > 0.3 p > 0.3 p > 0.3 . In the region around where we first pass the null test, and do not throw away an overly excessive fraction of the source sample (i.e. \delta_{\rm rexp} \in[0\farcs20, 0\farcs25] ) our recovered marginal value of Σ 8 \Sigma_8 Σ 8 is at most 0.006 (or 0.36 σ 0.36\sigma 0.36 σ ) smaller than in the case where we perform no masking at all (equivalent to the threshold of 1 ′ ′ 1'' 1 ′′ ). Thirdly, we can see that the null test p p p -value is correlated with the fraction of removed area, which begs the question of whether the improvement in our null test might simply be the result of reduced constraining power related to the smaller sample. To test this possibility, we must construct a new mask that removes the same amount of area from the survey at a given threshold, but which does not trace the underlying astrometric errors that we think are likely responsible for our B B B -mode signal. To do this, we constructed a randomised version of our astrometric metric, by randomising the individual masks (defined per square degree tile) between the 1347 tiles of the survey.
Figure E.1. Summary of our blinded B B B -mode investigation where we measure the COSEBIs E n / B n E_n/B_n E n / B n cosmic shear statistics, masking survey area where our astrometric metric, δ r e x p \delta_{\rm rexp} δ rexp , exceeds a series of threshold values (shown on the x x x -axis). The middle panel shows the p p p -value that quantifies the likelihood that the B B B -mode is consistent with zero. As we reduce the δ r e x p \delta_{\rm rexp} δ rexp threshold, the p p p -value increases above p ≤ 0.01 p \le0.01 p ≤ 0.01 (red shaded region), and the masked survey passes the B B B -mode null test. As the mask reduces the total number of sources (lower panel), the blue squares show the p p p -values from a random mask analysis, demonstrating that the improved B B B -mode is not driven by a reduction in area. In the upper panel we find that the astrometric masking has a benign impact on our Σ 8 \Sigma_8 Σ 8 cosmological constraints. The green data point and shaded regions shows the results from our fiducial analysis, demonstrating that failing data-vectors are consistent with our fiducial, all having less than 0.4 σ 0.4\sigma 0.4 σ changes in Σ 8 \Sigma_8 Σ 8 .
Figure E.1 (PDF p. 30).
Furthermore, we also flip the shuffled masks in both the RA and Dec directions, to remove any possible correlations between the constructed masks and the focal plane. We then repeated the analysis using the randomised mask, and found that the B B B -mode signal remained significant (see the blue points in Fig. E.1). As such, we conclude that our astrometric mask is removing a bone fide B B B -mode signal from the data vector, and not just removing sensitivity to the B B B -mode by lowering the constraining power.
As such, from these tests we can see that our astrometric statistic targets our B B B -mode signal, and that our recovered cosmology is fairly insensitive to our choice of the mask. We therefore opt to define a fiducial threshold in our test statistic that passes our B B B -mode test and makes some intuitive sense. To that end, we define the fiducial mask to be where \delta_{\rm rexp} > 0\farcs2 , where 0\farcs2 is chosen as a round number that is smaller than the THELI pixel size (0\farcs212 ). This corresponds to approximately 4% of the full KiDS-Legacy source sample, drawn from a previously unmasked area of 46.6 square degrees.
Cosmological parameters from ξ ± \xi_\pm ξ ±
Here, we present cosmological parameter constraints from ξ ± \xi_\pm ξ ± , for comparison with previous work and external datasets. We followed the KiDS tradition of using nine logarithmic bins in θ \theta θ , but constructed these over the same scale cuts as used by our fiducial E n / B n E_n/B_n E n / B n : θ ∈ [ 2 , 300 ] ′ \theta\in[2,300]' θ ∈ [ 2 , 300 ] ′ . We also removed the first θ \theta θ -bin from ξ − \xi_- ξ − , in effect limiting it to θ ∈ [ 4 , 300 ] ′ \theta\in[4,300]' θ ∈ [ 4 , 300 ] ′ , following previous KiDS analyses. For our theoretical predictions we must also apply this binning scheme, however we note that it is not sufficient to simply integrate ξ ± ( θ ) \xi_\pm(\theta) ξ ± ( θ ) in Eq. (15) over the bin width. Instead, we must weight ξ ± ( θ ) \xi_\pm(\theta) ξ ± ( θ ) by the weighted number of galaxy pairs, N p a i r ( θ ) N_{\rm pair}(\theta) N pair ( θ ) [see 10 for details] .
We used the finely binned ξ ± ( θ ) \xi_\pm(\theta) ξ ± ( θ ) and its corresponding number of galaxy pairs to make our measurements and our predictions,
ξ ± b i n n e d ( θ b i n ) = ∑ N p a i r ( θ ) ξ ± ( θ ) Δ ( θ b i n − θ ) ∑ N p a i r ( θ ) Δ ( θ b i n − θ ) , (39) \xi_\pm^{\rm binned}(\theta_{\rm bin}) = \frac{\sum N_{\rm pair}(\theta)\,\xi_\pm(\theta)\,\Delta(\theta_{\rm bin}-\theta)}{\sum N_{\rm pair}(\theta)\,\Delta(\theta_{\rm bin}-\theta)} \, , \tag*{(39)} ξ ± binned ( θ bin ) = ∑ N pair ( θ ) Δ ( θ bin − θ ) ∑ N pair ( θ ) ξ ± ( θ ) Δ ( θ bin − θ ) , (39)
where Δ ( θ b i n − θ ) \Delta(\theta_{\rm bin}-\theta) Δ ( θ bin − θ ) is the binning function, which is equal to 1 if θ \theta θ is within the boundaries of the bin labelled θ b i n \theta_{\rm bin} θ bin , and is 0 otherwise. When predicting ξ ± b i n n e d ( θ b i n ) \xi_\pm^{\rm binned}(\theta_{\rm bin}) ξ ± binned ( θ bin ) , we assumed that the distribution of the number of pairs and the value of ξ ± ( θ ) \xi_\pm(\theta) ξ ± ( θ ) were both essentially constant within each of the 1000 bins, and used the predicted values of ξ ± \xi_\pm ξ ± at the centres of the bins from Eq. (15) in Eq. (F.1).
Figure F.1 presents the whisker diagram for our ξ ± \xi_\pm ξ ± constraints under all analyses, and Fig. F.2 presents our posterior predictive distributions for ξ ± \xi_\pm ξ ± . Despite our concerns regarding the possible systematic contamination of our ξ ± \xi_\pm ξ ± signals, our cosmological parameter constraints from ξ ± \xi_\pm ξ ± are in good agreement with the results from our other statistics. We note, however, that our posterior predictive distributions (shown in Fig. F.2) are much closer to the significance limit (PTE < 0.01 < 0.01 < 0.01 ) than our other statistics (see results for the fiducial setup in Table F.1, and for the analysis variations in Table H.1). Finally, we note that our companion paper Stölzner et al. [128 ] finds strong evidence for a scale-dependent tension between cosmological parameters modelled with ξ ± \xi_\pm ξ ± . Furthermore, Oehl & Tröster (2024) found that the likelihood of ξ ± \xi_\pm ξ ± becomes non-Gaussian on large scales. This is evidence that our choice to present constraints with ξ ± \xi_\pm ξ ± as less reliable was warranted.
Table F.1. Constraints on S 8 S_8 S 8 (α = 0.50 \alpha= 0.50 α = 0.50 ) and Σ 8 \Sigma_8 Σ 8 (α \alpha α free) for ξ ± \xi_{\pm} ξ ± under our fiducial setup.
Figure F.1. Whisker diagram of Σ 8 \Sigma_8 Σ 8 cosmological results from ξ ± \xi_{\pm} ξ ± . The figure is annotated as in Figure 11 .
Figure F.2. Data vector and posterior predictive distributions for the fiducial run of our ξ ± \xi_{\pm} ξ ± . Scales in the grey shaded region are excluded from our fiducial analysis. TPDs are shown for 68% (red) and 95% (orange) percentiles of the posterior models.
Additional posterior constraints
In this appendix we provide additional posterior constraints, for various combinations of parameters and model setups.
Whisker diagram in S 8 S_8 S 8
In this section we show the S 8 S_8 S 8 equivalent of our main whisker-diagrams, for comparison with previous work. Figure G.1 shows the constraints for our E n E_n E n , CE, and ξ ± \xi_{\pm} ξ ± statistics, with all analysis variations included. The variability between the constraints is larger in this projection than in the Σ 8 \Sigma_8 Σ 8 projections, due to shifting of the posterior mode along the σ 8 \sigma_8 σ 8 , Ω m \Omega_m Ω m degeneracy. Otherwise, we note that all conclusions of this work are unchanged when analysing the S 8 S_8 S 8 space, rather than the Σ 8 \Sigma_8 Σ 8 space.
Redshift distributions
Figure G.2 shows the posterior constraints on our redshift distribution bias parameters, for chains analysed using various intrinsic alignment models, compared to the prior. The chains are all consistent with the prior, with only the low-z z z tomographic bins showing any tendency to shift away from the prior. Moreover, the posteriors are all fully consistent between the various analyses, demonstrating that the choice of intrinsic alignment model does not have a significant impact over the inferred redshift distribution biases.
Cosmological parameters
Figure G.3 shows the posterior constraints on our cosmological parameters, again for chains analysed using various intrinsic alignment models, and again compared to the prior. Only S 8 S_8 S 8 and Ω m \Omega_m Ω m can be reliably considered to be constrained, although there are interesting preferences in certain parameters.
Figure G.1. Whisker diagram of S 8 S_8 S 8 cosmological results presented in this work. Each panel shows the compilation of results for one of our two-point statistics, as annotated. In each panel, we show the fiducial constraints in the relevant statistic with coloured points (colour-coded by the statistic), and highlight the fiducial marginal HPDI with a band of the equivalent colour. We also highlight the marginal HPDI in S 8 S_8 S 8 from Planck-Legacy (shaded red). All variations to our analysis are annotated, and values plotted in dark blue. We annotate each analysis variation with the difference (quantified using the Gaussian approximation to the Hellinger tension) that the resulting marginal mode/HPDI has with respect to the fiducial (dark blue) or Planck (red).
Figure G.2. Posterior constraints on our redshift distribution bias parameters for chains utilising our various intrinsic alignment models. The filled contour shows the prior derived from our SKiLLS simulations. Smoothing kernels are matched between all contours (and marginal distributions), and are shown by the non-zero width of the σ z = 0 \sigma_z = 0 σ z = 0 contours/marginals.
Figure G.3. Posterior constraints on our cosmological parameters for chains utilising our various intrinsic alignment models. The filled contour shows the prior on each cosmological parameter.
Tabulated constraints
Table H.1. Constraints on Σ 8 \Sigma_8 Σ 8 under our analysis variations.
Post-unblinding analyses and changes
In this appendix we detail all analyses that were taken, and any changes that were made to the analysis and/or manuscript, after unblinding.
Comparison between Legacy and KiDS-1000
As discussed in Sect. 6, KiDS-Legacy presents a decrease in the marginal Hellinger tension with Planck of 2.2 σ 2.2\sigma 2.2 σ for Σ 8 \Sigma_8 Σ 8 (1.65 σ 1.65\sigma 1.65 σ for S 8 S_8 S 8 ). After the unblinding of KiDS-Legacy we investigated the differences between the Legacy and KiDS-1000 samples and analyses, to determine whether the change in cosmological parameters could be attributed to systematic or random variations. To do this, we reanalysed the KiDS-1000 and KiDS-Legacy shape catalogues with a series of analysis variations, which changed the sample used for cosmological inference and/or the analysis choices that were used in the inference. This required us, in each case, to rerun our N ( z ) N(z) N ( z ) estimation, our correlation function measurement, covariance computation, and posterior inference. We did not, however, remeasure the redshift distribution bias parameters or the multiplicative shape measurement biases for every setup.
Figure I.1. Variation in recovered E n E_n E n cosmological constraints for Σ 8 \Sigma_8 Σ 8 (left) and S 8 S_8 S 8 (right) when moving between KiDS-1000 and KiDS-Legacy analyses, relative to the constraints from Planck (red). Each point is annotated with the Hellinger tension between the marginal constraint and that of Planck .
Figure I.1 presents a graphical summary of the changes we found in our recovered cosmological constraints when analysing KiDS-1000 using Legacy analysis choices and calibration methods, and when analysing Legacy using subsets defined by the KiDS1000 footprint and photo-z z z baseline. Our investigation of the consistency between the data from Li et al. (2023a) and Legacy was focussed on two areas: the impact of Legacy analysis choices and calibration methods, and the impact of the new on-sky area available to Legacy. We quantified the effect of each analysis variation by recomputing the Hellinger tensions with respect to Planck , as this allows us to circumvent the complex covariance between various analyses of KiDS-1000 and KiDS-Legacy.
As a result of this investigation, we found no evidence of statistically significant systematic biases in the sample and analysis of KiDS-Legacy. We determined that the primary differences between the cosmological results of Li et al. (2023a) and Legacy are attributable to changes to our redshift calibration sample and methodology, and to statistical noise driven by the new area available to KiDS-Legacy outside the KiDS-1000 footprint. We have made no changes to the analysis methodology of Legacy as presented in the main text after this investigation, and thus conclude that this suite of tests do not represent a possible source of bias in our blinded analysis.
Analysing KiDS-1000 with Legacy calibration methods
We started by investigating the impact of the updated redshift estimation and calibration methods of Wright et al. [138 ] on the cosmological constraints of Li et al. (2023a). To do this we reran the entirety of the cosmic shear measurement from KiDS-1000, from raw shape catalogues to cosmological inference, incorporating our new redshift estimation and calibration methodologies (leveraging e.g. gold-weight and simulation construction with sample matching). In this analysis we computed revised redshift bias estimates using the original KiDS-1000 redshift calibration sample, but updating to SKILLS image simulations and gold-weight redshift distribution estimation. The result of this end-to-end reanalysis are provided in Fig. I.Figure 1 as ‘Legacy Nz est. & calib.’. We found that the updated redshift distribution calibration method brings the cosmological estimates from the Li et al. (2023a) sample into better agreement with Planck , at 2.08 σ 2.08\sigma 2.08 σ consistency in Σ 8 \mathrm{\Sigma}_8 Σ 8 and 1.42 σ 1.42\sigma 1.42 σ consistency in S 8 S_8 S 8 . Importantly, both of these measurements satisfy the typical KiDS threshold for consistency of 2.3 σ 2.3\sigma 2.3 σ (equivalent to a one-sided probability to exceed of 0.01).
Updating KiDS-1000 to the Legacy calibration sample
We also measured the contribution to the difference between KiDS-1000 and KiDS-Legacy caused by the differences to our analysis choices and calibration sample, including our new tomographic bin limits and our greatly expanded redshift calibration dataset. To do this we again performed a full end-to-end reanalysis of KiDS-1000, this time with our new tomographic bin limits (excluding the sixth bin), and new redshift calibration sample from KiDZ. In this analysis we assume the redshift distribution biases and multiplicative shear biases are the same as the fiducial KiDS-Legacy analysis. We find that the use of the larger KiDZ calibration sample also slightly reduces the offset with Planck , to 1.61 σ 1.61\sigma 1.61 σ consistency in Σ 8 \mathrm{\Sigma}_8 Σ 8 and 1.12 σ 1.12\sigma 1.12 σ consistency in S 8 S_8 S 8 (‘Legacy Tomo. & KiDZ (5 bins)’). This demonstrates that the KiDS-1000 sample, when analysed with our larger spectroscopic calibration sample and new redshift calibration methods (but without adding higher redshift sources, or using updated imaging, or new photo-z z z , or additional area) is consistent with the results from Planck and with the results from Legacy.
Legacy analysis in the KiDS-1000 volume
We explored the consistency between the constraints from KiDS-1000 and KiDS-Legacy by analysing the two samples in the same volume, and with the same analysis choices. We were able to do this fairly easily thanks to our analysis from Sect. I.1.2, which homogenised the tomography and N ( z ) N(z) N ( z ) estimation and calibration between the two samples. By restricting KiDS-Legacy to only the footprint of KiDS-1000, excluding the sixth tomographic bin, and using the NLA IA model, we were able to perform a comparison that is mostly free of statistical fluctuations. In fact there are still differences in the source samples due to changes in source extraction, image masking, source selection, and more, which we assumed add only a small amount of statistical noise to the comparison. We find that the Legacy sample in the same volume as KiDS-1000 is slightly closer to Planck than KiDS-1000, with a 1.17 σ 1.17\sigma 1.17 σ consistency in Σ 8 \mathrm{\Sigma}_8 Σ 8 and 0.85 σ 0.85\sigma 0.85 σ consistency in S 8 S_8 S 8 . This indicates that the sum total of changes (such as updated imaging calibrations, photo-z z z , masking, etc.) lead to a 0.44 σ 0.44\sigma 0.44 σ (0.27 σ 0.27\sigma 0.27 σ ) decrease in the Hellinger distance between Planck and KiDS; significantly less than the change introduced by our new calibration methods.
Sixth tomographic bin and NLA-M M M
We further investigated the influence of adding the sixth tomographic bin and our new NLA-M M M IA model to the analysis of Legacy in the footprint of KiDS-1000 (‘K1000 Tiles; NLA-M; 6 bins’), finding that these changes caused a significant reduction in our marginal uncertainties (as expected) but did not lead to a significant change in the estimated cosmological parameters. As such, our Hellinger distance with Planck increases with these additions, to 1.40 σ 1.40\sigma 1.40 σ consistency in Σ 8 \mathrm{\Sigma}_8 Σ 8 and 1.00 σ 1.00\sigma 1.00 σ consistency in S 8 S_8 S 8 .
Additional Legacy area
Finally, we added the data back into Legacy that resides outside the footprint of KiDS-1000, and recover our fiducial result, which exhibits a 0.73 σ 0.73\sigma 0.73 σ consistency in Σ 8 \mathrm{\Sigma}_8 Σ 8 and 0.74 σ 0.74\sigma 0.74 σ consistency in S 8 S_8 S 8 . This demonstrates that there is a non-negligible statistical noise component in the difference between the sources probed within the footprint of KiDS-1000 and over the final KiDS-Legacy area.
Expanded scale cuts
Post unblinding, we discovered a typographical error in the pipeline which performed the measurement of correlation functions with expanded scale cuts (Sect. 5.2.7). This error meant that the pre-unblinding measurements were made with uncalibrated shapes directly from our shape measurement code lensfit , rather than with the recalibrated shapes (Sect. 3.5). We reran the expanded scale cut analysis with the correct recalibrated shapes post-unblinding, finding that the consistency between the analyses with fiducial and expanded scale cuts improved from 0.3 σ 0.3\sigma 0.3 σ to 0.15 σ 0.15\sigma 0.15 σ with the use of the correct shapes. No conclusions were changed as a result of this error.
Iterative covariances
All fiducial chains were recomputed using our iterative covariance framework post-unblinding. This was a conscious pre-unblinding choice, as the iterative covariances require non-negligible CPU time for the computation of MAP and new covariances. We found that this process had a negligible impact on our constraints (see Sect. 5.2.3), and thus does not represent a possible source of bias in our blinded analysis.