© The Authors 2025

1. 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, ΛCDM. It is a strikingly simple model based on only a handful of numerical parameters (Planck Collaboration VI 2020; 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 (Λ, 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), 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; 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 LSS and their CMB-based predictions has revealed a tension in the $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{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₈ 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. 2016) 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. 2013), as well as from previous results from the ongoing ‘stage-III’ surveys by the Kilo-Degree Survey (KiDS; Asgari et al. 2021; Li et al. 2023a), the Dark Energy Survey (DES; Amon et al. 2022; Secco et al. 2022), the Hyper-Suprime Camera survey (HSC; Li et al. 2023b; Dalal et al. 2023), and a joint analysis by Dark Energy Survey and Kilo-Degree Survey Collaboration (2023). Other low-z LSS probes typically also share this trend but again with a varying level of significance (see Abdalla et al. 2022, 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 2024, for a graphical summary). This makes cosmic shear somewhat more sensitive than other low-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₈, rather than a more fundamental cosmological cause (Yoon et al. 2019; Yoon & Jee 2021; Amon & Efstathiou 2022; Preston et al. 2023).

A definitive answer to the reality of the S₈ 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 (Euclid Collaboration: Mellier et al. 2025), the Legacy Survey of Space and Time (LSST) conducted with the Vera C. Rubin Observatory (Ivezić et al. 2019), the NASA Roman Space Telescope (Spergel et al. 2015), and the Chinese Space Station Optical Survey (Gong et al. 2019). 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 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. 2024), 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-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_n-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. (2025) 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. (2025) describes the various redshift distribution and bias estimation developments made for KiDS-Legacy.

2. 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. (2015). 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.

2.1. Shear power spectrum

Our two-point statistic predictions are calculated from shear (cross-)power spectra of galaxy samples i and j, $C_{\mathrm{GG}}^{(ij)}(\ell )$, 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 1992; LoVerde & Afshordi 2008; Kilbinger et al. 2017) taking the form

$$\begin{array}{c}\hfill C_{\mathrm{GG}}^{(ij)}(\ell )={\displaystyle {\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},\mathrm{nl}}(\frac{\ell +1/2}{f_K(\chi )},z(\chi )),\end{array}$$(1)

where P_m, nl is the non-linear matter power spectrum and 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, χ_hor, in practice, the weak lensing kernel, $W_{\mathrm{G}}^{(i)}(\chi )$, limits the integral to the comoving distance of the furthest galaxy in our sample, expressed as

$$\begin{array}{c}\hfill W_{\mathrm{G}}^{(i)}(\chi )=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{f_K(\chi )}{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{hor}}}}\mathrm{d}\chi \prime n_{\mathrm{S}}^{(i)}(\chi \prime )\frac{f_K(\chi \prime -\chi )}{f_K(\chi \prime )},\end{array}$$(2)

where $n_{\mathrm{S}}^{(i)}(\chi )$ is the distribution of comoving distances of galaxies in source sample, i; in practice, we express this in terms of redshift, namely, N(z) dz = n_S(χ) dχ. Here, Ω_m is the total matter density parameter today, H₀ is the Hubble constant, c is the speed of light, and a(χ) is the scale factor normalised to unity today.

2.2. Matter power spectrum

In our cosmic shear analysis, we remained mostly sensitive to Fourier scales 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. 2000; Howlett et al. 2012) to estimate the linear matter power spectrum and the augmented halo model ‘HMCODE2020’ (Mead et al. 2021) 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. 2023; Aricò et al. 2021) and is twice as accurate as HALOFIT from Takahashi et al. (2012), which has been used in a number of previous cosmic shear analyses (see e.g. Amon et al. 2022; Hikage et al. 2019). Dark Energy Survey and Kilo-Degree Survey Collaboration (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₈, resulting in an underestimated S₈ value. Similarly, an earlier version of HMCODE, HMCODE2015 (Mead 2015), which was used in previous KiDS analyses, is ∼30% less accurate than HMCODE2020. This can subsequently produce small differences between cosmological constraints estimated using the two codes: Tröster et al. (2021) 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₈ 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²Ω_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.

2.3. 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/Mpc ≲ k ≲ 0.5 h/Mpc), BAOs produce oscillations in the power spectrum. Over a wide range of intermediate scales (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/Mpc. At even smaller scales (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_AGN) ∈ [7.1,8.3]). This feedback model is then calibrated against a suite of hydrodynamical simulations from BAHAMAS (McCarthy et al. 2017), 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 (Schaller et al. 2025), 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 (2015), Aricò et al. (2020) presented perturbative models that act directly on the outputs of dark-matter only simulations and can be empirically calibrated (Grandis et al. 2024; Schneider et al. 2022). We did not utilise these models for our cosmological analyses here, as these models were published over the course of the analysis of KiDS-Legacy or 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.

2.4. 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

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

where $C_{\mathrm{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:

$$\begin{array}{cc}\hfill C_{\mathrm{II}}^{(ij)}(\ell )& ={\displaystyle {\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}}^{(ij)}(\frac{\ell +1/2}{f_K(\chi )},z(\chi )).\hfill \end{array}$$(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 (Hirata & Seljak 2004):

$$\begin{array}{cc}\hfill C_{\mathrm{GI}}^{(ij)}(\ell )& ={\displaystyle {\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}_{\delta \mathrm{I},X}^{(j)}(\frac{\ell +1/2}{f_K(\chi )},z(\chi )),\hfill \end{array}$$(5)

with $C_{\mathrm{IG}}^{(ij)}=C_{\mathrm{GI}}^{(ji)}$. In Eqs. (4) and (5) the power spectra P_II, X and 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, where the Limber approximation is sufficient even for more compact kernels (Leonard et al. 2023).

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 (Joachimi et al. 2008). 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.

2.4.1. The non-linear linear alignment (NLA) model

The so-called NLA model (Bridle & King 2007) assumes that the intrinsic ellipticity of a galaxy depends linearly on its surrounding tidal quadrupole (Catelan et al. 2001). 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 (Hirata & Seljak 2004). This leads to the following expressions for the IA power spectra,

$$\begin{array}{cc}\hfill P_{\delta \mathrm{I},\mathrm{NLA}}(k,z)& =-A_{\mathrm{IA}}\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \\ \hfill P_{\mathrm{II},\mathrm{NLA}}(k,z)& =A_{\mathrm{IA}}^2{\left(\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}\right)}^2{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \end{array}$$(6)

where C₁ρ_cr ≈ 0.0134 is a constant and D(z) is the linear growth factor normalised to unity at z = 0. The model is empirically extended to non-linear scales by employing the full, non-linear matter power spectrum 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 2010). The model has a single free parameter, namely, the global, dimensionless amplitude 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_IA is well-constrained by tomographic cosmic shear data, we chose a wide top-hat prior for 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 2023).

2.4.2. Additional redshift dependence (NLA-z)

There is currently no observational evidence of redshift evolution in the amplitude of IA (Joachimi et al. 2011; Fortuna et al. 2021a), 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_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_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. 2019), we chose an extended model that is linear in terms of the scale factor,

$$\begin{array}{cc}\hfill P_{\delta \mathrm{I},\mathrm{zdep}}^{(i)}(k,z)& =-(A_{\mathrm{IA}}+B_{\mathrm{IA}}[\frac{{⟨a⟩}^{(i)}}{a_{\mathrm{piv}}}-1])\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \\ \hfill P_{\mathrm{II},\mathrm{zdep}}^{(ij)}(k,z)& =(A_{\mathrm{IA}}+B_{\mathrm{IA}}[\frac{{⟨a⟩}^{(i)}}{a_{\mathrm{piv}}}-1])\hfill \\ \hfill & \times (A_{\mathrm{IA}}+B_{\mathrm{IA}}[\frac{{⟨a⟩}^{(j)}}{a_{\mathrm{piv}}}-1]){\left(\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}\right)}^2{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \end{array}$$(7)

where ⟨a⟩⁽ⁱ⁾ is the average scale factor in tomographic bin i, evaluated using the estimated N(z). The parameter A_IA thus is the IA amplitude at a_piv ≈ 0.769, which corresponds to the commonly chosen pivot of z = 0.3 for IA studies (Joachimi et al. 2011). We set the prior range of the second parameter, B_IA, by fitting our linear amplitude dependence to the compilation of KiDS and SDSS measurements studied in Fortuna et al. (2021a, 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_IA with mean −3.7 and standard deviation 4.3. For direct comparisons between amplitudes estimated between this model and others, we define

$$\begin{array}{c}\hfill A_{\mathrm{IA},\mathrm{total}}^{\mathrm{NLAz},(i)}={\widehat{A}}_{\mathrm{IA}}+{\widehat{B}}_{\mathrm{IA}}(\frac{{⟨a⟩}^{(i)}}{a_{\mathrm{piv}}}-1),\end{array}$$(8)

where $\widehat{·}$ denotes the best-fitted (i.e. maximum a posteriori) value of each parameter.

2.4.3. 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 2004), accounting for the fact that alignments are only measured at local over-densities of the matter distribution where galaxies reside; and (ii) quadratic alignments as expected for tidal torquing effects (Crittenden et al. 2001; Catelan et al. 2001).

These contributions are subsumed into the tidal alignment, tidal torquing (TATT) model (Blazek et al. 2019), each with a free amplitude parameter. As can be seen from Figure 1 of Blazek et al. (2019), 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. 2019). 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 model is then expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_{\delta \mathrm{I},\mathrm{scaledep}}(k,z)& =P_{\delta \mathrm{I},\mathrm{NLA}}(k,z)\hfill \\ \hfill & +A_{\mathrm{IA}}{b}_{\mathrm{src}}\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}[A_{0|0E}(k,z)+C_{0|0E}(k,z)],\hfill \\ \hfill P_{\mathrm{II},\mathrm{scaledep}}(k,z)& =P_{\mathrm{II},\mathrm{NLA}}(k,z)\hfill \\ \hfill & +A_{\mathrm{IA}}^2{b}_{\mathrm{src}}^2{\left(\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}\right)}^2{A}_{0\mathrm{E}|0\mathrm{E}}(k,z)\hfill \\ \hfill & +2A_{\mathrm{IA}}^2{b}_{\mathrm{src}}{\left(\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}\right)}^2[A_{0|0\mathrm{E}}(k,z)+C_{0|0\mathrm{E}}(k,z)],\hfill \end{array}\end{array}$$(9)

where the expressions for A_0|0E, C_0|0E, and A_0E|0E are given in Blazek et al. (2019). 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_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] to ensure it encompasses NLA (b_src = 0) and reasonable values expected for the (unknown) source galaxy bias (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,

$$\begin{array}{c}\hfill A_{\mathrm{IA},\mathrm{total}}^{\mathrm{scaledep}}={\widehat{A}}_{\mathrm{IA}}.\end{array}$$(10)

2.4.4. 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 (Samuroff et al. 2023; Johnston et al. 2019; Heymans et al. 2013; Georgiou et al. 2025; see also Mandelbaum et al. 2011; Tonegawa et al. 2025; McCullough et al. 2024 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 (Joachimi et al. 2011; Fortuna et al. 2025). 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. (2018) 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 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⟩ within a tomographic bin. This leads to

$$\begin{array}{cc}\hfill P_{\delta \mathrm{I},\mathrm{massdep}}^{(i)}(k,z)& =-A_{\mathrm{IA}}{f}_{\mathrm{r}}^{(i)}{\left(\frac{{⟨M_{\mathrm{h}}⟩}^{(i)}}{M_{\mathrm{h},\mathrm{piv}}}\right)}^{\beta }\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \\ \hfill P_{\mathrm{II},\mathrm{massdep}}^{(ij)}(k,z)& =A_{\mathrm{IA}}^2{f}_{\mathrm{r}}^{(i)}{f}_{\mathrm{r}}^{(j)}{\left(\frac{{⟨M_{\mathrm{h}}⟩}^{(i)}}{M_{\mathrm{h},\mathrm{piv}}}\right)}^{\beta }{\left(\frac{{⟨M_{\mathrm{h}}⟩}^{(j)}}{M_{\mathrm{h},\mathrm{piv}}}\right)}^{\beta }\hfill \\ \hfill & \times {\left(\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(z)}\right)}^2{P}_{\mathrm{m},\mathrm{nl}}(k,z),\hfill \end{array}$$(11)

where we adopt the pivot halo mass M_h, piv = 10^13.5 h⁻¹ M_⊙ from Fortuna et al. (2025). 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_IA and β as our prior for cosmic shear inference. Their posterior is well approximated by a bivariate Gaussian with A_IA = 5.74 ± 0.29, β = 0.44 ± 0.03, and a correlation coefficient of −0.59. Details on how we determine 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

$$\begin{array}{c}\hfill A_{\mathrm{IA},\mathrm{total}}^{\mathrm{massdep},(i)}={\widehat{A}}_{\mathrm{IA}}{f}_r^{(i)}{\left(\frac{{⟨M_{\mathrm{h}}⟩}^{(i)}}{M_{\mathrm{h},\mathrm{piv}}}\right)}^{\widehat{\beta }}.\end{array}$$(12)

2.5. 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-/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. (2019).

Our cosmological analysis here leverages two-point statistics: the complete orthogonal sets of E/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).

2.5.1. 2PCFs (ξ_±)

The shear two-point correlation functions (ξ_±, Kaiser 1992) are functions of the tangential and cross components of shear (γ_t and γ_×, respectively) defined with respect to the line connecting pairs of galaxies with angular separation θ (see e.g. Bartelmann & Schneider 2001),

$$\begin{array}{c}\hfill {\xi }_{\pm }(\theta )=⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩(\theta )\pm ⟨{\gamma }_{\times }{\gamma }_{\times }⟩(\theta ).\end{array}$$(13)

The measured 2PCFs are binned in θ to reduce measurement noise and compress the data. In practice, we measure ξ_±(θ) for 1000 logarithmically spaced θ-bins between θ_min = 2′ and θ_max = 300′ using TREECORR¹ (Jarvis et al. 2004). The inputs to TREECORR are galaxy catalogues, which include galaxy ellipticities (ε₁ and ε₂), galaxy positions (right ascension and declination), and galaxy weights (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 ξ_±(θ) are then estimated using

$$\begin{array}{c}\hfill {\widehat{\xi }}_{\pm }^{(ij)}(\overline{\theta })=\frac{{\sum }_{\mathit{ab}}{w}_a{w}_b[{\epsilon }_{\mathrm{t},a}^{\mathrm{obs}}{\epsilon }_{\mathrm{t},b}^{\mathrm{obs}}\pm {\epsilon }_{\times ,a}^{\mathrm{obs}}{\epsilon }_{\times ,b}^{\mathrm{obs}}]{\mathrm{\Delta }}_{\mathit{ab}}^{(ij)}(\overline{\theta })}{{\sum }_{\mathit{ab}}{w}_a{w}_b(1+m_a)(1+m_b){\mathrm{\Delta }}_{\mathit{ab}}^{(ij)}(\overline{\theta })},\end{array}$$(14)

where ${\mathrm{\Delta }}_{\mathit{ab}}^{(ij)}(\overline{\theta })$ is the bin selection function (which limits the sums to galaxy pairs of separation within the angular bin labelled by $\overline{\theta }$ and the tomographic bins i and j), and ${\epsilon }_{\mathrm{t}}^{\mathrm{obs}}$ and ${\epsilon }_{\times }^{\mathrm{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 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:

$$\begin{array}{cc}\hfill {\xi }_{\pm }^{(ij)}(\theta )& ={\displaystyle {\int }_0^{\infty }}\frac{\ell \mathrm{d}\ell }{2\pi }[C_{\epsilon \epsilon ,\mathrm{E}}^{(ij)}(\ell )\pm C_{\epsilon \epsilon ,\mathrm{B}}^{(ij)}(\ell )]{\mathrm{J}}_{0/4}(\ell \theta ),\hfill \end{array}$$(15)

where $C_{\epsilon \epsilon ,\mathrm{E}/\mathrm{B}}^{(ij)}(\ell )$ are the E/B-mode angular power spectra for redshift bins i and j, and J_0/4(x) are the zeroth- and fourth-order Bessel functions of the first kind. We expect cosmic shear to produce only E-mode signals up to first order (Schneider et al. 1998; Hilbert et al. 2009), while data systematic errors can in principle create both E and 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. 2002; Hoekstra et al. 2002; Heymans et al. 2005; Kilbinger et al. 2013).

As ξ_± mix E and 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 ℓ probed by the ξ_± over a range of angular separations. The J₀ filter in Eq. (15) means that ξ₊ mixes a very wide range of scales (1 ≲ ℓ ≲ 10⁴) and J₄ makes ξ₋ particularly sensitive to smaller scales (10 ≲ ℓ ≲ 10⁵). This breadth can create a problem when using ξ_± 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 ξ_± in our primary cosmological analyses. Nevertheless, results for ξ_± are provided for completeness and comparison with previous work, in Appendix F.

Fig. 1. Integrands of the transformation between the angular power spectrum and En/Bn (Eq. 17, panel ‘a’), CE, B (Eq. C.5, panel ‘b’), and ξ± (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} for En and θ ∈ {2.0, 300.0} arcmin for ξ±). For CE, B we show all eight bins used in KiDS-Legacy, defined between ℓ = 100 and ℓ = 1500. Note that, unlike previous KiDS analyses, our new treatment of T(θ) means that the CE, B initially go to zero within their desired ℓ ranges.

2.5.2. COSEBIs (E_n/B_n)

The COSEBIs (Schneider et al. 2010) were originally designed to cleanly separate E- and B-modes over a finite θ-range. This was achieved through the use of filter functions T_±n(θ), which form complete bases, creating well-defined COSEBIs E- and B-modes (E_n/B_n). COSEBIs can be calculated or measured via a linear transformation of ξ_±,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{E}}_n& =\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\mathrm{d}\theta \theta [T_{+n}(\theta ){\xi }_{+}(\theta )+T_{-n}(\theta ){\xi }_{-}(\theta )],\hfill \\ \hfill {\widehat{B}}_n& =\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\mathrm{d}\theta \theta [T_{+n}(\theta ){\xi }_{+}(\theta )-T_{-n}(\theta ){\xi }_{-}(\theta )].\hfill \end{array}\end{array}$$(16)

The COSEBIs modes 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. 2012). This means that COSEBIs naturally compress the data into well-defined modes, once θ_min and θ_max are defined, without any need for subsequent data compression. In previous KiDS analyses (e.g. Asgari et al. 2021), 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₈. 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,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E_n& ={\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\ell \ell }{2\pi }{C}_{\mathrm{EE}}(\ell )W_n(\ell ),\hfill \\ \hfill B_n& ={\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\ell \ell }{2\pi }{C}_{\mathrm{BB}}(\ell )W_n(\ell ),\hfill \end{array}\end{array}$$(17)

where the weight functions W_n(ℓ) are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill W_n(\ell )& ={\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\mathrm{d}\theta \theta T_{+n}(\theta ){\mathrm{J}}_0(\ell \theta )\hfill \\ \hfill & ={\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\mathrm{d}\theta \theta T_{-n}(\theta ){\mathrm{J}}_4(\ell \theta ).\hfill \end{array}\end{array}$$(18)

Figure 1 shows two example integrands from Eq. (17). W_n(ℓ) are highly oscillatory weight functions with a limited range of support over ℓ. This allows us to truncate our integrals in Eq. (17) to a finite range in ℓ, as the E_n integrand in Eq. (17) is practically zero outside 10 ≲ ℓ ≲ 10³.

2.5.3. Band powers (C_E, B)

Band power spectra (hereafter ‘band-powers’) refer to binned angular power spectra measured by integrating over ξ_±(θ) (Schneider et al. 2002; van Uitert et al. 2016; Joachimi et al. 2021). Ideally, we want a finite number of bins, labelled as l, where each bin is only sensitive to a narrow range of ℓ-values between ℓ_lo, l and ℓ_up, l. We also want the E- and B-mode band-powers to be pure, namely, no leakage from B modes into E modes and vice versa. For this ideal case, we can express the following:

$$\begin{array}{c}\hfill {\mathcal{C}}_{\mathrm{E}/\mathrm{B},l}^{(ij)}:=\frac{1}{{\mathcal{N}}_l}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\ell \ell S_l(\ell )C_{\epsilon \epsilon ,\mathrm{E}/\mathrm{B}}^{(ij)}(\ell ),\end{array}$$(19)

for redshift bins i and j, where S_l(ℓ) is the response function for bin l. Here we have normalised the integral with

$$\begin{array}{c}\hfill {\mathcal{N}}_l={\displaystyle {\int }_0^{\infty }}\mathrm{d}\ell {\ell }^{-1}{S}_l(\ell ).\end{array}$$(20)

In this work, we consider the ideal case of a top-hat response function between ℓ_lo, l and ℓ_up, l which simplifies the normalisation to

$$\begin{array}{c}\hfill {\mathcal{N}}_l=log({\ell }_{\mathrm{up},l})-log({\ell }_{\mathrm{lo},l}),\end{array}$$(21)

and produces a band-power that traces ℓ²C(ℓ) at the logarithmic centre of the bin. We define eight band-powers logarithmically binned between ℓ_min = 100 and ℓ_max = 1500. In practice, measuring band-powers with perfect top-hat functions is infeasible, as ξ_±(θ) is only measured for a finite θ-range. This has the adverse effect that our band-powers will inevitably mix E and B modes². The mixing increases for narrower ℓ-bins and when the available range of θ-scales is smaller (Asgari & Schneider 2015). With this definition for band-powers, we forgo perfect E-/B-mode separation in favour of more control over the ℓ-values that contribute to our analyses. We present the detailed derivation of our band-power windows ${\mathcal{C}}_{\mathrm{E},l}$ in Appendix C, and show the integrands of ${\mathcal{C}}_{\mathrm{E},l}$ in Fig. 1. Note, however, that our integrands behave differently to those in Asgari et al. (2021), due to our different treatment of the edges of T(θ) (Appendix C) and due to our different scale-cuts.

2.6. Analytic covariance modelling

We presented the covariance modelling for KiDS-Legacy and the ONECOVARIANCE code³ in Reischke et al. (2025). In general, the methodology for the cosmic shear part still follows closely the covariance modelling of KiDS-1000 (Joachimi et al. 2021). 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, generally consists of four contributions:

$$\begin{array}{c}\hfill \mathbf{CC}={\mathbf{CC}}_{\mathrm{G}}+{\mathbf{CC}}_{\mathrm{NG}}+{\mathbf{CC}}_{\mathrm{SSC}}+{\mathbf{CC}}_{\mathrm{mult}},\end{array}$$(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_mult. Asgari et al. (2021), Secco et al. (2022) included the multiplicative shear bias in the theory prediction instead and marginalised over in the sampling process. If the residual uncertainties on m_a are small, these two approaches are equivalent. The Gaussian contribution can further be split into three parts:

$$\begin{array}{c}\hfill {\mathbf{CC}}_{\mathrm{G}}={\mathbf{C}}_{\mathrm{sva}}+{\mathbf{C}}_{\mathrm{sn}}+{\mathbf{C}}_{\mathrm{mix}},\end{array}$$(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. 2023) mocks to test the impact of variable depth and the overall agreement. These mocks are discussed in more detail in Appendix D. In Reischke et al. (2025), 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.

3. 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. (2020). 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) construction, shear calibration, data vector construction, covariance construction, cosmological inference, and post-processing. Details of each processing step are provided alongside the public code base⁴.

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 D and Yan et al. 2025). The weight is defined as described in Section 7 of Wright et al. (2024), and has the range $w\in [0,{\sigma }_{\mathrm{pop}}^{-2}]$, where σ_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. (2023c).

Fig. 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).

3.1. KiDS and VIKING

The KiDS-Legacy lensing sample was drawn from the fifth and final data release of the KiDS survey⁵, KiDS-DR5, described in Wright et al. (2024). It includes 1347 square-degree tiles observed in the u, g, r, and i filters with the VST (Capaccioli et al. 2005) and OmegaCAM (Kuijken 2011) between 2011 and 2019, as well as the 2009–2018 VIKING survey VISTA/VIRCAM observations (Edge et al. 2013) of the same area in Z, Y, J, H, and 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. 2019), which was the basis for the KiDS-1000 analyses (Asgari et al. 2021; Heymans et al. 2021; Tröster et al. 2021). It represents a substantial increase in survey area (by 34%) and in 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. 2021). The deeper 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. (2024). The astrometry reference catalogue for the THELI reduction moved from the Sloan Digital Sky Survey (SDSS Ahumada et al. 2020) and the 2-micron All Sky Survey (2MASS Skrutskie et al. 2006) to Gaia DR2 (Gaia Collaboration 2018), and the astrometric solutions and photometric zero-point determinations in the ASTRO-WISE 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. 2023c). Details and tests of these new implementations are given in Wright et al. (2024), 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.

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 deg², with n_eff = 8.79 per arcminute within the photometric redshift limits of KiDS-Legacy.

3.2. 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. 2007), the DEEP2 Galaxy Redshift Survey (DEEP2; Newman et al. 2013) 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. 2005) 14h field, and the VIMOS Public Extragalactic Redshift Survey (VIPERS; Guzzo et al. 2014) W1 and W4 fields (a full list is given in Wright et al. 2024). 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. 2022) survey, 2-degree Field Lensing Survey (Blake et al. 2016), WiggleZ survey (Blake et al. 2008), and the Dark Energy Spectroscopic Instrument (DESI Collaboration 2024).

3.3. Blinding

We used the same blinding methodology as the one implemented in Asgari et al. (2021), first presented in Kuijken et al. (2015). This process involves the double-blind generation of threerealisations 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σ change in the inferred S₈. The true catalogue is equally likely to be the one that yields the highest, middle, or lowest S₈. 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 afterunblinding.

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.

3.4. Redshift distributions and calibration

Redshift distributions (hereafter ‘N(z)’) used in KiDS-Legacy were estimated and calibrated leveraging both colour-based direct calibration using self-organising maps (SOMs Kohonen 1982) and spatial cross-correlations (CC) between tracer spectroscopic samples and the KiDS-Legacy sources. For a detailed description of the N(z) construction and validation process used in this work, we refer to Wright et al. (2025). In brief, the 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 times (where we arbitrarily chose 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 spectroscopiccompilation.

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. (2023c). For additional SOM redshift distribution calibration and cross-correlation redshift distribution calibration, we used the MICE2 simulated lightcone (Fosalba et al. 2015a; Crocce et al. 2015; Fosalba et al. 2015b; Carretero et al. 2015; Hoffmann et al. 2015) adapted to KiDS noise levels and sample selections (van den Busch et al. 2020).

Figure 3 shows the fiducial redshift distributions used in this cosmological analysis. These N(z) include the appropriate δz shifts in each of their mean, as inferred using our simulations and described in Wright et al. (2025). As such, the N(z) presented are shown as-used in the generation of the theory/models during parameter inference. In Wright et al. (2025) we quantified the impact of a series of different analysis choices for N(z) estimation and calibration, demonstrating that they have no significant impact on the final N(z) that we used for cosmological inference. The figure also presents the off-diagonal elements of the correlation matrices for δz, used to construct our correlated δ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.

Fig. 3. Left: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. The N(z) have been shifted according to the mean bias estimated using our SKiLLS simulations, making each N(z) a correct representation of the effective N(z) that is used in the cosmological analyses. Fiducial biases are provided in Table 1. Right:δz correlation matrix for our fiducial N(z) (δzSOM, lower triangle) and the correlation matrix inferred from our cross correlation N(z) on the data (DzCC, upper triangle).

3.5. 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. (2023c) 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. 2007, 2013; Fenech Conti et al. 2017), 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. 2024). 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. 2022), although this is determined to be of negligible importance to the computation of multiplicative bias (Li et al. 2023c) and redshift distribution bias (Wright et al. 2025) 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 Wright et al. 2025). Afterwards, these samples were also used to estimate the multiplicative shear measurement bias, 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 and δ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-terms’) in the dataset. We adopted the first-order systematics model from Heymans et al. (2006):

$$\begin{array}{c}\hfill {\epsilon }_k^{\mathrm{obs}}=(1+m_k)({\epsilon }_k^{\mathrm{int}}+{\gamma }_k)+{\alpha }_k{\epsilon }_k^{\mathrm{PSF}}+c_k,k\in \{1,2\},\end{array}$$(24)

where k denotes the individual ellipticity components, γ 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^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 Li et al. 2023c; Wright et al. 2024) and modelled PSF ellipticities.

Figure 4 shows the α_k and 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. 2023c). The residual α_k are all negligible, with the α₂ having a slightly larger residual amplitude overall (α₂ ≈ 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) 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: ${\widehat{c}}_k={⟨{\epsilon }_k⟩}_w$. The $\widehat{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-terms are small, but not necessarily zero.

Fig. 4. PSF leakage measured in the KiDS-Legacy lensing sample after recalibration of shapes, computed per tomographic bin and ellipticity component k, and following the first-order systematics model in Eq. (24).

As a further quality test, we computed the non-tomographic 2D c-terms in camera coordinates per ellipticity component, which is shown in Fig. 5. The 2D 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\stackrel{\prime }{.}2$), and the computed c-term scaled by the (also smoothed) measured ellipticity dispersion in each bin. The residual c-term can be seen to have amplitudes |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₁ was identified in Hildebrandt et al. (2020) 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. (2021), the amplitude of this signal (≲1% of the ellipticity dispersion) is negligible in the context of our cosmic shear measurements.

Fig. 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 and 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. (2020).

4. 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 Section 5. We performed three primary tests: an analysis of the Paulin-Henriksson et al. (2008) systematics model, which quantifies the impact of PSF modelling residuals on correlation functions; an analysis of Bacon et al. (2003) statistics, which test specifically for contamination of the galaxy shape measurements by the PSF; and an analysis of B modes, which tests for a systematic signal that cannot be generated by cosmic shear.

4.1. Paulin-Henriksson et al. (2008) model

To validate the fidelity of our shape measurement catalogue, we first compute the Paulin-Henriksson et al. (2008) 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:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ⟨{\epsilon }_{\mathrm{obs}}{\epsilon }_{\mathrm{obs}}⟩\simeq & (1+2\left[\frac{\overline{\delta T_{\mathrm{PSF}}}}{T_{\mathrm{gal}}}\right])\langle {\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}{\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}\rangle \hfill \\ \hfill & +{\left[\overline{\frac{1}{T_{\mathrm{gal}}}}\right]}^2⟨({\epsilon }_{\mathrm{PSF}}\delta T_{\mathrm{PSF}})({\epsilon }_{\mathrm{PSF}}\delta T_{\mathrm{PSF}})⟩\hfill \\ \hfill & +2{\left[\overline{\frac{1}{T_{\mathrm{gal}}}}\right]}^2⟨({\epsilon }_{\mathrm{PSF}}\delta T_{\mathrm{PSF}})(\delta {\epsilon }_{\mathrm{PSF}}{T}_{\mathrm{PSF}})⟩\hfill \\ \hfill & +{\left[\overline{\frac{1}{T_{\mathrm{gal}}}}\right]}^2⟨(\delta {\epsilon }_{\mathrm{PSF}}{T}_{\mathrm{PSF}})(\delta {\epsilon }_{\mathrm{PSF}}{T}_{\mathrm{PSF}})⟩,\hfill \end{array}\end{array}$$(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 Section 3.3 of Giblin et al. (2021), so that ⟨ab⟩ denotes ξ_±. Here, ε_obs is the observed ellipticity of a source, ${\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}$ is the noiseless and unbiased measurement of source ellipticity (which we approximate as being equal to the modelled shear: ${\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}={\gamma }_{\mathrm{theory}}$), and T is the R² object size estimate.

Figure 6 shows the various contributions to the additive systematic, $\delta {\xi }_{+}^{\mathrm{sys}}(\theta )$, from the Paulin-Henriksson et al. (2008) systematics model. The four terms from Eq. (25), shown in various colours, cause ⟨ε_obsε_obs⟩ to deviate from ⟨ε_perfectε_perfect⟩. The total systematic δξ₊^sys:

Fig. 6. Contributions to the additive systematic $\delta {\xi }_{+}^{\mathrm{sys}}(\theta )$ from the Paulin-Henriksson et al. (2008) 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 ξ+. 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.

$$\begin{array}{c}\hfill \delta {\xi }_{+}^{\mathrm{sys}}(\theta )⟨{\epsilon }_{\mathrm{obs}}{\epsilon }_{\mathrm{obs}}⟩-\langle {\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}{\epsilon }_{\mathrm{obs}}^{\mathrm{perfect}}\rangle ,\end{array}$$(26)

is given by the summation of these four terms⁶ and can be compared to the yellow band which encloses 10% of the uncertainty on the cosmic shear correlation function ξ₊, 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.

4.2. Bacon et al. (2003) statistics

We also used the method from Bacon et al. (2003) to infer the degree of PSF leakage into the two-point correlation function measurement using

$$\begin{array}{c}\hfill {\xi }_{\pm }^{\mathrm{sys}}=\frac{{⟨{\epsilon }^{\mathrm{obs}}{\epsilon }^{\mathrm{PSF}}⟩}^2}{⟨{\epsilon }^{\mathrm{PSF}}{\epsilon }^{\mathrm{PSF}}⟩}.\end{array}$$(27)

In Fig. 7, we show the value of ξ₊^sys in each tomographic bin autocorrelation, as a fraction of the theory prediction of the cosmic shear signal from our best-fit ξ_± cosmology (${\xi }_{+}^{\mathrm{\Lambda }\mathrm{CDM}}$). Results from the cross-correlation bins are quantitatively similar to the autocorrelation. As a reference, we also show 10% of the 1σ uncertainty on the cosmic shear signal as reported by the ONECOVARIANCE (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).

Fig. 7. Ratio of PSF contamination in the cosmic shear correlation function, ${\xi }_{+}^{\mathrm{sys}}$, to the predicted best-fit cosmic shear correlation function, ${\xi }_{+}^{\mathrm{\Lambda }\mathrm{CDM}}$. The statistic is shown relative to 10% of the 1σ uncertainty on our fiducial ξ+ 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;  ∀ j > i are quantitatively similar to the i − i signal).

4.3. B modes

Our third null test examined the significance of 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-mode signals of cosmological and astrophysical origin ought to be negligible. As such, any significant detection of B modes implies a contamination of the data vector by systematic effects, which may similarly (but undiagnosably) also impact the 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 modes (estimated via the COSEBIs 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-mode signal is presented in Fig. 8. The p-values for the measured B modes across all tomographic bin combinations are 0.04 in the case of 6 B_n modes, and 0.09 in the case of 20 modes, and thus both pass our required threshold of 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 modes, has a p-value below 0.01. However, with 21 tests (one per bin combination) it is expected that one will occasionally find 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-value, which passes.

Fig. 8. COSEBIs B modes measured in KiDS-Legacy. Each panel is annotated with the p-value of the B-mode signals for twenty modes (black) and our fiducial 6 modes (blue). The total 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-mode signals are highly correlated across modes within a tomographic bin, so readers are cautioned against so-called ‘χ-by-eye’.

5. 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. 2025) 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. 2022), and Planck.

5.1. Fiducial constraints

The fiducial analysis of KiDS-Legacy used two statistics (E_n, C_E) measured over $2\stackrel{\prime }{.}0\le \theta \le 300\prime$, an NLA-M intrinsic alignment model, correlated δz prior marginalisation, emulated COSMOPOWERP(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 2023) 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. 2023; Dark Energy Survey and Kilo-Degree Survey Collaboration 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 σ₈, S₈, and Ω_m. Marginal constraints over both S₈ and Σ₈ = σ₈(Ω_m/0.3)^α (with the exponent fitted to the Ω_m vs. σ₈ 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 (2023, Sect. 2.6) present a discussion of the benefits and detriments of each metric. We also report the χ² 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. (2021).

Fig. 9. Fiducial Ωm, S8 (left), and Ωm, σ8 (right) constraints for our fiducial two-point statistics (En and CE, 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σ and (outer) 2σ boundary levels.

Data vectors for the fiducial E_n and C_E are given in Fig. 10, alongside translated posterior distributions (TPDs; see Köhlinger et al. 2019), for the first application of these consistency metrics to cosmic shear datasets) drawn from the fiducial posterior computed for each statistic.

Fig. 10. Data vector and TPDs for the fiducial analysis of En (upper panel) and CE (lower panel). TPDs are shown as polygons spanning the 68th (red), and 95th (orange) percentiles of the posterior models. The En 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₈ = 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.{815}_{-0.021}^{+0.016}$, and the maximum a posteriori and PJ-HPD estimate is $S_8=0.{811}_{-0.015}^{+0.022}$. The consistency between the measurements is a reflection of the Gaussianity of the marginal posterior, and of the N-dimensional posterior (in the case of the PJ-HPD). This is 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 (1991). Testing of the HPDI method indicates that, for our posterior distributions in S₈, 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₈ constraints from different statistics. The variance in S₈ between our statistics is seemingly large, with our C_E and E_n differing by 0.55σ. However, examination of the 2D posterior distribution in Ω_m and S₈ (Fig. 9) indicates that the optimal degeneracy direction for KiDS-Legacy is somewhat different than the assumed α = 0.5 parametrised by S₈. 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₈ parameter (as translation along the degeneracy can give the false impression of internal inconsistency between final parameter constraints).

When fitting for the optimal α under the more general Σ₈ = σ₈(Ω_m/0.3)^α, we found the best fitting α is 0.58 for E_n and 0.60 for C_E. These significant deviations from α = 0.5 can cause offsets between derived values of S₈, as is seen between our C_E and 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 in particular, we find a significant ∼40% reduction in the relative marginal uncertainty between our S₈ and Σ₈ parameters. As such, we opt to preferentially present results in the Σ₈ parametrisation in this manuscript, especially when quantifying the impact of analysis choices on recovered constraints. We note, however, that the different α values between the different statistics mean that Σ₈ values recovered using each statistic are not directly comparable to one another.

5.2. 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 Σ₈ (that is, for the optimal Ω_m, σ₈-degeneracy) results presented in this work. The figure includes our fiducial results in each statistic, with constraints computed using the three methods previously 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 (Σ₈^fid − Σ₈^var)/σ_Σ8^var. Differences are computed between the Σ₈ values of each setup and the fiducial. We also annotate the fiducial constraints with their Hellinger distance (see e.g. Appendix G.1 of Heymans et al. 2021) to the constraints from the public, fiducial Planck-Legacy analysis (Planck Collaboration VI 2020). Finally, we present a similar figure for S₈ in Appendix G, and tables of all constraints are presented in Appendix H.

Fig. 11. Whisker diagram of Σ8 cosmological results presented in this work. The best-fitting α 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 from Planck-Legacy (shaded red). All variations to the analysis are annotated, with the difference relative to the fiducial result quantified using the metric (Σ8fid − Σ8var)/σΣ8var (listed in blue). The Hellinger tension between the fiducial analyses and Planck is annotated in red.

5.2.1. Modelling 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 Fig. 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 Σ₈ of 0.03σ. We consider this level of change to be consistent with chain-to-chain variations, and therefore negligible.

5.2.2. Gaussian covariance

We tested the impact of our covariance modelling on the recovered value of Σ₈ in Fig. 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 Σ₈ of between 12% (C_E) and 19% (E_n), and no change in the value of Σ₈ with respect to the fiducial (ΔΣ₈ = 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.

5.2.3. 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. 2018), 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) 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 (Reischke et al. 2017; Kodwani et al. 2019).

5.2.4. 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 Fig. 11.

First, we see that the impact of changing our IA models on the marginal constraints of S₈ and Σ₈ are mild. For our C_E statistic, we found maximal differences of 0.17σ across all model variations. Also, E_n values are more sensitive to our chosen IA model, showing maximal differences of 0.54σ (for the NLA and NLA-k models and Σ₈). 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 model: we see a consistent reduction in constraining power (compared to the fiducial NLA-M model) of 12 − 13% for Σ₈ across all statistics. For our NLA and NLA-k models, the reduction in constraining power compared to the fiducial is slightly smaller, at ∼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 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 χ² of the modified IA model is consistent with the fiducial (changes in PTE at the orderof 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_IA, total values (see Eqs. (8), (10), (12)). Figure 12 shows the value of A_IA, total inferred from our 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.

Fig. 12. Constraints on the total intrinsic alignment amplitude from the various intrinsic alignment models, per tomographic bin (where appropriate), computed using the En posteriors. Bands and error bars represent the 1σ posterior constraint on each model. Our NLA-M model (red) is shown per tomographic bin, with the point positioned at the mean photo-z of the bin.

For our NLA and NLA-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 and NLA-M models.

For our fiducial NLA-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)⁷. As such, the NLA-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_IA, total in our NLA-z and NLA-M models. Both show the expected increase in alignment amplitudes with increasing redshift, although the NLA-z model is at slightly lower amplitudes overall. We attribute this difference to the sixth tomographic bin: as the NLA-z model requires that all bins follow the same (monotonic) trend, any reduction in A_IA, total between the fifth and sixth tomographic bins will require a compensatory reduction in alignment amplitude across all bins.

5.2.5. 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 (σ_IA) and redshift and shear calibration uncertainty (σ_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 Σ₈ in the absence of marginalisation over these observational and systematic nuisance parameters are tighter by 10% (20%) for our E_n (C_E) measurements, and are shifted with respect to the fiducial by 0.28σ (0.14σ). This suggests that KiDS-Legacy is statistics-limited.

5.2.6. 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 ΔΣ₈ ≤ 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. 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 Fig. 11 the result of calibrating our redshift distributions using simulations from MICE2. These are the simulations that were used for KiDS-1000 N(z) calibration, but differ here in that we applied our KiDS-Legacy calibration methodology to these simulations. Wright et al. (2025) 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. As our bin six 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) and that from MICE2 is Δμ_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₈. As seen in Fig. 11, our analysis using N(z) calibrated with the MICE2 simulations increases the central value of our Σ₈ constraints, by up to 1.25σ (for 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. 2025). Our cosmological constraints from this analysis are essentially unchanged in their central value for both E_n and C_E, showing a 0.01σ shift in the central value of Σ₈ (computed with the fiducial α 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 Σ₈ are 23% (29%) larger for our CC N(z)/δ(z) marginalisation than in the fiducial case for E_n (C_E). Nonetheless, the strong agreement between the estimates made with our CC and SOM N(z) provides confidence that redshift distribution biases in KiDS-Legacy are under control.

5.2.7. Scale cuts

In KiDS-1000, Asgari et al. (2021) presented cosmological constraints using a range of angular scales from $0\stackrel{\prime }{.}5\le \theta \le 300\prime$. Subsequent investigation by Dark Energy Survey and Kilo-Degree Survey Collaboration (2023), however, found that an optimal treatment of baryonic systematic effect mitigation (for E_n/B_n) was achieved using $2\stackrel{\prime }{.}0\le \theta \le 300\prime$. For KiDS-1000 (Asgari et al. 2021), this resulted in a shift in the marginal parameter constraints of S₈ at the level of 0.8σ. Reanalysis of KiDS-1000 by Li et al. (2023a), with improved shape measurement and calibration from Li et al. (2023c), showed a reduced sensitivity to the choice of scale cut, seeing a shift of 0.3σ. We also tested the influence of changing our scale cuts, shown in the marginal projection of Ω_m, S₈ in Fig. 13, finding that including smaller scales results in a slight (0.1σ) decrease in the marginal value of Σ₈ and S₈. We also saw a slight change in the preferred value of our baryonic feedback parameter (see Fig. 13) where there is a reduction in posterior probability mass at negligible feedback amplitudes (i.e. log₁₀(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.

Fig. 13. Comparison of the marginal constraints on Ωm and S8 (top) and log10(TAGN/K) (bottom) when using a wider range of scales ($0\stackrel{\prime }{.}5\le \theta \le 300\prime$, 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 log10(TAGN/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 mode null test, causing it to once again fail our 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 systematics contamination is encouraging.

5.2.8. 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σ in Σ₈.(see e.g. Asgari et al. 2021). We performed the same test here, and show the results in Fig. 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 Σ₈ when removing individual bins of 0.24σ for E_n, 0.17σ for 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 (Stölzner et al. 2025), 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.

6. Discussion

The results from KiDS-Legacy are consistent with the results of Planck Legacy, with a marginal Σ₈ consistency of 0.73σ for E_n and 1.01σ for C_E. The constraints from KiDS-Legacy are also of comparable constraining power to those from Planck when projected into Σ₈, with Planck having a marginal uncertainty of ±0.015 compared to Legacy’s 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σ, 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σ) 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 Fig. I.1.

We also explored the consistency between our new Legacy cosmological parameters with recent results from the literature. Figure 14 presents the 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 (2023), and 2PCF measurements from HSC-Y3 (Li et al. 2023b). 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σ in Σ₈ and 0.41σ in S₈. 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σ and 1.8σ in Σ₈ 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₈ (or in the more constraining Σ₈) 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 (Stölzner et al. 2025). Lastly, as mentioned before, we do not include ξ_± into our fiducial analysis. However, the constraints and resulting consistency metrics are very similar and is discussed in Appendix F.

Fig. 14. Comparison between KiDS-Legacy En cosmological constraints and recent stage-III cosmic shear survey results from the literature. We show the 2D posteriors for each analysis Ωm, S8 (left) as well as the marginal projections in Σ8 (α = 0.58, middle), and S8 (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σ (where X is either Σ8 or S8), 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) and the amplitude of baryon feedback. Figure 15 shows our constraint over P(k) from the fiducial E_n, estimated by sampling the full posterior of KiDS-Legacy and generating posterior predictive distributions for P(k) using CAMB and HMCODE2020. We show the P(k) constraint normalised to the mean of all posterior samples, with a 0.1σ region shaded around this central reference. This allows us to directly compare our P(k) constraint to the distribution of P(k) returned by two power spectrum emulators, BACCO (Aricò et al. 2021) and EUCLIDEMULATOR2 (Euclid Collaboration: Knabenhans et al. 2021), 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. Nonetheless, we see that the differences between the emulated power spectra and those returned by CAMB are within the 0.1σ credible interval of our posterior constraint on P(k), demonstrating that differences between the various 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) will become a significant fraction of the finaluncertainty.

Fig. 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) 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 En/Bn 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 En/Bn is shown in red, relative to a range of hydrodynamical models of baryonic feedback suppression from FLAMINGO(Schaye et al. 2023).

We also explored the importance of baryon feedback modelling on our cosmological constraints, by examining the posterior on our feedback amplitude parameter log₁₀(T_AGN/K) in the context of our 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/Mpc. Above this scale our posterior constraint on the feedback model introduces a 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/Mpc, our marginal posterior on P(k) excludes the more extreme FLAMINGO feedback models. Our MAP estimate of log₁₀(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₁₀(T_AGN/K) is unconstrained (see Fig. 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, 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₁₀(T_AGN/K) = 7.3, while the MAP remains at the boundary (see Fig. 13). In the fiducial case our marginal posterior on log₁₀(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₁₀(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.

7. Summary

We present a cosmological analysis of the completed Kilo-Degree Survey (KiDS). Utilising data from KiDS Data Release 5 (Wright et al. 2024), 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. 2025). 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, 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. 2023c), 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), utilising conservative θ ∈ [2.0, 300.0] arcmin scales and a new physically-motivated intrinsic alignment (IA) model. Our fiducial E_n measurement constrains $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}=0.{815}_{-0.021}^{+0.016}$, which is in full agreement (Hellinger distance of 0.76σ) with the marginal prediction $S_8=0.{834}_{-0.016}^{+0.016}$ from the Planck Legacy cosmic microwave background experiment. Our fiducial measurement represents a ∼15% improvement in constraining power over the previous results from KiDS. However, we also found a significant change in optimal degeneracy between σ₈ and Ω_m, with Σ₈ ≡ σ₈(Ω_m/0.3)^α and α = 0.58 being preferred by our E_n measurement. We found ${\mathrm{\Sigma }}_8=0.{821}_{-0.016}^{+0.014}$, which is a ∼32% improvement in optimal constraining power compared to previous KiDS analyses. In this projection, the corresponding Hellinger distance to Planck is 0.73σ. We further computed the cosmological parameters using band power spectra (C_E) 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 (Stölzner et al. 2025) 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₁₀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.

Acknowledgments

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 RGSR2222268 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 20227RNLY3 “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 analyses of KiDS-Legacy.

References

Appendix A: 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. (2022). The original COSMOPOWER emulator was trained to reproduce the linear matter power spectrum, P_m, lin(k), and the perturbations caused by non-linear evolution, P_m, nl(k)/P_m, lin(k), given a set of relevant input parameters. We updated this process to emulate residuals log₁₀[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-dimensional⁸ 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 GitHub⁹, 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. (2022), across the wide range of scales that we emulate (k ∈ [10⁻⁵, 20] 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.

Appendix B: 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 and the red galaxy fraction 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, we used the output of the photometric redshift estimation from the template-fitting code BPZ (Benítez 2000). 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 encodes the best-fit templates and their admixture in steps of 0.1. We chose 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. 2011) surveys (i.e. for relatively bright and low-redshift galaxies), this T_B cut is capable of isolating the red sequence.

We found that in the first five tomographic bins 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 SKiLLS. The statistical error on f_r is very small. We repeated the measurement varying T_B by ±0.1 and obtained changes in f_r of less than 0.02. Therefore, we decided to consider f_r as fixed in the IA model.

To compute halo masses, we made use of the r-band luminosities derived as part of the KiDS Data Release 5 (Wright et al. 2024). van Uitert et al. (2015) 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

$$\begin{array}{cc}\hfill {log}_{10}{M}_{\mathrm{h}}^{(i)}[h^{-1}{M}_{\odot }]& =\gamma (\frac{1}{1+{⟨z⟩}^{(i)}}-\frac{1}{1+z_{\mathrm{piv}}})\hfill \\ \hfill & +\delta {log}_{10}\left(\frac{{⟨L_r⟩}^{(i)}}{L_0}\right)+\mathrm{\Delta }M,\hfill \end{array}$$(B.1)

where we chose z_piv = 0.3417, such that the parameters γ and ΔM are uncorrelated to a good approximation. An additional parameter δ describes the dependence on the average luminosity ⟨L_r⟩, normalised by the reference luminosity L₀ = 10¹¹(h/0.7)⁻² L_⊙. van Uitert et al. (2015) showed that the luminosity dependence evolves only mildly with redshift, so that we determine δ as the inverse variance-weighted average scaling over the four samples they considered. Thus, we obtained γ = 1.25 ± 0.33, δ = 1.42 ± 0.16, and Δ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₁₀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 a reduced error due to luminosities more compatible with the range fitted by van Uitert et al. (2015) 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. (2015). 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₁₀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 (Schneider & Bridle 2010). 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. 2022). For a given mass and redshift, early-type satellite fractions are typically twice as large as those of late types (Mandelbaum et al. 2006), 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. (2025) estimated a satellite fraction of 30 % and possibly higher (see also van Uitert et al. 2016) in a subset of the GAMA-KiDS sample and subsequently excluded it from their fit. Johnston et al. (2019) 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 (Fortuna et al. 2021b).

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 (Fortuna et al. 2021b) will be a natural choice to accurately represent different IA mechanisms and scalings on large and small scales.

Appendix C: Band-powers derivation

To arrive at the correct response functions for our band-powers, we first write ${\mathcal{C}}_{\mathrm{E}/\mathrm{B},l}$ as a function of ξ_±(θ), in the same manner as for the COSEBIs in Eq. (16),

$$\begin{array}{c}\hfill {\mathcal{C}}_{\mathrm{E}/\mathrm{B},l}=\frac{\pi }{{\mathcal{N}}_l}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\theta \theta T(\theta )[g_{+}^l(\theta ){\xi }_{+}(\theta )\pm g_{-}^l(\theta ){\xi }_{-}(\theta )],\end{array}$$(C.1)

where g_±^l(θ) are the band-power filter functions for bin l and T(θ) is a function that limits the range of θ-scales that contribute to the integral. If ξ_± is known for all scales then T(θ) = 1 everywhere. We now write g_±^l(θ) in terms of the idealised response function,

$$\begin{array}{c}\hfill g_{\pm }^l(\theta )={\displaystyle {\int }_0^{\infty }}\mathrm{d}\ell \ell S_l(\ell ){\mathrm{J}}_{0/4}(\ell \theta ).\end{array}$$(C.2)

If S_l(ℓ) is a top-hat then,

$$\begin{array}{cc}\hfill g_{+}^l(\theta )& =\frac{1}{{\theta }^2}[\theta {\ell }_{\mathrm{up},l}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{up},l})-\theta {\ell }_{\mathrm{lo},l}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{lo},l})],\hfill \\ \hfill g_{-}^l(\theta )& =\frac{1}{{\theta }^2}[{\mathcal{G}}_{-}(\theta {\ell }_{\mathrm{up},l})-{\mathcal{G}}_{-}(\theta {\ell }_{\mathrm{lo},l})],\hfill \end{array}$$(C.3)

with,

$$\begin{array}{c}\hfill {\mathcal{G}}_{-}(x)=(x-\frac{8}{x}){\mathrm{J}}_1(x)-8{\mathrm{J}}_2(x).\end{array}$$(C.4)

Inserting for ξ_± from Eq. (15) into Eq. (C.1) we arrive at,

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\mathrm{E},l}^{(ij)}& =\frac{1}{2{\mathcal{N}}_l}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\ell \ell [W_{\mathrm{EE}}^l(\ell )C_{\epsilon \epsilon ,\mathrm{E}}^{(ij)}(\ell )+W_{\mathrm{EB}}^l(\ell )C_{\epsilon \epsilon ,\mathrm{B}}^{(ij)}(\ell )],\hfill \\ \hfill {\mathcal{C}}_{\mathrm{B},l}^{(ij)}& =\frac{1}{2{\mathcal{N}}_l}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\ell \ell [W_{\mathrm{BE}}^l(\ell )C_{\epsilon \epsilon ,\mathrm{E}}^{(ij)}(\ell )+W_{\mathrm{BB}}^l(\ell )C_{\epsilon \epsilon ,\mathrm{B}}^{(ij)}(\ell )],\hfill \end{array}$$(C.5)

where the weight functions are,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill W_{\mathrm{EE}}^l(\ell )& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}\theta \theta T(\theta )[{\mathrm{J}}_0(\ell \theta )g_{+}^l(\theta )+{\mathrm{J}}_4(\ell \theta )g_{-}^l(\theta )],\hfill \\ \hfill W_{\mathrm{EB}}^l(\ell )& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}\theta \theta T(\theta )[{\mathrm{J}}_0(\ell \theta )g_{+}^l(\theta )-{\mathrm{J}}_4(\ell \theta )g_{-}^l(\theta )],\hfill \end{array}\end{array}$$(C.6)

with W_EE^l(ℓ) = W_BB^l(ℓ) and W_EB^l(ℓ) = W_BE^l(ℓ). Since ξ_± are measured between θ_min and θ_max, T(θ) has to be defined such that no information beyond this range is included. To avoid ringing, we apodise T(θ) with a Hahn window and smooth its transitions to zero:

$$\begin{array}{c}\hfill T(\theta )=\{\begin{array}{cc}\hfill 0;& x<x_{\mathrm{lo}}-\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hfill {cos}^2\left[\frac{\pi }{2}\frac{x-(x_{\mathrm{lo}}+{\mathrm{\Delta }}_x/2)}{{\mathrm{\Delta }}_x}\right];& x_{\mathrm{lo}}-\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{lo}}+\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hfill 1;& x_{\mathrm{lo}}+\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{up}}-\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hfill {cos}^2\left[\frac{\pi }{2}\frac{x-(x_{\mathrm{up}}-{\mathrm{\Delta }}_x/2)}{{\mathrm{\Delta }}_x}\right];& x_{\mathrm{up}}-\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{up}}+\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hfill 0;& x\ge x_{\mathrm{up}}+\frac{{\mathrm{\Delta }}_x}{2}.\hfill \end{array}\end{array}$$(C.7)

Here x = log(θ), Δ_x is the log-width of the apodisation, x_lo − Δ_x/2 = log(θ_min) and x_up + Δ_x/2 = log(θ_max). We require Δ_x > log(θ_up/θ_lo), and in practice have used Δ_x = 0.5 for all KiDS-Legacy 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_lo = log(θ_min) and x_up = log(θ_max). This change was made as it was deemed more intuitive/appropriate that the band power windows be contained primarily within the requested ℓ ranges after apodisation, rather than before.

Appendix D: 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 (see Reischke et al. 2025) in more detail.

The log-normal simulations are based on GLASS¹⁰ (Generator for Large-Scale Structure, Tessore et al. 2023) with a detailed description given in von Wietersheim-Kramsta et al. (2025) as well as Section 8.1 in Reischke et al. (2025). 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 SALMO¹¹ (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 (Heydenreich et al. 2020), 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_TG, (see Yan et al. 2025, 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. (2025) demonstrated that the TG weights are uncorrelated with 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-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_eff at each pixel.

Fig. D.1. Cartesian spatial map (corresponding Nside = 1024, showing a local flat projection of the healpix map) of the TG weights, wTG, throughout a single field of our KiDS-Legacy-like mock catalogue. Larger TG weights correspond to a higher local source density.

Fig. D.2. Dependence of the galaxy density, neff, on TG weights, wTG in the KiDS-Legacy-like mock catalogue. The data points represent the mean neff of ten equi-populated bins in wTG. 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 neff calculated from the galaxy samples with variable depth per tomographic bin, while the dashed horizontal lines show the values of neff for the respective galaxy samples without any spatial variations in the observational depth.

The local density variations seen here also affect the redshift distributions. This issue is addressed by reapplying the photometric redshift calibration using SOMs (see Wright et al. 2025, 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. (2025) with 10 percent agreement for the variances, very similar to what was found in previous analysis (Joachimi et al. 2021). 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.

Appendix E: B-mode analysis

During the initial phases of the blinded KiDS-Legacy analysis, our B-mode null-test analyses (see Sect. 4.3) failed our required significance threshold (p ≥ 0.01). This led to a dedicated effort to investigate the possible causes of the 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-mode signal.

As KiDS-1000 did not show significant 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-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-mode significance in the same footprint as KiDS-1000, finding that the KiDS-Legacy data in this footprint still exceeded the B-mode significance requirement.

Subsequent investigation of the B-mode origin focussed on PSF contamination (via Paulin-Henriksson statistics; see Sect. 4), one- and two-dimensional 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-mode signal. Finally, we investigated the lowest-level possible cause of the anomalous B-mode signal: systematic failures in the astrometric solution used in the THELI calibration (with Gaia) of r-band images. The global astrometric accuracy of KiDS-Legacy is presented in Wright et al. (2024), 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-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_int) measured the systematic difference between the centroid of all targets extracted from these images. This was achieved by matching all sources (within 2″) extracted from ASTRO-WISE and THELI co-added r-band images. With these matches, we computed the radial separation r between each source’s centroid in the two images and then calculated the median of these separations in 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 the THELIr-band images, which are used for shape measurement. This metric (δ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″). For the ASTRO-WISE static reference, we applied selections on SOURCE EXTRACTOR (Bertin & Arnouts 1996) output variables

$$\begin{array}{c}\hfill \mathtt{FLUX\_RADIUS}<4.0\end{array}$$(E.1)