© The Authors 2025

1. Introduction

In contrast to the continuous background expansion of the Universe, the structure growth in the matter density field is less certain, with uncertainties in theoretical models that vary by 10% or more in the non-linear regime through their dependence on baryon-galaxy feedback; possible deviations from purely cold, stable, and interaction-free dark matter; or, perhaps, modifications in the standard model of gravity (e.g., Laguë et al. 2024; Bucko et al. 2024; Mauland et al. 2024; Ferreira et al. 2024; Salcido et al. 2023; Schneider et al. 2022; Harnois-Déraps et al. 2015; Smith & Markovic 2011; Jing et al. 2006, and references therein). Currently the theory uncertainties are limiting studies of cosmological parameters because they often rely on measures of the cosmic structure, such as the matter power spectrum, P_δ(k, z), at different redshifts, z = :1/a − 1, and spatial (comoving) wave number, k (e.g., Pranjal et al. 2025; García-García et al. 2024; Secco et al. 2022; Li et al. 2023; Asgari et al. 2021; Heymans et al. 2013). To shed more light on the physics imprinted on non-linear scales, direct measurements of a model-free P_δ(k, z) with a minimum of assumptions promise to be a valuable model test. That this is already feasible for z ≲ 1 and 10⁻² ≲ k/(h Mpc⁻¹)≲10 with Stage III galaxy surveys when exploiting the weak gravitational lensing effect, such as with KiDS-1000 (Kuijken et al. 2019, 2015), is shown in this work. This foreshadows exciting applications to future lensing data, such as those by Euclid (Euclid Collaboration: Mellier et al. 2025), the Vera C. Rubin Observatory Legacy Survey of Space and Time (The LSST Dark Energy Science Collaboration 2018), or the Nancy Grace Roman Space Telescope (Spergel et al. 2015).

The coherent distortion of distant galaxy images, just sources hereafter, by the weak shearing of light bundles passing through intervening foreground structure probes the matter power spectrum (for a review, see, e.g., Kilbinger 2015; Schneider 2006; Bartelmann & Schneider 2001). More specifically, and putting negligible higher-order corrections aside (Hilbert et al. 2009), the second-order correlations, ξ_±(θ), at lag θ in the cosmic shear field are simply a linear projection of P_δ(k, z) throughout the look-back light cone up to the distance of the sources (Preston et al. 2024; Mezzetti et al. 2012; Simon 2012; Bacon et al. 2005; Pen et al. 2003; Tegmark & Zaldarriaga 2002; Seljak 1998). Therefore, by using ξ_±^(ij)(θ) between tomographic bins, i and j, for different characteristic distances of sources, the correlation data can be deprojected to recover, within limitations, the original P_δ(k, z). More conveniently, as in Simon (2012, S12), we focus on a transfer function, f_δ(k, z) = P_δ(k, z)/P_δ^fid(k, z), with respect to a reference model, P_δ^fid(k, z), to highlight deviations from the presumed reference inside an average z-bin. For this reference, we chose a purely cold dark matter (CDM) model to probe for deviations, f_δ(k, z)≠1, that may hint at missing physics in a basic, well-understood dark matter scenario or at unidentified systematic errors in the data. However, at the end of the day, the kind of reference is a deliberate choice.

Required for the deprojection, on the other hand, is the knowledge of the projection (lensing) kernel and the inclusion of the intrinsic alignment (IA) of sources in the correlation function, an essential part of modern cosmic-shear analyses (Hirata & Seljak 2004; Crittenden et al. 2001; Croft & Metzler 2000). The lensing kernel is fully described by the radial distribution of sources inside the tomographic bins; the average matter density in today’s Universe, Ω_m, relative to the critical density, ρ_crit = 3H₀²(8πG_N)⁻¹; and the expansion rate $E(a):=H_0^{-1}\dot{a}/a$, where H₀ = 100 h km s⁻¹ Mpc⁻¹ is the Hubble constant and G_N is Newton’s gravitational constant. The kernel parameters Ω_m and E(a) are confined to percentage precision by cosmological experiments already (E²(a)≈Ω_m[a⁻³ − 1]+1 for a flat universe), especially owing to CMB experiments (Planck Collaboration VI 2020; Hinshaw et al. 2013), but also by combining probes from the closer Universe (e.g., Heymans et al. 2021; Abbott et al. 2018). Nevertheless, in the refinement of the method in S12 for KiDS-1000, we marginalised here over the small background cosmology uncertainties, and over those in the source distributions (Hildebrandt et al. 2021). Another refinement is the required inclusion of the IA, and its uncertainty, in the deprojection procedure by the widely employed non-linear alignment (NLA) model (Joachimi et al. 2011; Bridle & King 2007), informed by the IA constraints in Asgari et al. (2021).

Despite the simplicity of the well-defined projection, the noise level in ξ_±^(ij)(θ) makes the recovery of the matter power spectrum still challenging because the deprojection has to undo a convolution in radial and transverse direction, producing strongly correlated, oscillating noise in an unstable reconstruction. This complication could be mitigated by data volumes substantially larger than KiDS-1000, beating down oscillating noise, or by additional assumptions on the redshift dependence or shape of the power spectrum in fitting an analytical model with few parameters to the data (Perez Sarmiento et al. 2025; Ye et al. 2024; Truttero et al. 2025; Broxterman & Kuijken 2024; Preston et al. 2024, 2023; Pen et al. 2003; Seljak 1998). Avoiding analytical models, we instead propose constraining f_δ(k, z) averages in different k- and z-bins, as in S12 but additionally filtered by a classic Tikhonov regularisation, for instance as applied in the mathematically related problem of deconvolving noisy images (e.g., Murata & Takeuchi 2022). In contrast to S12, the Tikhonov regularisation, installed on our Bayesian statistical model as prior, penalises strongly oscillating power spectra for noisy data, giving an advantage to solutions of f_δ(k, z) that, on the one hand, are smoother in the k-direction, but, on the other hand, use no prior information on their z-dependence. We show with verification data that this regularisation indeed stabilises the reconstruction and extracts from KiDS-1000 data useful constraints on the transfer function f_δ(k, z), either averaged over the full redshift range or averaged within three separate redshift bins to broadly probe a redshift evolution. Furthermore, for an efficient sampling of the posterior constraints of the matter power spectrum, we describe and test a Markov chain Monte Carlo (MCMC) code tailored to deal with the numerous 60 (or more) degrees-of-freedom of the binned, band power P_δ(k, z); its practical details are given in Appendix A.

The outline of the paper is as follows. Section 2 reviews our weak lensing formalism and shows that the projection of P_δ(k, z) into ξ_±^(ij)(θ), both binned, can be reduced to one projection matrix, to be reused as long as the same lensing kernel and IA parameters, including source redshift distributions, are employed. Section 3 defines the statistical model for our Bayesian analysis, including a prior for the Tikhonov regularisation. We summarise the KiDS-1000 data for our cosmic shear analysis in Sect. 4 and report our results in Sect. 5. Section 5 also reports, as verification of our reconstruction method, the results of a mock analysis based on KiDS-1000-like data produced by ray-tracing N-body data. Another verification test, this time based on the data vector of an analytical model subject to random noise, is presented in Appendix A. We discuss our results and conclusions on the reconstructed matter power spectrum in Sect. 6. Notably, our figures use comoving wave numbers, k, for spatial scales of the matter power spectrum, emphasised by the ‘c’ in the unit [k]=h cMpc⁻¹.

2. Weak lensing formalism

This work is an application of the well-established theory of cosmic shear. For a review and its mathematical foundation, we refer to Kilbinger (2015) or Schneider (2006), and only briefly summarise our formalism here.

2.1. Cosmic shear

Weak gravitational lensing by fluctuations in the large-scale structure of the foreground matter density, ${\delta }_{\mathrm{m}}={\rho }_{\mathrm{m}}/{\overline{\rho }}_{\mathrm{m}}-1$, distorts images of background galaxies. Their average shear distortion in direction θ of the sky shall be expressed by the complex γ(θ) = γ₁(θ)+iγ₂(θ) for an ensemble of galaxies with probability distribution function (PDF) p_χ(χ) in comoving distance, χ. To lowest order, and sufficiently accurate in practical applications, the Fourier transform $\stackrel{\sim }{\gamma }(\mathit{\ell })$ is related to the linear projection of fluctuations along θ,

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}\frac{\mathrm{d}\chi }{a(\chi )}\overline{W}(\chi )f_{\mathrm{K}}(\chi ){\delta }_{\mathrm{m}}[f_{\mathrm{K}}(\chi )\mathit{\theta },\chi ],\end{array}$$(1)

through the convolution (ℓ ≠ 0)

$$\begin{array}{c}\hfill \stackrel{\sim }{\gamma }(\mathit{\ell })={\displaystyle {\int }_{{\mathbb{R}}^2}}{\mathrm{d}}^2\theta \gamma (\mathit{\theta })exp(-\mathrm{i}\mathit{\theta }·\mathit{\ell })=\frac{{({\ell }_1+\mathrm{i}{\ell }_2)}^2}{{\ell }^2}\stackrel{\sim }{\kappa }(\mathit{\ell }),\end{array}$$(2)

where c is the vacuum speed of light. In the above equations, $\stackrel{\sim }{\kappa }(\mathit{\ell })$ is the Fourier transform of the convergence κ(θ), the scalar f_K(χ) is the (comoving) angular diameter distance for the curvature scalar K, the vector f_K(χ)θ denotes a separation in a tangential plane on the sky at distance χ, the term 1 + z = a⁻¹(χ) is our relation between redshift, z, and scale factor, a(χ), at the look back time of χ, and

$$\begin{array}{c}\hfill \overline{W}(\chi ):={\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{h}}}}\mathrm{d}{\chi }^{\prime }{p}_{\chi }({\chi }^{\prime })\frac{f_{\mathrm{K}}({\chi }^{\prime }-\chi )}{f_{\mathrm{K}}({\chi }^{\prime })}\end{array}$$(3)

is the lensing efficiency, cut off at the size of the observable Universe, χ_h. As approximation, we assumed a flat sky with Cartesian coordinates for θ but expect negligible inaccuracies for angular separations below several degrees.

Practical estimators of γ(θ) in the direction θ of a galaxy image use a convenient definition of galaxy ellipticity, ϵ, that is calibrated to be an unbiased estimator, ⟨ϵ⟩=γ, for a selected source population (e.g., Giblin et al. 2021).

2.2. Second-order statistics

Our analysis exploits the coherent shear distortion of galaxy images over the sky to infer the three-dimensional power spectrum, P_δ(k, z), of matter density fluctuations inside the light cone,

$$\begin{array}{c}\hfill ⟨{\stackrel{\sim }{\delta }}_{\mathrm{m}}(\mathit{k},z){\stackrel{\sim }{\delta }}_{\mathrm{m}}^{\ast }({\mathit{k}}^{\prime },z)⟩=:{(2\pi )}^3{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}^{\prime })P_{\delta }(k,z),\end{array}$$(4)

for a range of wave numbers and redshifts. The δ_D(x) denotes the Dirac delta function, and ${\stackrel{\sim }{\delta }}_{\mathrm{m}}(\mathit{k},z)$ is the Fourier coefficient of a fluctuation mode at redshift z inside the light cone. To obtain the power spectrum from γ(θ), we estimate the two-point correlation function of cosmic shear,

$$\begin{array}{c}\hfill {\xi }_{\pm }(\theta ):=\langle {\gamma }_{\mathrm{t}}({\mathit{\theta }}_1){\gamma }_{\mathrm{t}}({\mathit{\theta }}_2)\rangle \pm \langle {\gamma }_{\times }({\mathit{\theta }}_1){\gamma }_{\times }({\mathit{\theta }}_2)\rangle \end{array}$$(5)

at lag θ = |θ₂ − θ₁| from an ensemble of sources. Here, γ_t and γ_× denote the tangential and cross-shear components relative to the orientation of θ₂ − θ₁ = :θ e^iϕ on the sky, defined by γ_t + iγ_× := −e^−2iϕγ. In addition, using the hybrid extended Limber approximation in Fourier space (Kilbinger et al. 2017; Kaiser 1992), the correlation function of cosmic shear is a linear projection of P_δ(k, z),

$$\begin{array}{cc}& {\xi }_{\pm }(\theta )=\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\chi \mathrm{d}\ell \ell }{2\pi }\frac{{\overline{W}}^2(\chi )}{a^2(\chi )}{J}_{0,4}(\ell \theta )P_{\delta }(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},\chi )\hfill \end{array}$$(6)

through the lensing kernel $K(\chi )\propto {\mathrm{\Omega }}_{\mathrm{m}}^2{\overline{W}}^2(\chi )a^{-2}(\chi )$ – the minimal ingredient that has to be known to invert the projection. The expressions J_n(x) denote n^th-order Bessel functions of the first kind, of which J₀(x) has to be applied for ξ₊ and J₄(x) for ξ₋.

2.3. Tomographic lensing and intrinsic alignment of sources

Owing to the information loss in the projection Eq. (6) for a single sample of sources, a power spectrum averaged over the entire light cone could be obtained at best, as in Pen et al. (2003) or Schneider et al. (2002). More advanced is splitting the source samples by redshift in a tomographic analysis for a better statistical precision, reflected by improved cosmology constraints in tomographic analyses, or to achieve a (limited) redshift resolution for the recovered P_δ(k, z). Therefore, we split our source sample into five redshift bins, using non-overlapping ranges of photometric redshifts, namely the z_B-intervals (0.1, 0.3], (0.3, 0.5], (0.5, 0.7], (0.7, 0.9], and (0.9, 1.2], as done already in Hildebrandt et al. (2021). Figure 1 plots estimates of the resulting distributions p_z⁽ⁱ⁾(z) for each bin i. The tomographic analysis then correlates the shear signal between bins i and j, denoted by the superscript “(ij)” in ξ_±^(ij)(θ). Compared to Eq. (6), this increases the data vector size by a factor of 5(5 + 1)/2 = 15 but also adds information on the P_δ(k, z) evolution because the shear tomography probes the same foreground with sources at different z.

Fig. 1. Probability density distribution functions, pz(i)(z), of KiDS-1000 source galaxies within the five tomographic redshift bins (from i = 1 to i = 5): (0.1, 0.3], (0.3, 0.5], (0.5, 0.7], (0.7, 0.9], and (0.9, 1.2]. These estimates are from Hildebrandt et al. (2021).

This tomographic analysis requires special attention with regard to additional contributions to the shear signal from the intrinsic alignment (IA) of sources (e.g., Lamman et al. 2024). In the ideal absence of IA, the orientations of intrinsic source ellipticities are statistically independent among each other and to the cosmic shear signal. So-called ‘II’ contributions, however, originate from correlated orientations of physically close galaxies. In addition, galaxy shapes are aligned to the surrounding matter density, giving rise to a ‘GI’ signal (Hirata & Seljak 2004). While II contributions could in principle be reduced by correlating only bin combinations with little radial overlap, and the GI signal could be suppressed by a nulling technique (Joachimi & Schneider 2008), we followed the common approach of modelling the II and GI signal in ξ_±^(ij)(θ) to avoid information loss,

$$\begin{array}{c}\hfill {\xi }_{\pm }^{(ij)}(\theta ){=}^0{\xi }_{\pm }^{(ij)}(\theta ){+}^{\mathrm{I}\mathrm{I}}{\xi }_{\pm }^{(ij)}(\theta ){+}^{\mathrm{G}\mathrm{I}}{\xi }_{\pm }^{(ij)}(\theta ).\end{array}$$(7)

In this equation,

$$\begin{array}{cc}& {}^0{\xi }_{\pm }^{(ij)}(\theta ):=\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\chi \mathrm{d}\ell \ell }{2\pi }\frac{{\overline{W}}^{(i)}(\chi ){\overline{W}}^{(j)}(\chi )}{a^2(\chi )}{J}_{0,4}(\ell \theta )P_{\delta }(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},\chi )\hfill \end{array}$$(8)

is the ideal shear signal without IA. Covering both auto- (i = j) and cross-correlations (i ≠ j) of source samples, this expression is more general than Eq. (6). Herein, the lensing efficiency with superscript, ${\overline{W}}^{(i)}(\chi )$, refers to Eq. (3) but for the PDF p_χ⁽ⁱ⁾(χ) of the ith bin. Notably, the p_χ⁽ⁱ⁾(χ) are related to the p_z⁽ⁱ⁾(z) in Fig. 1 by p_χ⁽ⁱ⁾(χ) = p_z⁽ⁱ⁾[z(χ)] |dz/dχ| for dz/dχ = H₀ c⁻¹ E[(1 + z)⁻¹].

For the II and GI terms, we employed the NLA model by Joachimi et al. (2011), derived from the linear alignment model by Bridle & King (2007) but replacing the linear matter power-spectrum by the non-linear one. Clearly just tweaking the original linear IA model, the NLA is nevertheless sufficiently accurate to model contemporary lensing data (Harnois-Déraps et al. 2022), predicting

$$\begin{array}{cc}\hfill {}^{\mathrm{I}\mathrm{I}}{\xi }_{\pm }^{(ij)}(\theta )& :={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\chi \mathrm{d}\ell \ell }{2\pi }{F}^2(\chi )\hfill \\ \hfill & \times \frac{p_{\chi }^{(i)}(\chi )p_{\chi }^{(j)}(\chi )}{f_{\mathrm{K}}^2(\chi )}{J}_{0,4}(\ell \theta )P_{\delta }(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},\chi )\hfill \end{array}$$(9)

for II correlations and for the GI term

$$\begin{array}{cc}\hfill {}^{\mathrm{G}\mathrm{I}}{\xi }_{\pm }^{(ij)}(\theta )& :=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\chi \mathrm{d}\ell \ell }{2\pi }F(\chi )\hfill \\ \hfill & \times \frac{{\overline{W}}^{(i)}(\chi )p_{\chi }^{(j)}(\chi )+{\overline{W}}^{(j)}(\chi )p_{\chi }^{(i)}(\chi )}{a(\chi )f_{\mathrm{K}}(\chi )}{J}_{0,4}(\ell \theta )P_{\delta }(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},\chi ).\hfill \end{array}$$(10)

The amplitude of the IA signal scales in this model with

$$\begin{array}{cc}\hfill F(\chi )& =-A_{\mathrm{IA}}{C}_1{\rho }_{\mathrm{crit}}\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{D_{+}(\chi )}\hfill \\ \hfill & \approx -2.4\times {10}^{-2}\left(\frac{A_{\mathrm{IA}}}{3.0}\right)\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{0.3}\right){\left(\frac{D_{+}(\chi )}{0.5}\right)}^{-1}\hfill \end{array}$$(11)

and depends on distance χ (or redshift) only through the linear growth factor, D₊(χ) (by definition D₊ ≡ 1 at χ = 0). For our main result, we neglected further dependencies on redshift or the evolution of the average galaxy luminosity with redshift, similar to Asgari et al. (2021) who find in their cosmological analysis of KiDS-1000 little evidence for a more complex F(χ).

With respect to future applications, the work by Fortuna et al. (2021), and, more recently, by Preston et al. (2024) observe that the NLA model, even if overly simplistic, could be accurate enough to model IA in Stage IV survey data if a redshift dependence of the NLA parameters is accounted for. We briefly return to this topic in our discussion on conceivable systematic uncertainties in our analysis in Sect. 6.

2.4. Projection kernels of the matter power spectrum

The shear correlation function with IA terms, Eq. (7), is still linear in the matter power spectrum, hence a deprojection of ξ_±^(ij)(θ) into P_δ(k, z) is, as in S12, principally possible through linear minimum-variance estimators that depend only on the background fiducial cosmology, source distributions, and IA parameters. On the practical side, however, the deprojection is hampered by broad, partly similar lensing kernels, the relatively low signal-to-noise ratio, and a confined θ-range, rendering the estimators ill-conditioned (or singular). They are also biased when ignoring the θ-confinement which amounts to setting ξ^(ij)(θ)≡0 outside the θ range (Section 7 in Schneider et al. 2002). Similar to S12, we addressed these practicalities by boundary conditions P_δ(k, z)≥0 through priors in the framework of a Bayesian analysis, at most a couple of redshift bins for P_δ(k, z), and by varying P_δ(k, z) within a confined region in (k, z)-space only, assuming a reference power, P_δ^fid(k, z), otherwise. Additionally, it is sensible to constrain the relative deviations f_δ(k, z):=P_δ(k, z)/P_δ^fid(k, z) instead of P_δ(k, z) directly: Most of the evolution is probably already accounted for in P_δ^fid(k, z), and a slowly changing f_δ(k, z) within a broad z-bin is a reasonable quantity to be averaged. We describe the implementation details below.

In our set-up, the f_δ(k, z) shall be constant within cells of a regular mesh of (N_k + 1)×(N_z + 1) mesh points (k_m, χ_n), where k_m < k_m + 1, z_n < z_n + 1, and χ_n := χ(z_n); the intervals [z₁, z_Nz + 1) and [k₁, k_Nk + 1) define the z- and k-ranges where the power spectrum may differ from P_δ^fid(k, z), here for $k/(h{\mathrm{Mpc}}^{-1})\in [{10}^{-2},20)$ and z ∈ [0, 2). The matter power spectrum for Eq. (7) thus equals

$$\begin{array}{cc}\hfill P_{\delta }(k,\chi )& =P_{\delta }^{\mathrm{fid}}(k,\chi )(1+{\displaystyle \sum _{n,m=1}^{N_z,N_k}}{H}_{\mathit{mn}}(k,\chi )[f_{\delta ,mn}-1]),\hfill \end{array}$$(12)

where

$$\begin{array}{c}\hfill H_{\mathit{mn}}(k,\chi ):=\{\begin{array}{cc}1\hfill & ,\mathrm{if}k\in [k_m,k_{m+1})\mathrm{and}\chi \in [{\chi }_n,{\chi }_{n+1})\hfill \\ 0\hfill & ,\mathrm{otherwise}\hfill \end{array}\end{array}$$(13)

is a two-dimensional top-hat function. The average of f_δ(k, z) within a cell (or band) is denoted by the coefficient f_δ, mn. The mesh in k-direction is equi-spaced on a log-scale, Δk = k N_k⁻¹ ln(k_Nk + 1/k₁), while the z-direction uses for one variant z_n ∈ {0, 0.3, 0.6, 2}, for N_z = 3, and z_n ∈ {0, 2} in other variant with one wide bin, N_z = 1, when averaging f_δ(k, z) over the entire redshift range. The number of k-bins is always N_k = 20.

Using this f_δ, mn-representation of P_δ(k, z), Eq. (7) is now recast into

$$\begin{array}{c}\hfill {\xi }_{\pm }^{(ij)}(\theta )={\displaystyle \sum _{n,m=1}^{N_z,N_k}}{X}_{\pm }^{(ij)}(\theta ;m,n)f_{\delta ,mn}+{\xi }_{\pm ,\mathrm{fid}}^{(ij)}(\theta ),\end{array}$$(14)

for the projection matrix

$$\begin{array}{cc}\hfill X_{\pm }^{(ij)}(\theta ;m,n)& :=\frac{1}{2\pi {\theta }^2}{\displaystyle {\int }_{{\chi }_n}^{{\chi }_{n+1}}}\mathrm{d}\chi \left({}^0{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi ){+}^{\mathrm{II}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi ){+}^{\mathrm{GI}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi )\right)\hfill \\ \hfill & \times {\displaystyle \underset{k_m{f}_{\mathrm{K}}(\chi )\theta }{\overset{k_{m+1}{f}_{\mathrm{K}}(\chi )\theta }{\int }}}\mathrm{d}ss{J}_{0,4}(s)P_{\delta }^{\mathrm{fid}}(k[s,\chi ,\theta ],\chi ),\hfill \end{array}$$(15)

the integral kernels

$$\begin{array}{cc}\hfill {}^0{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi )& :=\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}\frac{{\overline{W}}^{(i)}(\chi ){\overline{W}}^{(j)}(\chi )}{a^2(\chi )};\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {}^{\mathrm{II}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi )& :=F^2(\chi )\frac{p_{\chi }^{(i)}(\chi )p_{\chi }^{(j)}(\chi )}{f_{\mathrm{K}}^2(\chi )};\hfill \\ \hfill {}^{\mathrm{GI}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi )& :=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}F(\chi )\hfill \end{array}$$(17)

$$\begin{array}{cc}& \times \frac{{\overline{W}}^{(i)}(\chi )p_{\chi }^{(j)}(\chi )+{\overline{W}}^{(j)}(\chi )p_{\chi }^{(i)}(\chi )}{a(\chi )f_{\mathrm{K}}(\chi )},\hfill \end{array}$$(18)

and a constant offset

$$\begin{array}{cc}\hfill {\xi }_{\pm ,\mathrm{fid}}^{(ij)}(\theta )& :=\frac{1}{2\pi {\theta }^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}}\mathrm{d}\chi \left({}^0{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi ){+}^{\mathrm{II}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi ){+}^{\mathrm{GI}}{\mathcal{K}}_{\gamma \gamma }^{(ij)}(\chi )\right)\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{\infty }{\int }}}\mathrm{d}ss{J}_{0,4}(s)P_{\delta }^{\mathrm{fid}}(k[s,\chi ,\theta ],\chi )\hfill \end{array}$$(19)

from (k, z) regions that are unaffected by the choice of f_δ, mn. The offset employs the definition

$$\begin{array}{c}\hfill P_{\delta }^{\mathrm{fid}}(k,\chi ):=\{\begin{array}{cc}{P}_{\delta }^{\mathrm{fid}}(k,\chi )\hfill & ,\mathrm{if}{\displaystyle \sum _{m,n=1}^{N_z,N_k}}{H}_{\mathit{mn}}(k,\chi )=0\hfill \\ 0\hfill & ,\mathrm{otherwise}\hfill \end{array},\end{array}$$(20)

and the foregoing equations abbreviate $k[s,\chi ,\theta ]:=\frac{s+\theta /2}{f_{\mathrm{K}}(\chi )\theta }$.

Our code computed the projection matrix X_±^(ij)(θ; m, n) once for a series of θ-bins, enabling a quick prediction of ξ_±^(ij)(θ) when Monte Carlo sampling the posterior PDF of f_δ, mn. An efficient way to numerically calculate X_±^(ij)(θ; m, n) and ${\xi }_{\pm ,\mathrm{fid}}^{(ij)}(\theta )$ is given in Appendix A of S12, albeit with the little modification (2f_K[χ])⁻¹ in what is now ${\overline{k}}_j:={(2f_{\mathrm{K}}[\chi ])}^{-1}+({\widehat{k}}_j+{\widehat{k}}_{j+1})/2$ due to the extended hybrid Limber approximation adopted here. As further practical footnote on code implementation, we performed the above integrals with integration variable a instead of χ, as well as p_z(z) instead of p_χ(χ), for which the following transformations are notable, given a general function g(χ) and H(a) =  : H₀ E(a),

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\chi }_1}^{{\chi }_2}}\mathrm{d}\chi g(\chi )=\frac{c}{H_0}{\displaystyle {\int }_{a_2}^{a_1}}\frac{\mathrm{d}a}{a^{2}E(a)}g[\chi (a)]\end{array}$$(21)

since |dχ/da|=c/(a²H[a]),

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\chi }_1}^{{\chi }_2}}\mathrm{d}\chi p_{\chi }^{(i)}(\chi )g(\chi )={\displaystyle {\int }_{a_2}^{a_1}}\frac{\mathrm{d}a}{a^2}{p}_z^{(i)}[z(a)]g[\chi (a)]\end{array}$$(22)

since |dz/dχ|=H(z)/c, and hence

$$\begin{array}{cc}& {\displaystyle {\int }_{{\chi }_1}^{{\chi }_2}}\mathrm{d}\chi p_{\chi }^{(i)}(\chi )p_{\chi }^{(j)}(\chi )g(\chi )\hfill \\ \hfill & =\frac{H_0}{c}{\displaystyle {\int }_{a_2}^{a_1}}\frac{\mathrm{d}aE(a)}{a^2}{p}_z^{(i)}[z(a)]p_z^{(j)}[z(a)]g[\chi (a)].\hfill \end{array}$$(23)

When integrating over p_z⁽ⁱ⁾(z), we adopted histograms with top-hat functions, without interpolation, exactly as indicated in Fig. 1. The very first bin, left-aligned towards z = 0, has p_z⁽ⁱ⁾(z) = 0. Probably obvious through these integrals and the pre-factors in Eqs. (8)–(10), the projections of P_δ(k, z) into ξ^(ij)(θ) are independent of H₀ as long as [k]=h Mpc⁻¹ (comoving), whereas the ξ_±^(ij)(θ) scale with Ω_m² since F(χ)∝Ω_m.

3. Bayesian inference of the three-dimensional matter power-spectrum

With the linear relation between the three-dimensional (3D) power spectrum and the two-point shear correlation functions described, we now outline the statistical model and the numerical sampler for the posterior probability density of f_δ, mn within a Bayesian framework (e.g., Gelman et al. 2003, for a review). In comparison to S12, our updated methodology introduces the Tikhonov regularisation as Bayesian prior, in effect inferring a k-smoothed power spectrum, and a Hamiltonian MCMC for improved computational speed for N_z × N_k ∼ 60 model variables.

3.1. Statistical model

The statistical model of f_δ, mn employs a compact notation for the tomographic data, model variables and parameters, all presented here. The vector of model variables π = (f_δ, mn|1≤n≤N_z,1≤m≤N_k) compiles the bin-averaged f_δ(k, z) on a (N_k + 1)×(N_z + 1) mesh, and the data vector d contains the binned ξ_±^(ij)(θ_k) for a sequence of N_θ = 9 angular separations, θ_k. Equation (14) predicts the components of d for a given (flat) background cosmology, IA parameters, and source redshift distributions, altogether compressed into q = {Ω_m, A_IA, p_z⁽¹⁾(z),…,p_z⁽⁵⁾(z)} hereafter. The prediction is denoted by the model vector $\mathit{m}(\mathit{\pi },\mathit{q})={\mathsf{X}}_{\mathit{q}}\mathit{\pi }+{\mathit{\xi }}_{\mathit{q}}^{\mathrm{fid}}$ with elements arranged in the same order as those in d, using the projection matrix X_q which contains the coefficients X_±^(ij)(θ; m, n) for fixed q. The elements of ${\mathit{\xi }}_{\mathit{q}}^{\mathrm{fid}}$ are the offsets ${\xi }_{\pm ,\mathrm{fid}}^{(ij)}({\theta }_k)$, also for fixed projection parameters q. The statistical information on π is expressed by the Bayesian posterior PDF of π conditional on d,

$$\begin{array}{c}\hfill P(\mathit{\pi }|\mathit{d})={\mathcal{E}}^{-1}(\mathit{d}){\mathcal{L}}_{\mathit{q}}(\mathit{d}|\mathit{\pi })P_{\mathrm{hat}}(\mathit{\pi })P_{\tau }(\mathit{\pi }),\end{array}$$(24)

adopting a Gaussian likelihood

$$\begin{array}{c}\hfill -2ln{\mathcal{L}}_{\mathit{q}}(\mathit{d}|\mathit{\pi })={[\mathit{d}-\mathit{m}(\mathit{\pi },\mathit{q})]}^{\mathrm{T}}{\mathsf{C}}^{-1}[\mathit{d}-\mathit{m}(\mathit{\pi },\mathit{q})]\end{array}$$(25)

with error covariance C of the measurement d, and two prior probability densities, P_hat(π) and P_τ(π), to regularise the sampling in f_δ-parameter space. The normalisation, $\mathcal{E}(\mathit{d})={\displaystyle \int \mathrm{d}}\mathit{\pi }{\mathcal{L}}_{\mathit{q}}(\mathit{d}|\mathit{\pi })P_{\mathrm{hat}}(\mathit{\pi })P_{\tau }(\mathit{\pi })\ne 0$, is not of interest for the Monte Carlo sampler, and hence we set $\mathcal{E}(\mathit{d})\equiv 1$ without impacting the results. We describe the two prior PDFs in the following.

3.2. Tikhonov regularisation prior

The Tikhonov regularisation prior, P_τ(π), greatly improves the constraints, pruning implausible solutions for f_δ, mn. This is necessary because the combination of heavy smoothing and added noise in ξ_±^(ij)(θ) has the unpleasant side-effect of producing strong oscillations, correlated errors, in f_δ(k, z) when inverting the tomographic signal. Similar to 3D-lensing mass reconstructions (Hu & Keeton 2002), the oscillations are reduced by applying regularisation conditions. To this effect, the Tikhonov regularisation down-weights oscillating f_δ(k, z), preferring solutions of f_δ, mn smoothed in k-direction,

$$\begin{array}{c}\hfill -\frac{1}{\tau }ln{P}_{\tau }(\mathit{\pi })={\displaystyle \sum _{n=1}^{N_z}}{\displaystyle \sum _{m=1}^{N_k-1}}{(f_{\delta ,mn}-f_{\delta ,(m+1)n})}^2.\end{array}$$(26)

At the same time, no prior constraints are imposed on the overall amplitude of f_δ, mn or the difference signal between z-bins because only signal differences from the same z-bin, n, appear inside the prior density. The weight of the prior relative to the data likelihood ${\mathcal{L}}_{\mathit{q}}(\mathit{d}|\mathit{\pi })$ is controlled by the Tikhonov parameter τ, demonstrated in Fig. 2 for a simulated analysis with three redshift bins, N_z = 3.

Fig. 2. Impact of the Tikhonov regularisation used to suppress oscillating solutions fδ(k, z). Shown are, for a noise-free mock data vector and the KiDS-1000 error covariance, the posterior constraints (68% credible regions) on fδ averaged over the redshift bin Z1 = [0, 0.3] (left), Z2 = [0.3, 0.6] (middle), and Z3 = [0.6, 2] (right) with and without regularisation (dark orange τ = 5.0 with median as dashed line or light orange τ = 0). To boost scale-dependence and evolution, fδ(k, z) is here defined relative to the power spectrum at fixed redshift, probing the relative structure growth since z = 1. The solid line is the median posterior fδ for 100× reduced measurement errors, providing a nearly noise-free reference (noise/100) that averages the growth over the redshift bin while still exhibiting artefacts near the edges. The dotted lines are the theoretical $P_{\delta }^{\mathrm{fid}}(k,\overline{z})/P_{\delta }^{\mathrm{fid}}(k,z=1)$ for one specific $\overline{z}$ chosen to most closely match the solid lines, indicating the redshift with highest weight in the average.

This intentionally extreme scenario in Fig. 2 highlights the impact of the regularisation with a truly z-evolving and scale-dependent f_δ, mn by choosing a constant reference P_δ^fid(k, z) for all z-bins: a theory power spectrum at z = 1. Unlike the actual KiDS analysis in Sect. 5, the f_δ, mn therefore now quantify the structure growth relative to z = 1. To assess the quality of the reconstruction and our ability to detect the z- and k-dependence of f_δ, mn, we analysed a mock vector d = m(π, q) with parameters in Table 1 and KiDS-1000-like noise, C, yielding posterior 68% credible regions with regularisation (τ = 5, dark orange) and without (τ = 0, light orange). The figure compares these credible regions to virtually noise-free data, unaffected by the Tikhonov regularisation, depicted as solid lines (“noise/100”). The solid lines are, to account for the z-weighting of f_δ(k, z) in a reconstruction, the posterior medians for a data vector with unrealistically high signal-to-noise ratio (S/N) of C′ = 10⁻⁴ × C as error covariance and no Tikhonov regularisation. But even for such high S/N data, artefacts appear towards large scales, close to k ∼ 0.02 h Mpc⁻¹, and at the high-k end, near k ∼ 10 h Mpc⁻¹, thus towards k where the tomographic data poorly constrain the power spectrum. To mitigate these artefacts, the additional dotted lines depict for each z-bin a theoretical $P_{\delta }^{\mathrm{fid}}(k,\overline{z})/P_{\delta }^{\mathrm{fid}}(k,z=1)$ at one specific redshift, $\overline{z}$, that matches the solid lines most closely, and which shall be the ‘true’ average f_δ, mn for this experiment. Undoubtedly, a Tikhonov prior, τ = 5, substantially reduces the statistical uncertainty in the reconstruction (sizes of credible regions for τ = 0 versus τ = 5). On the downside, however, the regularisation might be overly restrictive by straightening out oscillations actually present in the data, as already partly present for τ = 5 (dotted lines versus dashed lines, the posterior median). By running more tests with different τ, not shown here, we consider τ = 5 a good compromise that catches gradual trends with k while improving the precision in the reconstruction, as in the left panel. We reiterate that the degree of evolution and scale-dependence in this test, due to the static reference power spectrum, is extreme compared to what is expected in the KiDS-1000 analysis. If the reference is close to the true matter power at all z, f_δ, mn will be close to unity.

A cruder alternative to Tikhonov filtering (τ = 0) for smoothing the reconstruction involves drastically reducing the number of k-bins, such as to N_k = 5. This approach, however, significantly blurs f_δ(k, z) and their scale-dependent features, effectively smoothing it to the fixed size of the larger k-bins, resulting in a distinctly lower resolution compared to Tikhonov filtering. In contrast, Tikhonov filtering applies adaptive smoothing, using smaller smoothing kernels where the S/N is higher and broader kernels elsewhere. This method’s resolution is limited only by the size of the numerous, smaller k-bins, for which we chose N_k = 20.

3.3. Positivity priors

Already applied for Fig. 2, we restricted valid solutions to positive f_δ, mn (since P_δ(k, z)≥0), by using uniform, top-hat prior PDFs for $f_{\delta ,mn}\in [0,f_{\delta ,\max }]$,

$$\begin{array}{cc}& -{\sigma }_{\mathrm{f}}^{2}ln{P}_{\mathrm{hat}}(\mathit{\pi })={\displaystyle \sum _{n,m=1}^{N_z,N_k}}(f_{\delta ,mn}^2\mathrm{H}[-f_{\delta ,mn}]\hfill \\ \hfill & +{[f_{\delta ,mn}-f_{\delta ,\max }]}^2\mathrm{H}[f_{\delta ,mn}-f_{\delta ,\max }]),\hfill \end{array}$$(27)

where H(x) is the Heaviside step function. Enforcing positive solutions for the matter power spectrum is beneficial in reducing both the statistical errors and oscillations in the reconstruction, as already reported in S12. The positivity priors adopt wide intervals, $f_{\delta ,\max }=100$, compared to expected values of f_δ, mn ∼ 1. In addition, for numerical convenience, the top-hat priors have soft edges, σ_f = 10⁻², to relax issues with low acceptance rates and undefined gradients in the Monte Carlo sampler of P(π|d) near the edges.

3.4. Monte Carlo sampler

For sampling the posterior PDF, Eq. (24), we applied a MCMC technique that represents P(π|d) by n_mcmc = 10⁴ points (π_i, w_i) with statistical weights, w_i, for given projection parameters q. The MCMC sampler is a Hamiltonian Monte Carlo algorithm, known to be efficient even for high-dimensional models (N_z × N_k ∼ 60) by proposing new, quickly decorrelating sampling points with high acceptance rate (Gelman et al. 2003). The sampler, however, requires as input the gradient ${\mathrm{\nabla }}_{\mathit{\pi }}ln{\mathcal{L}}_{\mathit{q}}(\mathit{d}|\mathit{\pi })$ and the gradient of the logarithmic prior densities – all easily available in our case, foremost because of the linear projection $\mathit{m}(\mathit{\pi },\mathit{q})={\mathsf{X}}_{\mathit{q}}\mathit{\pi }+{\mathit{\xi }}_{\mathit{q}}^{\mathrm{fid}}$. We demonstrate the appropriate sampler performance in Appendix A and provide more implementation details there. Due to the efficient decorrelation of the MCMCs, we shortened the burn-in phase for the KiDS-1000 analysis, n_burnin = 5 × 10³, in comparison to Appendix A, and used shorter chain lengths because of the subsequent merger of many chains when marginalising over q, outlined in the next section.

3.5. Marginal errors including uncertainties in IA and lensing kernel parameters

The parameters for the lensing kernel and IA are neither exactly known nor intended to be constrained by our shear data. Therefore, to account for their uncertainties in f_δ, mn, the analysis combines many chains with n_runs = 500 realisations of q drawn from a baseline error model. The error baseline adopted a flat cosmology, Ω_Λ = 1 − Ω_m, and a random sample of the posterior joint distribution of (Ω_m, A_IA) from the recent 3 × 2 pt cosmological analysis in Heymans et al. (2021); their marginal errors are ${\mathrm{\Omega }}_{\mathrm{m}}=0.{305}_{-0.012}^{+0.012}$ and $A_{\mathrm{IA}}=+1.{04}_{-0.30}^{+0.28}$ (68% CI). The 3 × 2 pt experiment combines the KiDS-1000 cosmic shear constraints with those from galaxy-galaxy lensing and, importantly, the galaxy clustering in the partly overlapping, spectroscopic surveys of both the 2-degree Field Lensing Survey and the Baryon Oscillation Spectroscopic Survey (Blake et al. 2016; Alam et al. 2015) to break the σ₈ − Ω_m degeneracy in the shear data. This substantially reduces the uncertainty of Ω_m in the lensing kernel while, at the same time, using only z ≲ 1 data. Furthermore, for error realisations of the source distributions, we shift the histograms in Fig. 1 according to p_z⁽ⁱ⁾(z)→p_z⁽ⁱ⁾(z + δ_zⁱ), for random draws of δ_zⁱ confined by the error model in Hildebrandt et al. (2021), their Figure 6 and Table 3; the root-mean-square (RMS) error of δ_zⁱ is σ_zⁱ ≈ 1.3 × 10⁻² for all tomographic bins, except for i = 3 where σ_zⁱ = 2 × 10⁻². Within this Monte Carlo process, we fixed the reference P_δ^fid(k, z) for all chains to that in Table 1 to avoid conflicting definitions of f_δ, mn while varying Ω_m.

After all MCMCs with randomised lensing kernels and IA parameters were available, we combined the n_runs runs by selecting from each chain j a number of n_merge = 10³ sampling points π_ij with probability w_ij/∑_iw_ij. Selected points were put back into the sample j, to be possibly selected again. By doing so, the selected n_merge points sample the posterior of chain j when equally weighted. Mixing the MCMC points from all n_run chains for the final Monte Carlo sample then contains n_merge × n_runs = 5 × 10⁵, equally weighted sampling points π_ij of the marginalised posterior distribution of f_δ, mn.

4. Data

We used the public data from the fourth Data Release (DR4) of the Kilo-Degree Survey (Kuijken et al. 2019), commonly referred to as the KiDS-1000 data¹. The images underlying the data were taken by the high-quality VST-OmegaCAM (Kuijken 2011), covering a total area of about 1006 deg² in four optical filters (ugri). After masking, the effective area of the KiDS-1000 data is 777.4 deg². Its overlap with the VIKING survey (Vista Kilo-degree Infrared Galaxy survey, Edge et al. 2013), which observes in five near-infrared bands, ZYJHK_s, allows a better control over redshift uncertainties (Hildebrandt et al. 2021).

4.1. Source catalogue

The following summarises the production of the shear catalogue, carried out by the KiDS-1000 Consortium. The images were processed using the THELI (Erben et al. 2005) and ASTRO-WISE (Begeman et al. 2013) pipelines, source ellipticities were estimated by lensfit (Miller et al. 2013; Fenech Conti et al. 2017). Based on the nine-band catalogue, individual photometric redshift estimates, z_B, were obtained using BPZ (Benítez 2000), dividing the sources into five tomographic bins (Fig. 1).

To calibrate the redshift distribution of each tomographic bin, Wright et al. (2020a) construct a self-organising map (SOM) linking the nine-band photometry of the sources to a representative spectroscopic sample. Galaxies without matching spectroscopic counterparts, or for which the z_B are catastrophically different from the redshift of the matched spectroscopic sample, are expunged from the catalogue (Wright et al. 2020b). The mean and the scatter, δ_zⁱ, on the bias of the means of the p_z⁽ⁱ⁾(z), estimated by the SOM, are numerically estimated using a suite of KiDS-like mock observations (Hildebrandt et al. 2021)². The resulting ‘gold sample’ contains ∼2.1 × 10⁷ sources with well-calibrated photometric redshift distributions; Giblin et al. (2021) gives more details. Table 2 in the Appendix summarises the main parameters relevant for our analysis. Cosmological parameter constraints from the KiDS-1000 tomographic data, with the methodology in Joachimi et al. (2021), are presented in Asgari et al. (2021), exclusively using cosmic shear, in Heymans et al. (2021), using 3 × 2pt data, and in Tröster et al. (2021), for a beyond ΛCDM analysis using 3 × 2pt data.

4.2. Data vector and covariance of measurement noise

Based on the tomographic shear data, the KiDS-1000 DR4 provides estimates of the binned shear correlation function, ξ_±^(ij)(θ), for every combination of tomographic bins, (ij), employing the popular computer software TreeCorr (Jarvis et al. 2004; Jarvis 2015). We used these data without modification. For the ξ_± estimator details, we refer to Asgari et al. (2021), Section 2.1 and Appendix C, and only note that we analysed the ξ_±^(ij)(θ) binned into N_θ = 9 angular log-bins for the range $\theta =0\stackrel{\prime }{.}5$–300′. Another notable technical detail, the data vector for the DR4 was obtained by rebinning a previously finely binned vector, instead of rerunning TreeCorr. Further, Asgari et al. (2021) exclude scales smaller than 4′ for ξ₋ due to high systematic uncertainties in the theoretical matter power-spectrum on small scales, whereas our model-free measurement of P_δ(k, z) included scales down to $0\stackrel{\prime }{.}5$.

The covariance matrix of noise, C in Eq. (25), for the estimator of ξ_±^(ij)(θ) used the analytical model by Joachimi et al. (2021), Appendix E. This model accounts for the intrinsic shape scatter of sources, for both Gaussian and non-Gaussian cosmic variance, mostly relevant towards larger angular scales, for the super-sample covariance due to fluctuation modes greater than the survey footprint, and for the calibration uncertainty in the multiplicative shear bias. As for ξ_±^(ij)(θ), the covariance matrix is part of the package of additional DR4 data products.

4.3. Systematic shear error

In addition to intrinsic alignments, a residual systematic error in the shear signal through the data reduction or the instrument may also be potentially present in ξ_±^(ij). Detecting B-modes in the shear field, specifically $\stackrel{\sim }{\kappa }(\mathit{\ell })\ne {\stackrel{\sim }{\kappa }}^{*}(-\mathit{\ell })$ as defined in Eq. (2), is considered a strong indicator of systematic errors, because B-modes generated by weak gravitational lensing are negligible. Clearly separating B-modes from E-modes in a ξ_±^(ij) with finite support is, however, unfeasible, but other related statistics are known to separate both modes (and ambiguous modes), most prominently COSEBIs (Schneider et al. 2010). For KiDS-1000, logarithmic COSEBIs capturing scales ℓ ≲ 10³ do not indicate B-modes on a 95% confidence level (Asgari et al. 2021; Giblin et al. 2021), several other tests for systematic errors are negative, and the cosmological analysis with the E-mode COSEBIs yields results consistent with a ξ_±^(ij) analysis. Therefore, a significant bias of parameters in a cosmological analysis due to a systematic shear error is unlikely.

Notably, however, the COSEBIs cosmological analysis alone is insensitive to small angular scales of ℓ ∼ 10⁴, included in our analysis through θ ≲ 4′ in ξ₋^(ij). In addition, even a small systematic shear error in ξ_±^(ij) might cause artificial deviations of P_δ(k, z) from the reference P_δ^fid(k, z), i.e., f_δ, mn ≠ 1, if the error is of the order of the difference signal between the measured ξ_±^(ij) and the expected signal in Eq. (14) for f_δ, mn ≡ 1. The 1σ–2σ window left open for B-modes by the logarithmic COSEBIs therefore leaves the possibility that deviations in KiDS-1000 might partly be related to a systematic shear error in the data. We return to this caveat in Sect. 6.

5. Results

This section reports our results for the relative power spectrum, f_δ, mn, in two variants: an entirely average, non-evolving scenario (one broad redshift bin) and an evolving scenario with three independent redshift bins. In addition, we verified our pipeline with results from a mock analysis based on ray-traced N-body data in an entirely ΛCDM framework.

5.1. Reference ΛCDM matter power-spectrum

For the KiDS-1000 analysis presented here, the reference power spectrum, P_δ^fid(k, z), is a revised halofit model with the parameters in Table 1 (Takahashi et al. 2012; Smith et al. 2003). These parameters are maximum-posterior values taken from Table C.1 (3 × 2 pt, joint) in Heymans et al. (2021) with one modification: the value for σ₈ was lowered from 0.76 to 0.72 to achieve the best fit to KiDS-1000 data that averages out to ⟨f_δ, mn⟩_mn ≈ 1 over all m and n (the next Sect. 5.2 provides more details). This σ₈ reduction is presumably required because the P_δ(k, z) model in Heymans et al. (2021) includes neutrino suppression and baryonic feedback, missing here, and is, even without neutrinos and baryons, already systematically lower by roughly 5%–10% within k ≈ 0.1–10 h Mpc⁻¹ compared to halofit. We refer to Section 4 in Mead et al. (2015) for a discussion on the latter aspect, and reiterate that our P_δ^fid(k, z) choice is arbitrary – in our analysis, a purely ΛCDM model that matches the KiDS-1000 data. The results for f_δ(k, z) in the following sections consequently probe the deviations of the true P_δ(k, z) relative to the best-fitting ΛCDM reference.

5.2. Evidence for deviation from reference

A hypothesis test against null models with constant $f_{\delta ,mn}={\overline{f}}_0$, for all m and n, and fixed projection parameters q shows that the KiDS-1000 data support deviations from the reference power spectrum in either k or z, shown by the results in Fig. 3. For each ${\overline{f}}_0$ on the x-axis in this figure, 10⁵ noise realisations of d_null were used in the null model to predict the cumulative distribution C[t_null] of the test statistic $t_{\mathrm{null}}:=-2lnP(f_{\delta ,mn}={\overline{f}}_0|{\mathit{d}}_{\mathrm{null}})$. The y-axis shows the probability, p := C[t_d], that the null model exceeds the $t_d:=-2lnP(f_{\delta ,mn}={\overline{f}}_0|\mathit{d})$ in the actual data (solid line) or, for verification, in random null-model data (dashed line). Starting with the reference power spectrum (${\overline{f}}_0=1$), as indicated by the yellow vertical line, the KiDS-1000 data have p = 0.012 and the verification data p = 0.78. Therefore, while the verification data are well within the expectation of the null hypothesis, the KiDS-1000 data are inconsistent with reference model on a 1 − p = 98.8% confidence level. Varying ${\overline{f}}_0$, however, would not improve the test result because the KiDS-1000 curve already peaks at ${\overline{f}}_0\approx 1$: the reference power spectrum, P_δ^fid(k, z), is the best fit out of a family of models with constant f_δ, mn. The only way to improve the fit would be to make f_δ, mn varying with either scale, k (this means m), or redshift, z (this means n).

Fig. 3. Test of data against a null model with identical $f_{\delta ,mn}={\overline{f}}_0$ for all m and n (Nz = 3). Shown on the y-axis is the probability of $-2lnP(f_{\delta ,mn}={\overline{f}}_0|\mathit{d})$ being greater than in the null model. The dotted line is the result for one noisy verification data vector that has ${\overline{f}}_0=1$; the solid line is for the KiDS-1000 data. The set-up in Table 1 and the source distributions in Fig. 1 without errors in the projection parameters are assumed.

In fact, the peak at ${\overline{f}}_0\approx 1$ was our deliberate choice for KiDS-1000, achieved by lowering σ₈ for P_δ^fid(k, z) to 0.72 compared to the best-fitting value of 0.78 in Heymans et al. (2021) in order to move the peak from ${\overline{f}}_0\approx 0.9$ to its final location ${\overline{f}}_0\approx 1$. This way, the f_δ, mn in our data are normalised to the average ⟨f_δ, mn⟩_mn ≈ 1. The failed null test still stands, however, and is evidence for the existence of a better fit by introducing either a k- or z-dependence of f_δ, mn, or both. This is explored in the following two more complex scenarios.

5.3. Average deviations over full redshift range

The first scenario introduces a k-dependence for f_δ, mn without redshift evolution (N_z = 1), reported with 68% CI about the median in the left panel of Fig. 4. In other words, the redshift evolution of P_δ(k, z) is assumed to be that of the reference power spectrum. Towards non-linear scales k ≳ 0.1 h Mpc⁻¹, there is no indication for a deviation from the reference exceeding the 68% CI, in particular the suppression of power does not fall below f_δ ≈ 0.8, nor do the data support a boost of power beyond f_δ ≈ 1.2. Towards linear scales, below k ≈ 0.1 h Mpc⁻¹, is a slight, yet insignificant, rise towards f_δ ≈ 1.5. This might be a sign that halofit with higher σ₈ provides a slightly better fit on linear scales. Yet, on average, variations relative to P_δ^fid(k, z) stay within about 20% (68% CI) down to k ∼ 10 h Mpc⁻¹.

Fig. 4. Reconstructed matter power spectrum in KiDS-1000 for τ = 5.0 and Nk = 20 in three variants. The errors marginalise over uncertainties in the lensing kernel and IA. Shown are the posterior 68% CI about the median (lines). Left panel: Transfer function fδ(k, z) averaged over the entire redshift range Z1 = [0, 2] and relative to the reference, Pδfid(k, z), with parameters listed in Table 1. Middle panel: Same as the left panel, but for three independent transfer functions averaged separately within the ranges Z1 = [0, 0.3], Z2 = [0.3, 0.6], and Z3 = [0.6, 2]. Right panel: Dimensionless power spectrum, Δ2(k, z) = 4π k3Pδ(k, z) = 4π k3 fδ(k, z) Pδfid(k, z), interpolated to the centres of three redshift bins Z1, Z2, and Z3.

5.4. Deviations for three separate redshift bins

For more insight in a possible z-evolving f_δ in the second scenario, we probed the average f_δ(k, z) inside three separate redshift bins, Z₁ = [0, 0.3], Z₂ = [0.3, 0.6], and Z₃ = [0.6, 2]. Their marginalised posterior constraints (68% CI about the median) are plotted as three credible regions in the middle panel of Fig. 4, one region for each redshift bin. The lowest bin, Z₁, is fully consistent with f_δ = 1, on k-average f_δ = 1.15 ± 0.28, while the middle bin, Z₂, on average f_δ = 0.57 ± 0.27, indicates a power suppression, reaching the depth f_δ = 0.3 ± 0.2 at k ∼ 0.5 h Mpc⁻¹, and the highest bin, Z₃, overall prefers a boosted power of f_δ = 2.22 ± 0.81. The signal suppression in the middle redshift bin and the boost in the highest bin cancel each other when averaging f_δ over the whole redshift range, resulting in the f_δ = 0.99 ± 0.20 in the left panel.

The close description of the ξ_±^(ij)(θ) data points by the flexible f_δ model is illustrated by the posterior predictive of the model in Fig. B.1 (Appendix). Shown in the various panels are our KiDS-1000 data points (black with 1σ errors bars) for θ ξ₋^(ij)(θ) (lower left panels) and θ ξ₊^(ij)(θ) (upper right panels) in comparison to the posterior model constraints (68% and 95% CIs in dark and light blue). The CIs in this plot do not marginalise over q; however, this has an effect of less than 10%. The red line corresponds to the ΛCDM reference model, this means f_δ, mn ≡ 1, obtained from a best fit of σ₈ to the data points (Sect. 5.1). Allowing for deviations f_δ, mn ≠ 1 in the three separate redshift bins Z₁ to Z₃ moves the CIs relative to the reference, thereby improving the fit, although the reference still stays within the 95% CI, as, for instance, for ξ_±⁽⁵⁵⁾(θ). Whilst an overall good match to the (correlated) data points, a conflict with the model is probably present for ξ₊⁽²²⁾(θ), and perhaps in ξ₊⁽¹²⁾(θ) to ξ₊⁽¹⁵⁾(θ), where the data prefer a higher amplitude than the model fitting all tomographic bins. Figure A.4 shows a similar plot with verification data, based on the reference model, as comparison. A higher value of $S_8:=\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}{\sigma }_8\approx 0.791$, as indicated by the green line for the N-body verification data (without IA, Sect. 5.5), would account for the higher ξ₊^(1j)-amplitude but, on the other hand, is rejected by the data due to mismatches at higher redshift (e.g., ξ₊⁽³⁵⁾ or ξ₊⁽⁴⁵⁾). This might indicate spurious systematic errors in the low-z bins, discussed in Sect. 6.

Expressed in terms of the actual power spectrum, ${\mathrm{\Delta }}^2(k,\overline{z}):=4\pi k^3{f}_{\delta }(k,\overline{z})P_{\delta }^{\mathrm{fid}}(k,\overline{z})$, interpolated to the z-bin centres, $\overline{z}$, the diverse redshift evolution of f_δ(k, z) in the KiDS-1000 data translates into a suppression of non-linear structure growth between 0.3 ≲ z ≲ 2. Namely, the $\overline{z}=0.15$ matter power spectrum in the right panel of Fig. 4 (magenta) has clearly more power than the two higher redshift bins (green for $\overline{z}=0.45$ and cyan for $\overline{z}=1.3$) for k ≳ 0.1 h Mpc⁻¹, while the two other power spectra are consistent with each other, maybe with the green rising over cyan near k = 10 h Mpc⁻¹. Consequently, with cyan and green being moved closer to each other, there is for k ≳ 0.1 h Mpc⁻¹ less structure growth detected between $\overline{z}=0.45$–1.3 than in the reference ΛCDM cosmology but, conversely, more between $\overline{z}=0.15$–0.45, pulling back the low-z bin to the reference power spectrum.

It is noteworthy that the statistical errors of f_δ, mn and hence Δ²(k, z) are correlated, best visualised by the correlation matrix in Fig. 5. Each pixel in the correlation matrix encodes the Pearson correlation coefficient, r, between two k-bins, either from same or from different z-bins. The matrix is organised in N_k × N_k blocks, where k increases from left to right and top to bottom. The three blocks on the diagonal, bottom left to top right (lowest z to highest z), are correlations within the same redshift bin, and the six off-diagonal blocks are cross-correlations between distinct z-bins. The enhancement of correlations within the blocks, making them stick out as tiles in a 3 × 3 mosaic, is a side-effect of the Tikhonov regularisation smoothing the f_δ, mn along the k-axis. Correlations are high for adjacent k-bins in the same z-bin, r ≈ 0.4–0.8, especially in the highest z-bin (top right), but quickly drop off with distance on the k-scale. Adjacent k-bins are also statistically dependent when from adjacent z-bins, albeit now anti-correlated with typically r ≈ −0.4 to −0.1. Not shown is the correlation matrix for a N_z = 1 reconstruction, left panel in Fig. 4, to save space; it is essentially a diagonal matrix where adjacent k-bins have |r|≲0.1, except at the lower and upper boundary of the k-range where r ≈ 0.5–0.6.

Fig. 5. Correlation matrix of fδ, mn errors for the KiDS-1000 analysis with three redshift bins, corresponding to the middle and right panels in Fig. 4. The three blocks on the diagonal from bottom left to top right are the correlation matrices for the z-bins Z1 = [0, 0.3], Z2 = [0.3, 0.6], Z3 = [0.6, 2], respectively; off-diagonal blocks show correlations between distinct z-bins; pixels within the 20 × 20 blocks are correlations between the k-bins, ascending from left to right inside a block.

5.5. Verification test with ray-traced mock data

To add more realism to the code verification than in Appendix A, especially with respect to the f_δ result for N_z = 3, we performed a more stringent test with N-body ray-traced data (Takahashi et al. 2017). As in Burger et al. (2024), we subdivided the 108 full-sky shear fields from the simulated data according to the KiDS-1000 sky area and, for each tomographic bin, combined 30 shear grids for ascending source redshifts according to the redshift distributions in Fig. 1. This yields an ensemble of n = 1944 mock KiDS-1000 surveys. From these, source positions were picked to match the angular positions, shear weights, and shape noise to the measured counterparts in the public KiDS-1000 ellipticity catalogue (Giblin et al. 2021). In contrast to Appendix A, however, intrinsic alignments of sources are not included here, thus A_IA = 0. Similar to the actual KiDS-1000 data vector, we used TreeCorr to compute ξ_±^(ij)(θ) from the mock shear catalogues.

Driven first by the question whether our code can recover f_δ(k, z) = 1 from a simulated (unrealistically) high S/N data vector when using exactly the same cosmological parameters as in the simulation, we averaged all n realisations of the data vector, denoted by ${\overline{\mathit{d}}}^{\mathrm{sim}}$, and used the individual realisations, d_i^sim, to estimate the (standard) error covariance of ${\overline{\mathit{d}}}^{\mathrm{sim}}$,

$$\begin{array}{c}\hfill {\mathsf{C}}^{\mathrm{sim}}=\frac{1}{n(n-1)}{\displaystyle \sum _{i=1}^n}({\mathit{d}}_i^{\mathrm{sim}}-{\overline{\mathit{d}}}^{\mathrm{sim}}){({\mathit{d}}_i^{\mathrm{sim}}-{\overline{\mathit{d}}}^{\mathrm{sim}})}^{\mathrm{T}}.\end{array}$$(28)

This covariance matrix was turned into an unbiased estimator of the inverse covariance using Hartlap et al. (2007). The average ${\overline{\mathit{d}}}^{\mathrm{sim}}$ and C^sim are basically a single measurement in a hypothetical survey with n times the angular area than KiDS-1000, albeit underestimating the cosmic variance error. For an initial inspection, the solid green lines SIM-THS17 in Fig. B.1 (Appendix) compares ${\overline{\mathit{d}}}^{\mathrm{sim}}$ to the real KiDS-1000 data (black data points with error bars): The simulated tomography correlations are mostly within the 1σ noise scatter of the KiDS-1000 data points, yet exhibit systematically more signal compared to the reference model (red lines) by occasionally exceeding the 95% CIs, probably due to the higher fluctuation amplitude, S₈ ≈ 0.791 in the simulation compared to S₈ ≈ 0.765 in KiDS-1000, and, less relevantly, the missing IA terms in the simulation. Now, using ${\overline{\mathit{d}}}^{\mathrm{sim}}$ and C^sim as input, the data points in the top panel of Fig. 6 show the code results for f_δ(k, z) in three redshift bins. To test the code’s ability to recover the power spectrum in the simulation, f_δ(k, z) = 1, the reference P_δ^fid(k, z) was set here to the same cosmological parameters as in the N-body data, namely Ω_m = 0.279, Ω_Λ = 1 − Ω_m, Ω_b = 0.046, h = 0.7, n_s = 0.97, and σ₈ = 0.82. Since the lensing kernel is exactly known in this experiment, no marginalisation of q was done. We find the following.

Fig. 6. Verification test of the analysis code for fδ(k, z) using: an N-body simulated data vector, averaged over n = 1944 realisations of the KiDS-1000 data without IA; the code and data set-up for the KiDS-1000 analysis has Nk = 20 k-bins and three redshift bins Z1 = [0, 0.3], Z2 = [0.3, 0.6], and Z3 = [0.6, 2]. The reference power, Pδfid(k, z), is here identical to the cosmology of the simulation (Takahashi et al. 2017). The errors (68% CI) do not marginalise over a projection kernel uncertainty. Top panel: Simulated covariance of the average data, comprising n times the area of KiDS-1000. Deviations from the expected fδ(k, z) = 1 towards higher z and k ≳ 1 h Mpc−1 are related to the finite grid pixel size in the simulation (pixelation bias). The lines are additional reconstructions using noise-free data, based on Eq. (8) and halofit without pixelation bias. Bottom panel: Same as top panel, but with statistical errors as in a single KiDS-1000 survey.

The reconstructed power spectrum in Fig. 6, top panel, for the lowest redshift bin (filled squares) is overall consistent with f_δ, mn = 1 for the range k = 0.01 − 20 h Mpc⁻¹ within the 68% CI, indicated by the error bars. But two significant deviations emerge: first, in the higher redshift bins beyond k ≈ 5 h Mpc⁻¹, and, second, smaller deviations around 0.01 h Mpc⁻¹. The high-k suppression of the signal is a real feature in the simulation data due to the pixelated shear grids of the simulation (pixel size $0\stackrel{\prime }{.}4$), significantly biasing low the values for ξ₋^(ij) for θ ≲ 5′. The lowest redshift bin is presumably also affected, but scales of a few arcmins correspond here to a characteristic k too high to be visible in the plot: the typical wave number, k_eff = 2π/(f_K(χ_d) θ_1′), of an angular scale θ_1′ = 1′ on a lens plane at distance χ_d corresponds to $k_{\mathrm{eff}}\approx 33.5,12.0,5.2h{\mathrm{Mpc}}^{-1}$ for lens planes located at the centres of the bins Z₁, Z₂, and Z₃, respectively. The additional lines in the plot support the idea of a pixelation bias: they are reconstructions of noise-free data, based on the simple Limber-Kaiser projection in Eq. (8) – crucially, not showing the systematic signal drop for Z₂, or for Z₃ around k = 1 h Mpc⁻¹. Instead, however, the lines exhibit artefacts as zig-zag oscillations of amplitude Δf_δ ≈ 0.1, strongly increasing towards the edges of the plotted k-range and probably related to the limitations of the deprojection method (Sect. 3.1). The artefacts are partly mirrored in the data point scatter and overlap with the pixelation bias for Z₃ (yet not for Z₂) at large k.

The second, weaker ∼1.5σ–2σ deviations appear around k = 0.01 h Mpc⁻¹ for Z₂ and Z₃. They, too, are probably to some extend present in the N-body data: Figure 19 in Takahashi et al. (2017) report a 10% deviation from halofit near k = 0.01 h Mpc⁻¹. Nevertheless, we doubt that this explains the full effect and point out that the plots compare an average data vector in a simulation of finite volume, still subject to noise, to the cosmic average. Especially on large scales, cosmic variance noise is relevant and under-estimated in C^sim. In addition, the noise-free data (solid lines) suggest some excess signal on large scales due to edge artefacts. In summary, the bottom panel of Fig. 6 is an excellent reconstruction of the relative power with Δf_δ ≈ 0.1 accuracy, or better, for k ≲ 10 h Mpc⁻¹ and z ≲ 1.

Another inference of f_δ(k, z) with the simulated data is shown in the bottom panel of Fig. 6 – this time for a KiDS-1000-like error covariance, n × C^sim. By using the average ${\overline{\mathit{d}}}^{\mathrm{sim}}$ as input, the (essentially linear) reconstruction in the figure represents the ensemble mean of f_δ(k, z) posteriors in three redshift bins. In this average reconstruction, the posterior errors are similar to the ones in the KiDS-1000 reconstruction, middle panel of Fig. 4, increasing in size for higher z, although now the f_δ, mn are consistent with f_δ = 1 throughout (68% CI). The suppression of power by pixelation bias is here fully contained within the statistical errors, only vaguely hinted at in the green Z₂ posterior region of Fig. 6, bottom panel, which falls just below f_δ = 1 around k = 10 h Mpc⁻¹. Therefore, also for KiDS-1000-like noise levels, the code recovers the correct power spectrum in the N-body data. The highest precision of f_δ = 1.00 ± 0.27 is achieved for Z₁, the middle bin Z₂ has f_δ = 0.81 ± 0.31, and Z₃f_δ = 1.29 ± 0.63 (medians and 68% CIs, no q errors).

For comparison, Fig. A.5 (Appendix) shows an ensemble of reconstructions with individual, 16 randomly chosen d_i^sim, illustrating the possible variations in the inferred f_δ(k, z). Each panel is a prediction of a KiDS-1000 reconstruction in a ΛCDM scenario. Here, we observe a small tendency of the posterior median in Z₃ to be above f_δ = 1 and the median in Z₂ to be below f_δ = 1, despite both having f_δ ≈ 1 for k ≲ 5 h Mpc⁻¹. This tendency is also visible in Fig. 6 of average reconstructions, bottom panel, with the posterior median always being above f_δ = 1 for Z₃, or below for Z₂, and probably has two reasons. First, the Z₃ bin is poorly constrained, giving a broad marginalised posterior PDF skewed towards higher f_δ values, shifting up the centres (posterior medians) of our CIs; this is illustrated by Fig. A.3 in the Appendix. Second, errors in Z₂ and Z₃ are anti-correlated, preferring a lower value for Z₂, if the f_δ at similar k in Z₃ is scattered upwards. But, on a broad note, one should not expect the maximum of the marginalised PDF to coincide with the maximum of the full posterior (for the non-marginalised parameters) because this requires that the full posterior PDF must obey certain symmetries about its maximum, as, for instance, present in a multivariate Gaussian posterior – excluded here because of our skewed, non-Gaussian marginalised posterior distributions. For increased S/N in the data, however, the long tails of skewed posteriors are suppressed, and the median shift is no longer discernible in the top panel of Fig. 6.

5.6. Impact of uncertainty in IA and lensing kernel

The individual impact by uncertainties in the projection parameters q on the reconstructed f_δ, mn is of interest for improving future analyses or to gauge the relevance of errors in the projection kernel. We therefore compare the marginalised f_δ, mn errors in the middle panel of Fig. 4 to the statistical errors obtained without marginalisation, or to the errors where all parameters q are fixed but for one that is marginalised over.

Figure 7 shows the results in these different scenarios as function of scale k for three z-bins, indicated in the lower right of each panel; a 100% error is as large as in a scenario where all q uncertainties are marginalised over. The solid black lines use a fixed q without any marginalisation, setting a lower limit for all scenarios (but is also subject to numerical noise owing to the limited number of realisations, n_merge = 500). A lower limit of 90% to the black line for all k means that uncertainties due to source shape noise, cosmic variance, and calibrated shear bias, all accounted for in the noise covariance C, are responsible for at least 90% of the total error in the reconstruction, while the addition of lensing kernel uncertainties and IA errors amounts to 10% or less of the total. This addition is explored in more detail by the dashed lines that marginalise over redshift errors only, the dotted-dashed lines for only Ω_m uncertainties in the lensing kernel, and the dotted lines for A_IA errors only. In summary, IA errors have a subdominant impact, about several per cent, reflected by the closeness of the dotted line to the solid line. The largest impact of the Ω_m uncertainty is in the lowest z-bin, whereas the uncertainty in the source z-distributions is the dominating contributor in the middle and highest z-bin for k ≲ 1 h Mpc⁻¹. Therefore, better calibrating photo-z errors in future analyses yields improvements of up to 10% in the total statistical error.

Fig. 7. Percentage fraction in total marginal posterior error of fδ, mn due to uncertainties in the lensing kernel and IA parameter. Shown is the statistical error relative to the full marginal error (RMS variance of posterior) as a function of k in three redshift bins without marginalisation (solid lines), with marginalisation over errors in the source redshift distributions only (dashed lines), Ωm in the lensing kernel only (dash-dotted lines), and for the AIA error only (dotted lines).

6. Discussion

In this study we presented, verified, and applied to KiDS-1000 data a revised method for reconstructing the three-dimensional (band) matter power spectrum with tomographic shear data. Our approach relies on the specification of a projection kernel only, and is hence agnostic about an analytical model for P_δ(k, z), while being applicable for all cosmologies where the ξ_±^(ij)(θ) are projections of P_δ(k, z) along the line of sight. Numerically convenient, we express P_δ(k, z) =  : f_δ(k, z) P_δ^fid(k, z) in terms of a reference and probe deviations, f_δ, mn ∈ [0, 100], relative to it, putting the free variables f_δ, mn into the same dynamical range. Not required, however, is a reference close to the true power spectrum; for instance, in Fig. 2 we use a reference fixed in cosmic time to highlight the structure growth, or, on the extreme end, a constant P_δ^fid(k, z)≡P₀ would yield a f_δ, mn that directly represents the power spectrum. The rescaling of P_δ^fid(k, z) by a constant inside any z-bin, whilst redefining f_δ, mn, does not affect the inferred Δ²(k, z)∝f_δ(k, z) P_δ^fid(k, z); the reconstructed matter power spectrum is robust against the reference used. Choosing a reference towards a slowly evolving f_δ(k, z) is nevertheless preferable because the z-bin-average that the f_δ, mn represent depends on the specific noise properties and source z-distributions in the tomographic survey; a comparison to theoretical models will have to replicate the data-specific z-weighting more closely for a significantly varying f_δ(k, z) inside a wide z-bin. This is obvious in Fig. 2 for a z-independent reference, where more weight is given closer to the lower boundary of each z-bin: the average f_δ, mn is close to the true f_δ(k, z) near z_max = 0.13, 0.4, 0.7 for Z₁, Z₂, and Z₃; yet the bin centres are at z_mid = 0.15, 0.45, 1.30. Although a best-fitting ΛCDM model without baryonic feedback or neutrinos as reference is a good fit to the data already, as shown by the red lines relative to the black data points in Fig. B.1, the f_δ, mn values indicate where additional changes might be needed or how much deviation from the reference is acceptable until we conflict with the data. We proceed here with two variants of the analysis, one for an f_δ, mn-average over the full redshift range and one to probe the z-evolution of f_δ(k, z).

Our first variant in the left panel of Fig. 4 averages f_δ(k, z), denoted ${\overline{f}}_{\delta }(k)$ hereafter, over the entire z-range covered by KiDS-1000, with the most signal from essentially z ≲ 1. This analysis shows that in the non-linear regime, k ≈ 0.05–10 h Mpc⁻¹, the KiDS-1000 matter power-spectrum is consistent with the ΛCDM reference in Table 1 within 20% (68% CI). Therefore, there is no evidence for a variation, suppression or boost, exceeding 20% relative to the reference scenario with cold dark matter (excluding neutrino suppression and baryon feedback). This is also within the range reported by other lensing studies, discussed below, although we caution that the results are not always one-to-one comparable to ours because they are partly based on the premise that the lensing P_δ(k, z) has the PlanckS₈ (Planck Collaboration VI 2020), but is attenuated in the non-linear regime to fit the lensing measurements of significantly lower S₈ to explain the so-called S₈ tension. Such an attenuation, if present, would be addressed in our reference by a systematically lower S₈ instead of fully showing up in ${\overline{f}}_{\delta }(k)$, especially because our data lack constraining power in the unattenuated linear regime, k ≲ 0.01 h Mpc⁻¹. To level the field with analyses using the PlanckS₈ as the premise, we changed the reference to a power spectrum with Planck cosmology, Δ_Planck²(k, z), and recast Fig. 4, left panel, into Fig. 8 for the ratio ${\mathrm{\Delta }}^2(k,\overline{z})/{\mathrm{\Delta }}_{\mathrm{Planck}}^2(k,\overline{z})$ at $\overline{z}=0.5$, a redshift where the lensing method is roughly the most sensitive. Only from this Planck perspective do our results now clearly require a suppression, 1 − f_δ, mn, of the ΛCDM reference at k ≈ 3 h Mpc⁻¹ of up to 35%±15%±25% (68% and 95% CI). As indicated by the lines in the plot, it is conceivable to explain the tension with Planck by a variety of hypothetical mechanisms involving, among other things, extreme baryonic feedback, (sterile) neutrinos, or alternatives to non-interacting cold dark matter.

Fig. 8. Power spectrum ratio ΔKiDS2(k, z)/ΔPlanck2(k, z) at z = 0.5. Shown in grey are the 68% and 95% CI from our analysis. The KiDS-1000 spectrum ΔKiDS2(k, z) uses the fδ, mn in Fig. 4 averaged within one broad redshift bin (left panel). The Planck model for ΔPlanck2(k, z) is the TT, TE, EE+lowE best fit of a standard ΛCDM (Planck Collaboration VI 2020) extrapolated to z = 0.5 with HMCode2020 by Mead et al. (2021). The lines are model transfer functions to ΔPlanck2(k, z) invoking flavours of extreme baryonic feedback (HMCode2020: log10(TAGN/K) = 8.8, Abaryon = 3.13, ηbaryon = 0.603; BCM by Schneider & Teyssier 2015: log10(Mc/h−1 M⊙) = 15.0792, ηb = 0.5, ks = 55.0), neutrino suppression for cumulative neutrino masses ∑νmν (BCM), or a mixture of CDM and WDM (Kamada et al. 2016, mwdm: WDM particle mass, fwdm: WDM mass fraction).

Therefore, whereas a CDM power spectrum without suppression and low S₈ is on average a good description of the KiDS-1000 shear tomography, small modifications to CDM would be insufficient to reconcile KiDS-1000 with Planck: mechanisms with a substantial suppression (more than 20%) have to be invoked, in agreement with other recent studies analysing KiDS-1000 lensing data. For instance, the recent Schneider et al. (2022, S+22) combines KiDS-1000 data with observations of X-ray emission and the kinematic Sunyaev–Zeldovich (kSZ) effect in galaxy clusters to inform baryonic feedback parameters incorporated in their ΛCDM-flavoured P_δ(k, z). By initially using the KiDS-1000 data alone, the deviations from their model without baryonic feedback are inconclusive but remain, as in our measurement in the left panel of Fig. 4, within 10–20% for k = 1–10 h Mpc⁻¹ (68% CI). Only when adding the X-ray and kSZ data do their constraints strongly favour a power suppression of up to 30% (68% CI) at k ≈ 5 h Mpc⁻¹, somewhat more than our data support. Building only on lensing information, Preston et al. (2023, PEG23) analyse KiDS-1000 and data from the third-year Dark Energy Survey (DES Y3) to extend a ΛCDM model by one modification parameter, A_mod, pushing the non-linear P_δ(k, z), equivalent to A_mod = 1, towards the linear power spectrum, A_mod = 0. Their KiDS-1000 fit with uninformative priors on cosmological parameters and A_mod ∈ [ − 5, 5] is, as in Fig. 4, inconclusive, A_mod ≈ 0.9 ± 0.2 (their Figure 2, left panel). However, when assuming a Planck cosmology and the suppression-inclined prior A_mod ∈ [0.5, 1], PEG23 claim strong evidence for a suppression of ≈25%±6% at k ≈ 3 h Mpc⁻¹ (68% CI), more significant yet fully consistent with our Fig. 8 (we reiterate our relaxed prior f_δ, mn ∈ [0, 10²]). Furthermore, Broxterman & Kuijken (2024, BK24) fit a double power law P_δ(k, z)∝k^p (1 + z)^−m_BK to ξ_±^(ij)(θ) for KiDS-1000, DES Y3, and the third year Hyper Suprime-Cam Survey (HSC Y3), yielding a reasonable ΛCDM match without suppression in the range k = 0.01–5 h Mpc⁻¹ (68% CI) and a 2σ agreement among the three Stage III surveys. Compared to our results for N_z = 1, similarly averaging over the full redshift range, their constraints are weaker in the non-linear regime ($f_{\delta }\approx 0.8_{-0.2}^{+0.5}$), despite the same data, probably because of our tighter Ω_m prior for the projection kernel: the P_δ(k, z) amplitude in ξ_±^(ij)(θ) is degenerate with the pre-factor Ω_m² in the Eqs. (8)–(10). Finally, similar results are reported by authors that use only DES Y3 lensing data although the required suppression in the non-linear regime relative to Planck tends to be somewhat lower than for KiDS-1000 by a factor between one and two (Bigwood et al. 2024, PEG23, Ferreira et al. 2024, Aricò et al. 2023). Utilising lensing probes alone, Perez Sarmiento et al. (2025) demonstrates consistency with the Planck cosmology on linear scales, yet reveals a 20–25% suppression in the non-linear regime, by combining DES Y3 and CMB lensing data. In conclusion, other works support our finding of an average matter power-spectrum that at k ≳ 0.05 h Mpc⁻¹ is either, within 20% tolerance, a purely ΛCDM spectrum of low S₈ ≈ 0.73, or a spectrum with both higher PlanckS₈ ≈ 0.83 and significant suppression of up to 20%–30% in the non-linear regime relative to CDM.

Surprisingly, our KiDS-1000 analysis becomes more complex and discrepant with our best-fitting ΛCDM reference of low S₈ if the f_δ, mn, unlike in previous studies, are free to vary with redshift, achieved in our second analysis by increasing the number of z-bins to N_z = 3 broad z-bins. This is only practical by employing Tikhonov regularisation, or something else to this effect, as clearly illustrated by the dramatically shrinking sizes of the credible regions from τ = 0 to τ = 5 in Fig. 2. Applied to KiDS-1000, the 68% CI constraints in the middle panel of Fig. 4 then depict a remarkably diverse picture: only the f_δ, mn in the low-z bin, Z₁ = [0, 0.3], are now consistent with the reference matter power spectrum, averaged for all k-bins to ${\overline{f}}_{\delta }=1.15\pm 0.28$, whereas the power spectrum in Z₂ = [0.3, 0.6] is overall suppressed, ${\overline{f}}_{\delta }=0.57\pm 0.27$, and boosted for Z₃ = [0.6, 2], where ${\overline{f}}_{\delta }=2.22\pm 0.81$. The verification tests in Appendix A, for analytic data (noisy) vectors, and in Sect. 5.5, for mock N-body shear catalogues, assure us of an f_δ, mn accuracy of 10% or better within k ≈ 0.01–10 h Mpc⁻¹. The observed tendency of shifting the median posterior f_δ, mn of our 68% CIs towards ${\overline{f}}_{\delta }\approx 1.3$ for Z₃ and towards ${\overline{f}}_{\delta }\approx 0.8$–0.9 for Z₂, due to the broadened skewed posterior PDFs in these z-bins, actually goes in the direction of the KiDS-1000 result. It is nevertheless too small for the pronounced split between Z₂ and Z₃ credible regions in the real data, which also are, unlike the verification results, blurred by the IA and lensing kernel uncertainties (an extra ∼10% error, Fig. 7). For the baseline f_δ(k, z)≡1 in Fig. A.5, the split in the mock data realisations is not nearly as pronounced as for KiDS-1000, making the KiDS-1000 result very roughly a 1/16 ≈ 6% event, or less, for the baseline scenario.

To explore the possibility of spurious systematic errors possibly mimicking the evolving f_δ(k, z), we estimated in three scenarios the amplitude of bias needed to replicate our KiDS-1000 results if the ΛCDM reference were the true matter power spectrum. First, a spurious shear m-bias in the data is assumed to be the dominating source of systematic error. The required bias is then obtained through fitting $(1+m_i)(1+m_j){\xi }_{\pm }^{(ij),\mathrm{fid}}(\theta )$, predicted by the reference model, to the KiDS-1000 data vector by maximising the data likelihood, Eq. (25), with respect to the bias parameters m_i in the tomographic bins i. The corresponding (posterior median) f_δ, mn of the best-fitting model are plotted as data points with filled squares in Fig. 9 in comparison to KiDS-1000 (shaded regions). Evidently, the m-bias model closely traces the KiDS-1000 results, but, according to its best-fit parameters {m_i}={+1.00, +0.48, −0.023, −0.028, +0.035}, it needs an unrealistic m-bias in the first two source bins that exceeds the reported values in Asgari et al. (2021) by more than 30σ. The second scenario assumes spurious photo-z errors in p_z⁽ⁱ⁾(z) as dominating source of bias. Again, we fit a ${\xi }_{\pm }^{(ij),\mathrm{fid}}(\theta )$ to the KiDS-1000 data vector, now for the p_z⁽ⁱ⁾(z + δ_zⁱ) shifted by the bias parameters δ_zⁱ, while the reconstruction uses the original p_z⁽ⁱ⁾(z). The best-fit values are {δ_zⁱ}={−0.075, −0.079, +0.009, +0.018, −0.052}. As in the previous m-bias scenario, this model fits the KiDS-1000 results equally well (filled circles in Fig. 9), but requires a high bias (photo-zs are systematically low) in the first two and the last source z-bins. This bias level exceeds the estimated tolerances of photo-z errors in the KiDS-1000 gold sample by more than 5σ (Table 2). A third conceivable source of systematic error is an IA model bias. Too simplistic physically, the approximations in the NLA model are nevertheless flexible enough to fit measurements of IA as a means for calibrating the IA distortion of ξ_±^(ij)(θ) (Fortuna et al. 2021; Preston et al. 2024). It may be, however, that z-binning of Δ²(k, z) makes us sensitive to a z-dependent IA amplitude, missing in the NLA model and typically irrelevant in Stage III cosmological analyses. To test if a z-dependent IA amplitude A_IA(z):=A_IA′ (1 + z)^η (1 + z_piv)^−η, as in the eNLA model, potentially explains the f_δ, mn data after all, we fit ${\xi }_{\pm }^{(ij),\mathrm{fid}}(\theta )$ to KiDS-1000 by varying {A_IA′,z_piv, η} only, yielding the parameters {0.3, 0.19, 1.5} for the open diamond data points in Fig. 9. Crucially, the reconstruction algorithm is ignorant about A_IA(z) and assumes a constant A_IA. Indeed, then some of the split between Z₂ and Z₃ (especially the drop in Z₂) could be explained by the IA model bias. And, since z_piv = 0.3 and η = 1.5 are deemed consistent with KiDS-1000 uncertainties in Asgari et al. (2021), this hypothesis is a plausible explanation for at least part of the ${\overline{f}}_{\delta }$ redshift variation. More plausible still is a combination of A_IA(z) and a moderate photo-z bias, 2σ–4σ, indicated by the crosses in Fig. 9. Lastly, recalling the 1σ–2σ evidence for B-modes in Sect. 4.3, a systematic shear error might also be relevant here. However, we lack a KiDS-1000 model of these for ξ_±^(ij), and it is unclear if they plausibly contribute to the ${\overline{f}}_{\delta }\ne 1$ result. In summary, whereas (especially) m-biases and photo-z biases alone are unlikely sources of systematic errors that shape the result of evolving ${\overline{f}}_{\delta }$, the negligence of a z-dependent IA amplitude, A_IA, as in our reconstruction set-up, in combination with moderate photo-z biases might explain our N_z = 3 result. At this point, we are unable to exclude either possibility, or that of small systematic errors related to B-modes, but note that the discussed scenarios mainly target the first two tomographic bins, z_B ≤ 0.5, where the posterior predictive of the model also indicates some conflict with the data (Fig. B.1). Future reconstructions of the matter power spectrum inside separate redshift bins should also pay close attention to the IA modelling.

Fig. 9. KiDS-1000 results for fδ, mn, as in the middle panel of Fig. 4, in relation to alternative hypotheses using the reference cosmology Pδfid(k, z) and best-fitting bias parameters to match the KiDS-1000 data vector, ξ±(ij)(θ). The filled squares assume shear m-bias factors {mi}={+1.00, +0.48, −0.023, −0.028, +0.035}, the filled circles assume photo-z biases {δzi}={−0.075, −0.079, +0.009, +0.018, −0.052}, and the open diamonds employ the eNLA model with {AIA′,zpiv, η}={0.3, 0.19, 1.5}. The crosses assume a hybrid model of eNLA and a smaller z-bias: {AIA′,zpiv, η}={0.45, 0.25, 2.3} and {δzi}={−0.017, −0.058, +0.027, +0.025, −0.043}. The data points indicate the posterior median (CIs are similar to those of the shaded regions); the filled data points use alternating k-bins to avoid overlap.

A truly evolving f_δ(k, z), on the other hand, would offer insights into the physics of non-linear structure growth at k ≳ 0.1 h Mpc⁻¹ and z ≲ 1. In brief, the KiDS-1000 data indicate too little power for the epoch z ≈ 0.3–0.6 and too much power for z ≈ 0.6–2, compared to the purely-ΛCDM reference model that actually matches the f_δ, mn averaged in a single z-bin (68% CI). Phrasing this a bit more carefully, this statement applies to weighted averages inside three wide z-bins where most weight is given to Δ²(k, z) around z ≈ 0.13, 0.4, 0.7 for Z₁, Z₂, and Z₃, similar to Fig. 2. Even so, the average deviations in the Z₂ and Z₃ are, in fact, such that they make Δ²(k, z) inside the bins statistically consistent (right panel of Fig. 4), independent of the chosen reference – thus there is no significant detection of any growth from about z = 0.7 to 0.4. The only increase in power is detected towards lower z ≈ 0.13. This qualitatively matches the S₈-tension: cosmological probes that give more weight to z ≳ 0.7, such as CMB experiments, would observe a Δ²(k, z) with amplitude higher than in our reference, hence a higher S₈, whereas probes more sensitive to z ∼ 0.4 or lower, such as lensing experiments, would observe a Δ²(k, z) with S₈ equal to our reference or even lower. We are not aware of a similar conclusion in the lensing literature, although the recent BK24 report a vanishing growth rate at 1σ, their growth index m_BK, for KiDS-1000 and HSC Y3 (yet not for DES Y3). That finding for m_BK might be related to our result, depending on how the Z₁–Z₃ are effectively weighted in the BK24 analysis. Possible physical causes of the peculiar structure growth in KiDS-1000 compared to the reference model are, however, presumably non-trivial and beyond the scope of this paper. For example, a baryonic feedback model would have to be strongly evolving with redshift below z ∼ 1 (and non-monotonic) to match the KiDS-1000 f_δ, mn, typically not found in state-of-the art models (e.g., Schaller et al. 2025, and references therein).

Finally, the close connection of the S₈-tension to the KiDS-1000 result for Δ²(k, z) in N_z = 3 z-bins makes this an interesting case, either as evidence in favour of modifications towards the structure growth in the ΛCDM reference, or in favour of spurious systematic errors in a shear tomography analysis mimicking the tension. This evidence appears to be unnoticed in previous Stage III analyses – a full analytic model of Δ²(k, z) fitted to the shear tomography can give an unsuspecting match similar to our N_z = 1 analysis – so that a confirmation of the KiDS-1000 result with other lensing data is important. Furthermore, looking forward not too far into the future, Stage IV data will be available soon with roughly 20 times the survey area of KiDS-1000 and also deeper observations, with surveys such as the space-based missions Euclid and Roman, or the ground-based survey by Rubin. Of these the Euclid Wide Survey will provide tomographic measurements within 13 source z-bins up to z ∼ 2.5, statistical errors reduced by about one order of magnitude compared to Stage III, and an excellent control of systematic errors for shear and photo-z. We see no restriction to applying our reconstruction technique to Stage IV data, merely assuming a lensing kernel and an IA model (possibly updated and carefully calibrated). This promises measurements of the matter power spectrum at different redshifts to unprecedented precision, perhaps up to z ∼ 2 with N_z = 4 or more bins, critically testing models of the non-linear structure growth.

Acknowledgments

This paper went through the KiDS review process. We also want to thank Jeger Broxtermann for useful comments. LP acknowledges support from the DLR grant 50QE2002. PB acknowledges financial support from the Canadian Space Agency (Grant No. 23EXPROSS1) and the Waterloo Centre for Astrophysics. Some results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A- 3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples).

References