© The Authors 2026

1. Introduction

The Λ cold dark matter (ΛCDM) model has been remarkably successful in explaining a variety of cosmological observations, such as weak gravitation lensing (Wright et al. 2025b; Stölzner et al. 2025; Amon et al. 2022; Secco et al. 2022; Li et al. 2023; Dalal et al. 2023), cosmic microwave background measurements (CMB; Planck Collaboration VI 2020; Louis et al. 2025; Camphuis et al. 2025), baryon acoustic oscillations (BAOs; Alam et al. 2021; Abdul Karim et al. 2025), redshift space distortions (RSDs; Alam et al. 2021), and measurements of the expansion history with observations of Type Ia supernovae (SN Ia; Scolnic et al. 2022; Brout et al. 2022). Under the assumption of baryonic and CDM components, along with a dark energy component expressed by a constant energy density with negative pressure, which are connected via a spatially flat gravitational framework, it provides an exceptionally simple model describing the evolution of the Universe on large scales with just six parameters. Despite its success, the ΛCDM model faces a number of challenges both on the observational side, such as discrepancies in the expansion rate between probes of the early and late Universe (Di Valentino et al. 2025), and on the theoretical side, such as the lack of an explanation for the value of the cosmological constant by fundamental physics. This makes the study of cosmological models beyond ΛCDM an active field of research.

Studies of extended cosmological models beyond ΛCDM typically include variations of the neutrino mass (νΛCDM), a deviation from a spatially flat Universe (Ω_KCDM), dark energy models with a dynamical equation of state (wCDM), and modifications to the theory of gravity on large scales. In our companion paper Reischke et al. (2025a), we derived constraints on νΛCDM, Ω_KCDM, and wCDM models. In this work, we focus on scalar-tensor modified gravity theories, which introduce a scalar degree of freedom in addition to the tensor metric degrees of freedom. Here, the scalar field replaces the cosmological constant as the driver of the accelerated expansion of the Universe and also influences the growth of structure through its dynamics. By construction, these theories typically mimic ΛCDM or minimal extensions thereof on the background level, given the strong observational constraints on the expansion history of the Universe. In particular, we derived constraints on the Horndeski class of modified gravity models (Horndeski 1974; Nicolis et al. 2009; Deffayet et al. 2011; Kobayashi et al. 2011), which encompasses many scalar-tensor theories that have been studied in the literature, such as f(R) gravity, Brans-Dicke, Galileons, and Quintessence.

A popular strategy in studies of Horndeski gravity is an effective field theory (EFT) approach that fully describes the evolution of perturbations at the linear level on top of a given background cosmological model. This approach was adopted to forecast the possible constraints with Stage-IV surveys (Spurio Mancini et al. 2018; Reischke et al. 2019). Moreover, Spurio Mancini et al. (2019) carried out a combined cosmic shear, galaxy clustering, and galaxy-galaxy lensing analysis using data from the third data release of the Kilo-Degree Survey (KV450; Hildebrandt et al. 2020; Wright et al. 2020) within a modified gravity framework. Here, we follow up on this work with significant improvements in the methodology as well as the data. Our new analysis uses, for the first time, a parametrisation of Horndeski gravity using an inherently stable parameter basis. This is enabled by the MOCHI_CLASS code (Cataneo & Bellini 2024, hereafter CB24), an extension of the Einstein-Boltzmann solvers HI_CLASS (Zumalacárregui et al. 2017; Bellini et al. 2020) and CLASS (Lesgourgues 2011; Blas et al. 2011), which features a parametrisation of the Horndeski models in terms of the stable basis functions proposed by Kennedy et al. (2018). In contrast to other EFT parametrisations, this basis automatically selects models that are free from ghost and gradient instabilities by design. Furthermore, it activates modifications to general relativity (GR) only after a user-defined redshift while fixing the early time evolution to a chosen background cosmological model, such as ΛCDM. Additionally, we derived constraints on modified gravity models with an effective dark energy equation of state parametrised by w₀ and w_a.

We made use of cosmic shear data from the final, fifth data release of the Kilo-Degree Survey (KiDS-Legacy; Wright et al. 2024). Based on these data, Wright et al. (2025b) and Stölzner et al. (2025) derived constraints on flat ΛCDM with KiDS-Legacy data. The amount of matter clustering, described by the structure growth parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, was found to be consistent with observations of the CMB by Planck (Planck Collaboration VI 2020), resolving the tension seen in previous KiDS-Planck comparisons (Asgari et al. 2021; Heymans et al. 2021). This enabled us to perform a joint analysis of KiDS and Planck. Additionally, we employed BAO measurements from the Dark Energy Spectroscopic Instrument (DESI; Abdul Karim et al. 2025) and RSD measurements from the extended Baryon Oscillation Spectroscopic Survey (eBOSS Alam et al. 2021).

This paper is structured as follows: In Sect. 2, we describe our Horndeski modelling framework and in Sect. 3, we briefly summarise the weak lensing, BAO, RSD, and CMB data employed in this work. Sect. 4 presents the results of our analysis before concluding in Sect. 5. In Appendix A, we provide constraints on all cosmological parameters and Appendix B presents a study of the impact of the CMB lensing anomaly on modified gravity constraints. In Appendix C, we show the prior on the derived modified gravity parameters for each parametrisation considered in this work.

2. Methodology

2.1. Horndeski theory

We focus on the Horndeski theory of gravity, which is the most general scalar-tensor theory in four dimensions with second-order equations of motion. Its action is expressed as (Horndeski 1974; Nicolis et al. 2009; Deffayet et al. 2011; Kobayashi et al. 2011)

$$\begin{array}{c}\hfill S={\displaystyle \int {\mathrm{d}}^4}x\sqrt{-g}[{\displaystyle \sum _{i=2}^5}\frac{1}{8\pi G_N}{\mathcal{L}}_i[g_{\mu \nu },\varphi ]+{\mathcal{L}}_{\mathrm{m}}[g_{\mu \nu },{\psi }_{\mathrm{m}}]],\end{array}$$(1)

with

$$\begin{array}{cc}\hfill {\mathcal{L}}_2& =G_2(\varphi ,X),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\mathcal{L}}_3& =-G_3(\varphi ,X)\square \varphi ,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathcal{L}}_4& =G_4(\varphi ,X)R+G_{4X}(\varphi ,X)[{(\square \varphi )}^2-{\varphi }_{;\mu \nu }{\varphi }^{;\mu \nu }],\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\mathcal{L}}_5& =G_5(\varphi ,X)G_{\mu \nu }{\varphi }^{;\mu \nu }-\frac{1}{6}{G}_{5X}(\varphi ,X),\hfill \\ \hfill & \times [{(\square \varphi )}^3+2{\varphi }_{;\mu }^{\nu }{\varphi }_{;\nu }^{\alpha }{\varphi }_{;\alpha }^{\mu }-3{\varphi }_{;\mu \nu }{\varphi }^{;\mu \nu }\square \varphi ],\hfill \end{array}$$(5)

where g_μν is the metric tensor with its determinant denoted by g, ℒ_m is the matter Lagrangian, which is a function of the metric and the minimally coupled matter fields, ψ_m, R is the Ricci scalar, and G_i(ϕ,X) are arbitrary functions of the scalar field, ϕ, and its kinetic term, $X=-\frac{1}{2}{\varphi }_{;\mu }{\varphi }^{;\mu }$. Here, ϕ_;μ denotes the covariant derivative, G_iX = ∂G_i/∂X indicates the partial derivative, □ϕ = g^μν∇_μ∇_νϕ denotes the d’Alembert operator, and we use the Einstein summation convention. The specific choice of G_i(ϕ,X) corresponds to a selection of a single modified gravity model contained within the general class of Horndeski models (see e.g. Kobayashi 2019, for a review).

For any given cosmological background, perturbations can be parametrised in an EFT approach (Gubitosi et al. 2013; Bloomfield et al. 2013; Gleyzes et al. 2013; Bloomfield 2013; Piazza et al. 2014; Gergely & Tsujikawa 2014). In particular, as shown in Gleyzes et al. (2013) and Bellini & Sawicki (2014), the evolution of linear perturbations can be fully described by four arbitrary functions of time, α_i(t). In practice, the α_i(t) functions are commonly chosen to be proportional to the time-dependent fractional energy density of dark energy, Ω_DE(t), so that

$$\begin{array}{c}\hfill {\alpha }_i(t)={\widehat{\alpha }}_i{\mathrm{\Omega }}_{\mathrm{DE}}(t),\end{array}$$(6)

with the proportionality coefficients, ${\widehat{\alpha }}_i$. This choice selects a sub-class of Horndeski theories, which enforce GR at early times with modifications to GR becoming relevant as soon as dark energy contributes significantly to the background energy density, assuming the same time-dependence for all functions. However, as discussed by Gleyzes (2017), this simple parametrisation is sufficient to provide information on a large fraction of the Horndeski theory space. Therefore, it has been adopted in previous studies of modified gravity with KiDS data (Spurio Mancini et al. 2019).

Furthermore, a large range of popular dark energy models can be described by adopting particular parametrisations of the α functions (see Table 1 in Bellini & Sawicki 2014). Here, we briefly summarise the α functions and refer to Bellini & Sawicki (2014) for a complete description. The kineticity term, α_K, describes the kinetic energy of the scalar perturbation, which does not enter the equations of motion in the sub-horizon, quasi-static approximation; therefore, it is unconstrained by large-scale structure probes (Bellini et al. 2016; Alonso et al. 2017). The tensor speed excess, α_T, parametrises a deviation between the speed of light and the propagation speed of gravitational waves. This parameter is tightly constrained by the detection of the binary neutron star merger GW170817 and the corresponding gamma-ray burst GRB170817A, which determined the tensor speed to be close to the speed of light and thus constrained α_T to be vanishingly small at the present time (Abbott et al. 2017a,b; Amendola et al. 2018; Baker et al. 2017; Bettoni et al. 2017; Ezquiaga & Zumalacárregui 2017; Lombriser & Lima 2017; Sakstein & Jain 2017). However, outside the regime described by EFT, there remains the possibility of a different propagation speed of gravitational waves (Baker et al. 2023). The evolution of the run rate of the effective Planck mass, M_*², is described by

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{M}}=\frac{\mathrm{d}ln{M}_{\ast }^2}{\mathrm{d}lna},\end{array}$$(7)

which parametrises the rate of change of the cosmological strength of gravity. Finally, the braiding term, α_B, describes the mixing of the scalar field with the metric kinetic term, which leads to a coupling of scalar and metric perturbations.

One particular challenge for Horndeski gravity models comes from violations of the stability conditions, for which the perturbations cause the cosmological background to be non-viable. Physical instabilities can be classified as: (i) gradient instabilities, which occur when the square of the speed of sound becomes negative; (ii) ghost instabilities due to a wrong sign of the kinetic term of background perturbations; and (iii) tachyonic instabilities, arising from a negative effective mass squared of the scalar field perturbation (Hu et al. 2014a; Bellini & Sawicki 2014). Additionally, mathematical instabilities in the equation of motion for the scalar field fluctuations can lead to exponentially growing perturbations (Hu et al. 2014b). As a consequence, sampling the parameter space in Markov chain Monte Carlo (MCMC) analyses can become highly inefficient, with only a small fraction of samples being generated in stable regions of parameter space (Pèrenon et al. 2015). To circumvent these issues, Kennedy et al. (2018), Lombriser et al. (2019) introduced an inherently stable basis. In this parametrisation, the α functions are replaced by three derived functions and a constant. These are comprised by the de-mixed kinetic term for the scalar field perturbation,

$$\begin{array}{c}\hfill D_{\mathrm{kin}}={\alpha }_{\mathrm{K}}+\frac{3}{2}{\alpha }_{\mathrm{B}}^2,\end{array}$$(8)

its effective sound speed squared,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c_{\mathrm{s}}^2=\frac{1}{D_{\mathrm{kin}}}[(2-{\alpha }_{\mathrm{B}})(-\frac{H^{\prime }}{a{H}^2}& +\frac{1}{2}{\alpha }_{\mathrm{B}}+{\alpha }_{\mathrm{M}})\hfill \\ \hfill & -\frac{3({\rho }_{\mathrm{M}}+p_{\mathrm{M}})}{H^2{M}_{\ast }^2}+\frac{{\alpha }_{\mathrm{B}}^{\prime }}{\mathit{aH}}],\hfill \end{array}\end{array}$$(9)

the effective Planck mass, M_*², and an initial condition on the braiding, α_B0 (CB24). Here, H is the Hubble parameter, ρ_M and p_M denote the matter energy density and pressure, respectively, and primes denote the derivative with respect to conformal time, defined as τ = ∫_z^∞dz/H(z). The evolution of α_B then follows from the ordinary differential equation given in Eq. (27) in CB24.

Adopting D_kin, M_*², and c_s as basis functions, we followed the approach of CB24 and approximated the stable basis functions in the matter-dominated era either as power laws or constant values which approach the GR limit. For α_B ≈ 0, the braiding term can be approximated by the linearised early-time solution, ${\alpha }_{\mathrm{B}}^{(\mathrm{lin})}$ given in Eq. (41) in CB24. By minimising the difference $\mathrm{\Delta }{\alpha }_{\mathrm{B}}={\alpha }_{\mathrm{B}}^{(\mathrm{ode})}-{\alpha }_{\mathrm{B}}^{(\mathrm{lin})}$ in the matter-dominated era we then inferred the value of α_B, 0 for a given set of modified gravity parameters. We note that in synchronous gauge (i.e. the gauge adopted by MOCHI_CLASS), α_B = 2 corresponds to a singular point in parameter space (Lagos et al. 2018; Noller & Nicola 2019; Traykova et al. 2021). Therefore, this value sets an upper bound on α_B at any time during the evolution of perturbations.

While the absence of gradient and ghost instabilities is essential when constructing a viable model, tachyonic instabilities can still produce viable models (Zumalacárregui et al. 2017; Frusciante et al. 2019; Gsponer & Noller 2022). Therefore, we require

$$\begin{array}{c}\hfill D_{\mathrm{kin}}>0,M_{\ast }^2>0,c_{\mathrm{s}}^2>0\end{array}$$(10)

to select viable modified gravity models that are free from gradient and ghost instabilities. In practice, we replaced the Planck mass, M_*², with its deviation from the fiducial GR value, defined by ΔM_*² ≡ M_*² − 1, when sampling the parameter space. As shown in CB24, the stable parametrisation is of particular advantage in efficiently sampling viable theoretical models, which the α functions might incorrectly classify as unstable due to inaccurate numerical cancellations in the calculation of c_s² via Eq. (9), and in improving the numerical stability. Therefore, we employed the stable basis as our fiducial parametrisation of modified gravity in this work and computed the derived α functions for comparison with previous works. As functional form for the three free functions of the stable parameter basis, we adopted an approach that resembles the one proposed by Lombriser et al. (2019) and generalised the proportionality of the α functions to the dark energy density, given in Eq. (6), to the stable parameter basis. In particular, we parametrised the free functions D_kin(a) and ΔM_*²(a) via

$$\begin{array}{c}\hfill X_i(a)={\widehat{X}}_i\frac{{\mathrm{\Omega }}_{\mathrm{DE}}(a)}{{\mathrm{\Omega }}_{\mathrm{DE}}(a=1)},\end{array}$$(11)

where the amplitudes ${\widehat{X}}_i$ are the sampling parameters in our likelihood analysis, which will be denoted in the following sections by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$. We assume the effective sound speed, denoted by ${\widehat{c}}_{\mathrm{s}}^2$, to be constant over time in order to simplify the parameter space in the likelihood analysis due to the degeneracy between D_kin and c_s², which can be inferred from the effective Newtonian coupling and the gravitational slip in the quasi-static approximation, which are given below. Additionally, we compared our results with an analysis that uses the original alpha functions, instead of the stable basis.

The quasi-static approximation plays an important role in the study of modified gravity models. On scales much smaller than the sound horizon (i.e. c_s²k² ≫ a²H²), the time derivative of the scalar field perturbation can be neglected, and the effective Newtonian coupling and the gravitational slip are given by (CB24)

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{QSA}}(k,a)& =\frac{1}{M_{\ast }^2}\frac{{\mu }_{\mathrm{p}}+k^2{c}_{\mathrm{sN}}^2{M}_{\ast }^2{\mu }_{\infty }/a^2{H}^2}{{\mu }_{\mathrm{p}}+k^2{c}_{\mathrm{sN}}^2/a^2{H}^2},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\gamma }_{\mathrm{QSA}}(k,a)& =\frac{{\mu }_{\mathrm{p}}+k^2{c}_{\mathrm{sN}}^2{M}_{\ast }^2{\mu }_{\mathrm{Z},\infty }/a^2{H}^2}{{\mu }_{\mathrm{p}}+k^2{c}_{\mathrm{sN}}^2{M}_{\ast }^2{\mu }_{\infty }/a^2{H}^2},\hfill \end{array}$$(13)

with μ_∞, μ_Z, ∞, and μ_p given in Appendix D in CB24 and c_sN² ≡ D_kinc_s². This approximation allows for significant computational improvements and is activated selectively based on scale and time (see Sect. 2.2 in CB24, for more details about the implementation in the modelling code).

2.2. Weak lensing

In conformal Newtonian gauge, adopting a spatially flat background, and considering linear perturbations, the line element is expressed as

$$\begin{array}{c}\hfill \mathrm{d}{\mathit{s}}^2=-(1+2\mathrm{\Psi })c^2\mathrm{d}{t}^2+a^2(1-2\mathrm{\Phi })\mathrm{d}{\mathit{x}}^2,\end{array}$$(14)

where the Newtonian potential, Ψ, and the spatial curvature potential, Φ, are the Bardeen potentials (Bardeen 1980), which in GR and in the absence of anisotropic stress fulfil the condition Φ = Ψ. In a modified gravity scenario, on sub-horizon scales and adopting the quasi-static approximation, the linear Poisson equation in Fourier space becomes (Amendola et al. 2008):

$$\begin{array}{c}\hfill -k^2\mathrm{\Psi }=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^{2}a}{\delta }_{\mathrm{m}}\mu (k,a)\end{array}$$(15)

with

$$\begin{array}{c}\hfill \mathrm{\Phi }=\gamma (k,a)\mathrm{\Psi },\end{array}$$(16)

where δ_m denotes the matter density contrast and we introduced the effective Newtonian coupling, μ, which parametrises deviations in the gravitational constant, G, due to additional fields or modified dynamics, and the gravitational slip, γ, which describes the strength of anisotropic stress. In the GR limit, the effective Newtonian coupling and the gravitational slip become equal to unity.

Cosmic shear is sensitive to the Weyl potential Φ_W = (Φ + Ψ)/2, for which Eqs. (15) and (16) yield

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{W}}=-\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2{k}^{2}a}{\delta }_{\mathrm{m}}\mathrm{\Sigma }(k,a).\end{array}$$(17)

Here, Σ(k, a)≡μ(k, a)(1 + γ(k, a))/2 parametrises modifications to the GR lensing potential. Since, on small scales, there is strong observational evidence for GR (see e.g. Will 2014, for a review), we incorporated a phenomenological Vainshtein screening mechanism (Vainshtein 1972) in the non-linear regime to recover GR at high-k. We define

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{NL}}(k,a)=1+\mathcal{S}(k,a)[{\mathrm{\Sigma }}_{\mathrm{L}}(k,a)-1],\end{array}$$(18)

where the screening factor is given by

$$\begin{array}{c}\hfill \mathcal{S}(k,a)=exp[-{\left(\frac{k}{k_{\mathrm{V}}}\right)}^{n(a)}],\end{array}$$(19)

which converges to zero for k ≫ k_V, recovering the GR value of Σ_NL = 1. Here, k_V and n(a) define the screening scale and the slope of the screening, respectively. In practice, we fixed the slope of the screening scale to n = 1 and treated k_V as a free model parameter. This choice was based on tests on the normal-branch Dvali-Gabadadze-Porrati (nDGP; Dvali et al. 2000; Deffayet 2001) braneworld model, which showed that it provides a good approximation when applied to the nDGP power spectrum with no indication for a significant evolution of n over the redshift range that is relevant for this work. The power spectrum of the Weyl potential is then given by

$$\begin{array}{c}\hfill P_{\mathrm{W},\mathrm{NL}}^{\mathrm{MG}}(k,a)={\left(\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2{k}^{2}a}\right)}^2{\mathrm{\Sigma }}_{\mathrm{NL}}^2{P}_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}(k,a),\end{array}$$(20)

where $P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}$ denotes the non-linear matter power spectrum in modified gravity. Under the assumption of the extended Limber approximation (Kaiser 1992; LoVerde & Afshordi 2008; Kilbinger et al. 2017), we find the gravitational lensing power spectrum to be

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill C_{\mathrm{GG}}^{\mathit{ij}}(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}}\mathrm{d}\chi \frac{W_{\mathrm{G}}^i(\chi )W_{\mathrm{G}}^j(\chi )}{f_{\mathrm{K}}^2(\chi )}& {\mathrm{\Sigma }}_{\mathrm{NL}}^2(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},z(\chi ))\hfill \\ \hfill & \times P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},z(\chi )),\hfill \end{array}\end{array}$$(21)

where i and j denote tomographic bins, f_K, χ, and χ_H are the comoving angular diameter distance, the comoving radial distance, and the comoving horizon distance, respectively. The weak lensing kernel, $W_{\mathrm{G}}^i(\chi )$, is defined as

$$\begin{array}{c}\hfill W_{\mathrm{G}}^i(\chi )=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{f_{\mathrm{K}}(\chi )}{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}}\mathrm{d}{\chi }^{\prime }{n}_{\mathrm{S}}^i\left({\chi }^{\prime }\right)\frac{f_{\mathrm{K}}({\chi }^{\prime }-\chi )}{f_{\mathrm{K}}\left({\chi }^{\prime }\right)},\end{array}$$(22)

where n_Sⁱ(χ′) is the distribution of source galaxies per comoving distance. In practice, we obtain the theoretical prediction for the gravitational lensing power spectrum via Eq. (21) and compute Σ_L directly from the transfer functions T_Φ, T_Ψ, and T_m via

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{L}}=-\frac{c^2{k}^{2}a}{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}\frac{T_{\mathrm{\Phi }}+T_{\mathrm{\Psi }}}{T_{\mathrm{m}}},\end{array}$$(23)

which do not require the assumption of the quasi-static approximation.

To recover GR at small scales, we adopted an approach inspired by the reaction framework of Cataneo et al. (2019) and defined the non-linear matter power spectrum as

$$\begin{array}{c}\hfill P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}(k,a)=\mathcal{R}(k,a)P_{\mathrm{m},\mathrm{NL}}^{\mathrm{pseudo}}(k,a),\end{array}$$(24)

where $P_{\mathrm{m},\mathrm{NL}}^{\mathrm{pseudo}}(k,a)$ is the non-linear power spectrum in a reference cosmology, called the pseudo cosmology, which we compute in the Horndeski modified gravity framework. The reaction function ℛ is then chosen so that $P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}$ converges to the GR non-linear power spectrum on small scales. In general, we modelled the background evolution of dark energy for the linear matter power spectrum via the Chevallier-Polarski-Linder (CPL) parametrisation (Chevallier & Polarski 2001; Linder 2003), in which the dark energy equation of state is given by w(a) = w₀ + w_a(1 − a), and derived its non-linear evolution assuming a ΛCDM background via the augmented halo model HMCODE2020 (Mead et al. 2021). The phenomenological reaction is defined as

$$\begin{array}{c}\hfill \mathcal{R}(k,a)={[(1-\mathcal{S}(k,a))\sqrt{\frac{P_{\mathrm{m},\mathrm{NL}}^{w_0{w}_a}}{P_{\mathrm{m},\mathrm{NL}}^{\mathrm{pseudo}}}}+\mathcal{S}(k,a)]}^2,\end{array}$$(25)

where P_m, NL^w₀w_a denotes the non-linear matter power spectrum in GR assuming a w₀w_aCDM cosmology for the linear power spectrum and its non-linear correction. In this parametrisation, we find $P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}\approx P_{\mathrm{m},\mathrm{NL}}^{w_0{w}_a}$ for k ≫ k_V and $P_{\mathrm{m},\mathrm{NL}}^{\mathrm{MG}}\approx P_{\mathrm{m},\mathrm{NL}}^{\mathrm{pseudo}}$ for k ≪ k_V, and thus we recover GR on small scales.

The primary cosmic shear observable is the ellipticity correlation between galaxy pairs, given by the sum of the gravitational lensing power spectrum (GG), the intrinsic alignment of galaxies (II), and their cross term (GI):

$$\begin{array}{c}\hfill C_{ϵϵ}^{\mathit{ij}}(\ell )=C_{\mathrm{GG}}^{\mathit{ij}}(\ell )+C_{\mathrm{II}}^{\mathit{ij}}(\ell )+C_{\mathrm{GI}}^{\mathit{ij}}(\ell )+C_{\mathrm{IG}}^{\mathit{ij}}(\ell ).\end{array}$$(26)

In weak lensing studies, intrinsic alignments (IAs) are commonly treated as additional effects contaminating the observed cosmic shear signal. As shown in Reischke et al. (2022), IAs by themselves probe the gravitational potential and therefore can provide constraints on modified gravity. However, Reischke et al. (2022) showed that even the statistical power of upcoming Stage-IV surveys is not sufficient to derive constraints on modified gravity with IAs without further improvements on the IA modelling site. Therefore, we neglected the impact of modified gravity on the IA power spectra and adopted the fiducial expressions for $C_{\mathrm{II}}^{\mathit{ij}}(\ell )$ and $C_{\mathrm{GI}}^{\mathit{ij}}(\ell )$, given in Eqs. (4) and (5) in Wright et al. (2025b), which we computed from the non-linear matter power spectrum in modified gravity.

2.3. Baryon acoustic oscillations and redshift space distortions

Baryon acoustic oscillations are imprints of sound waves in the early Universe, which cause fluctuations in the matter density. Spectroscopic galaxy surveys observe the BAO feature both in the line-of-sight direction and the transverse direction and therefore provide measurements of the Hubble distance,

$$\begin{array}{c}\hfill D_{\mathrm{H}}(z)=\frac{c}{H(z)},\end{array}$$(27)

and the transverse comoving distance,

$$\begin{array}{c}\hfill D_{\mathrm{M}}(z)=\frac{c}{H_0\sqrt{{\mathrm{\Omega }}_{\mathrm{K}}}}sinh\left(\sqrt{{\mathrm{\Omega }}_{\mathrm{K}}}{\displaystyle {\int }_0^z}\mathrm{d}{z}^{\prime }\frac{H_0}{H(z^{\prime })}\right),\end{array}$$(28)

relative to the comoving sound horizon at the drag epoch, r_d. Here, H(z) follows from the Friedmann equation for the assumed background cosmology, which we adopted to be ΛCDM. The exception is given in Sect. 4.4, where we consider a dynamical dark energy background. Since the BAO measurements constrain background quantities, we did not expect them to provide constraints on our modified gravity model, which parametrises perturbations around the background. However, reconstructed BAO measurements are commonly derived under the assumption of GR, which can potentially bias the position of the BAO peak. As was shown by Bellini & Zumalacárregui (2015) and Pan et al. (2024), this shift is well below the statistical error budget, even for Stage-IV surveys.

Redshift space distortions (RSDs), on the other hand, are a probe of the growth of structure in the Universe and thus provide constraints on modified gravity models. Typically, RSD measurements probe fσ₈(z), where f denotes the growth rate of linear perturbations and σ₈ is the root-mean-square of matter density fluctuations in spheres of 8 Mpc/h radius. The growth rate is defined as

$$\begin{array}{c}\hfill f=\frac{\mathrm{d}lnD}{\mathrm{d}lna},\end{array}$$(29)

where the growth factor, D, in GR can be computed by solving (Heath 1977):

$$\begin{array}{c}\hfill D^{″}+aH{D}^{\prime }+\frac{3}{2}{a}^2{\rho }_{\mathrm{M}}D=0,\end{array}$$(30)

as implemented, for example, in the Boltzmann code CLASS. However, this equation no longer holds in modified gravity due to modifications to the Poisson equation (Zhao et al. 2009). Therefore, we computed the effective, scale-independent fσ₈(z) via

$$\begin{array}{c}\hfill f{\sigma }_8(z)=\frac{\mathrm{d}{\sigma }_8}{\mathrm{d}lna},\end{array}$$(31)

which does not rely on GR-specific assumptions in the computation of fσ₈(z). We note that the linear growth rate in general is scale-dependent. However, as discussed in Appendix A of Noller & Nicola (2019), the growth rate is found to be relatively scale-independent, except on ultra-large scales, for the sub-class of models analysed in this work. We verified explicitly that this finding holds for the stable parameter basis adopted in this work. Moreover, RSD measurements from galaxy surveys are typically based on GR assumptions and are validated using GR mock catalogues. In this work, we followed the approach of Alam et al. (2021) and assumed the adopted RSD measurements to be model-independent, which is justified given the approximately 10% level of precision for the RSD data. We note that this assumption will not hold for Stage-IV surveys, for which the statistical uncertainties become comparable to the modelling bias induced when fitting modified gravity cosmologies with GR-based RSD templates. Therefore, self-consistent RSD analyses under the assumed modified gravity model will be required (Bose et al. 2017).

2.4. Cosmic microwave background

Cosmic microwave background anisotropies are a key probe of the cosmological model. In combination with probes of the low-redshift Universe, CMB measurements essentially fix the cosmological model at early times (Planck Collaboration VI 2020), which allows us to derive constraints on the evolution of dark energy or modified gravity at late times. The most dominant effects of modified gravity on the CMB power spectrum are due to changes of the lensing potential (Acquaviva & Baccigalupi 2006; Carbone et al. 2013) and the integrated Sachs-Wolfe (ISW) effect (Sachs & Wolfe 1967). Since the ISW effect is a secondary anisotropy peaking at large angular scales, which are limited by cosmic variance, it is challenging to detect via direct measurements of the CMB temperature power spectrum. Therefore, it is typically constrained through cross-correlation measurements between the CMB and tracers of the large-scale structure (Boughn et al. 1998; Boughn & Crittenden 2001; Giannantonio et al. 2006; Stölzner et al. 2018). In this work, we account for the ISW effect in the theoretical prediction of the CMB power spectrum and leave a study of modified gravity models through ISW cross-correlations for future work.

As discussed in detail, for example, in Planck Collaboration VIII (2020), CMB lensing causes a smoothing of the CMB power spectrum, which can be computed from the theoretical prediction for the unlensed CMB power spectrum and the CMB lensing potential power spectrum (Seljak 1996; Lewis & Challinor 2006). Internal consistency tests of the CMB typically quantify whether the smoothing effect matches the expectation by allowing for a scaling of the lensing power spectrum via an amplitude parameter A_L, which in GR is expected to be equal to one (Calabrese et al. 2008). A deviation from the expected value then indicates either unaccounted systematics in the data or new physics beyond ΛCDM. This parameter is of particular importance when studying modified gravity models, which were found to be degenerate with lensing-induced smoothing effects, since CMB lensing is sensitive to the Weyl potential (Pogosian et al. 2022). Previous studies based on the public Planck PR3 plik likelihood found a preference for A_L > 1 at the 2σ level (Planck Collaboration VI 2020), while analyses based on the more recent Planck PR4 data reported A_L to be consistent with the ΛCDM expectation within 0.7σ (Tristram et al. 2024). In this work, we therefore chose to adopt A_L as a sampling parameter, which we marginalised over in our likelihood analysis. In Appendix B, we explore the impact of a fixed A_L = 1 on our modified gravity constraints.

2.5. Analysis setup

We computed the theoretical prediction in the Horndeski model of modified gravity for the linear matter power spectrum, its transfer functions, the CMB power spectrum, and the growth rate via the public MOCHI_CLASS¹ code (CB24) and modelled the effects of baryonic processes on the matter power spectrum with the augmented halo model HMCODE2020 (Mead et al. 2021). We performed the sampling of the parameter space via the NAUTILUS² sampler (Lange 2023), which is interfaced within the COSMOSIS³ framework (Zuntz et al. 2021). A list of cosmological and nuisance parameters and their priors is provided in Table 1. In the combined analysis of KiDS, DESI, eBOSS, and Planck data, we treated each dataset as independent. Therefore, we computed the joint likelihood by multiplying the individual likelihoods of each experiment.

We considered two analysis setups for the stable parametrisation of Horndeski gravity: First, we varied the de-mixed kinetic term of the scalar field perturbation and the deviation of the effective Planck mass from its fiducial value, ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, while keeping the effective sound speed, ${\widehat{c}}_{\mathrm{s}}^2$, fixed. Second, we varied ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, fixing the Planck mass to its fiducial value. Additionally, we considered the parametrisation in terms of the α functions as a direct comparison with previous studies. In this parametrisation, we sampled the braiding term and the run rate of the Planck mass, ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, for fixed values of the kineticity and tensor speed excess, as discussed in Sect. 2.1. We note that although the kineticity was shown to be unconstrained by large-scale structure probes, it enters the stability conditions and therefore places a theoretical bound on the underlying model. We explicitly tested the impact of the implicit prior by conducting an analysis for ${\widehat{\alpha }}_{\mathrm{K}}=0.01$ in addition to our fiducial choice of ${\widehat{\alpha }}_{\mathrm{K}}=1.0$. We found both analyses to be in very good agreement and therefore conclude that our cosmological constraints are not impacted by the chosen value of ${\widehat{\alpha }}_{\mathrm{K}}$. Furthermore, following CB24, we set the transition redshift for enabling the GR approximation scheme in the stable parametrisation in MOCHI_CLASS to z = 99, which ensures that the switch occurs deep in the matter-dominated regime (see Sect. 3.3 in CB24, for more details about the GR approximation scheme).

The main focus of this work is to study the impact of modified gravity on perturbations under the assumption of a ΛCDM background cosmology. Therefore, we considered a background dark energy model with a constant equation of state parameter, w(a) =  − 1, which is a common approach in modified gravity analyses (see for example Bellini et al. 2016; Noller & Nicola 2019; Spurio Mancini et al. 2019; Shah et al. 2026). This choice is motivated by our selection of cosmological datasets, which independently probe the background evolution and the structure growth. In Sect. 4.4, we consider a more general analysis setup, in which we free up the evolution of the background.

3. Data

We employed weak lensing measurements from the fifth and final data release of the Kilo-Degree Survey (Wright et al. 2024), dubbed KiDS-Legacy. The KiDS and the complementary VISTA Kilo-Degree Infrared Galaxy Survey (VIKING; Edge et al. 2013) combine optical and near-infrared imaging in nine photometric bands. We adopted the fiducial data vector of Wright et al. (2025b), in which the KiDS-Legacy lensing sample was divided into six approximately equi-populated bins between 0.1 < z_B ≤ 2.0 via their best-fit photometric redshift estimate, z_B, from the template-fitting code BPZ (Benítez 2000). The redshift distributions per tomographic bin were calibrated with a direct calibration method with deep spectroscopic surveys using self-organising maps, as outlined in Wright et al. (2025a). For the uncertainty modelling (Reischke et al. 2025b), we used the fiducial covariance and did not recompute it for the new best fit values, which was shown in Wright et al. (2025b) to have a negligible impact on the inferred cosmological constraints. The fiducial KiDS-Legacy cosmic shear analyses (Wright et al. 2025b; Stölzner et al. 2025) found the resulting cosmological constraints to be in agreement with CMB data as well as with other low redshift probes, such as BAO and RSD measurements and other weak lensing experiments. This enables a combined analysis of KiDS-Legacy with a broad variety of cosmological probes.

Our analysis pipeline is based on the fiducial KiDS-Legacy pipeline (Wright et al. 2025b), which is implemented in the public COSMOPIPE⁴ infrastructure. We employed complete orthogonal sets of E/B-integrals (COSEBIs, Schneider et al. 2010; Asgari et al. 2012) as our cosmic shear summary statistic and adopted the fiducial mass-dependent IA model (see Sect. 2.3.4 and Appendix B in Wright et al. 2025b), which parametrises the amplitude and the slope of the IA mass scaling via two nuisance parameters, A_IA and β, with a bivariate Gaussian prior inferred from the joint posterior in Fortuna et al. (2025) and accounts for the halo mass per tomographic bin via a multivariate Gaussian prior. Additionally, we adopted the fiducial KiDS-Legacy redshift distributions from Wright et al. (2025a) and the corresponding multivariate Gaussian prior on the shift in the mean of the redshift distribution per tomographic bin.

We adopted BAO measurements from the second data release of the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016, 2022). DESI targets four classes of extragalactic objects, covering a wide range of redshifts: a bright galaxy sample between 0.1 < z < 0.4 (Hahn et al. 2023), two samples of luminous red galaxies (LRGs) between 0.4 < z < 0.6 and 0.6 < z < 0.8 (Zhou et al. 2023), an emission line galaxy (ELG) sample between 1.1 < z < 1.6 (Raichoor et al. 2023), a combined ELG and LRG sample between 0.8 < z < 1.1, and a quasar sample between 0.8 < z < 2.1 (Chaussidon et al. 2023). Additionally, DESI obtained a BAO measurement from the Lyman-α forest of high-redshift quasars between 1.77 < z < 4.16. In this work, we made use of the DESI DR2 BAO likelihood of Abdul Karim et al. (2025), which is publicly available in the COSMOSIS standard library⁵.

In addition to DESI BAO observations, we employed RSD measurements from the Sloan Digital Sky Survey’s extended Baryon Oscillation Spectroscopic Survey (eBOSS). In particular, we made use of the ‘RSD-only’ data from Alam et al. (2021), which were inferred by marginalising over D_H and D_M in the full-shape measurements of the matter power spectrum, yielding measurements on the structure growth parameter fσ₈(z) from the SDSS DR7 (Ross et al. 2015; Howlett et al. 2015) and BOSS DR12 (Alam et al. 2017) main galaxy samples and the eBOSS DR16 LRG (Gil-Marín et al. 2020; Bautista et al. 2021), ELG (Tamone et al. 2020; Raichoor et al. 2021), and quasar (Neveux et al. 2020; Hou et al. 2021) samples. Following Alam et al. (2021), who found the correlation between eBOSS samples to be negligibly small, we assumed the RSD measurements to be independent and adopted an individual Gaussian likelihood for each tracer as listed in table III of Alam et al. (2021).

Finally, we made use of measurements of the CMB temperature and polarisation power spectra from Planck Collaboration VI (2020). In particular, we employed the public Planck high-ℓ TTTEEE, low-ℓ TT, and low-ℓ EE likelihoods (Planck Collaboration V 2020). While the fiducial Planck likelihood features a vector of 47 nuisance parameters accounting for astrophysical foregrounds and instrumental characteristics, we made use of the Plik_lite likelihood, which is pre-marginalised over all nuisance parameters except for the Planck absolute calibration and the lensing anomaly parameter. Although the pre-marginalisation was performed under the assumption of ΛCDM, Noller & Nicola (2019) showed that the inferred constraints on Horndeski parameters are consistent with the full likelihood. We explicitly confirmed that this finding holds for the Horndeski models considered in this work and thus employed the pre-marginalised likelihood for all analysis setups discussed in Sect. 4.

4. Results

4.1. Constraints on S₈

Our main results in the S₈-Ω_m plane, inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets, are shown in Fig. 1 for the three parametrisations of modified gravity considered in this work. Table 2 presents the marginal mode and the 68% highest posterior density interval (HPDI) for S₈ and Ω_m, the χ² of our cosmological model evaluated at the maximum a posteriori (MAP), the model probability to exceed (PTE) at the MAP, the difference in χ² between ΛCDM and modified gravity models, and the N_σ preference level for modified gravity inferred in a Bayesian suspiciousness test (Handley & Lemos 2019). In Appendix A, we provide the best-fit values and posterior contours for all model parameters.

Fig. 1. Posterior distribution in the S8-Ωm plane for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contours illustrate the posterior inferred in a Horndeski model parametrised by the stable basis parameters ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, and the green contour shows constraints from sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$. The orange contour shows the posterior in a Horndeski model using the ${\widehat{\alpha }}_{\mathrm{B}}$-${\widehat{\alpha }}_{\mathrm{M}}$ parametrisation. The black contour corresponds to the fiducial ΛCDM analysis. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

For the two stable parametrisations, we find the S₈ constraints to be consistent with the ΛCDM results with a Hellinger tension (see, e.g., Appendix G.1 of Heymans et al. 2021) of 0.54σ and 0.13σ, respectively. We find a slight reduction in constraining power with 1.2% and 1.1% precision measurements of S₈ compared to 1.0% precision under the assumption of ΛCDM and an increase in uncertainty on S₈-pagination by 18% and 6%, respectively. Parameterising modified gravity with the α functions, we again find the S₈ constraint to be consistent with the ΛCDM result at 1.06σ, although the uncertainty on S₈ increases by 31%, which corresponds to a measurement at 1.6% precision. We find a marginally better fit in modified gravity compared to ΛCDM as can be inferred from the best-fit χ² and the corresponding PTE, reported in the fourth and fifth column in Table 2. However, since the KiDS-Legacy data already provide an excellent fit to ΛCDM (Wright et al. 2025b), the model comparison between ΛCDM and modified gravity models only indicates a slight preference for modified gravity at the 1.4σ level, which we do not consider to be statistically significant.

We note that the weak lensing observable is sensitive to both the growth of matter fluctuations and the lensing potential. Therefore, a model that modifies μ and Σ can reproduce a given shear amplitude by compensating changes in the matter power spectrum, leading to a degeneracy between S₈ and the modified gravity functions. As shown in Sect. 4.3, the parametrisation in terms of the α functions allows for μ and Σ to vary largely independently, resulting in a weaker constraint on S₈. By contrast, in the stable parametrisation, our particular choice of time evolution for the basis functions imposes a correlation between μ and Σ, breaking the degeneracy with S₈. Therefore, the constraining power on S₈ for this particular parametrisation remains similar to that of ΛCDM.

4.2. Constraints on modified gravity parameters

Our parameter constraints on the Horndeski model parameters for the stable parametrisation in terms of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ under the assumption of a fixed sound speed are displayed in Fig. 2. Additionally, we display the derived posterior distribution of the braiding term at z = 0, α_B, 0. We conducted the analysis separately for the KiDS-Legacy, Planck, and DESI BAO + eBOSS RSD datasets, as indicated by the blue, green, and orange contours, respectively. The purple contour illustrates the posterior from a joint analysis. We find that both cosmic shear and CMB data yield competitive constraints on modified gravity parameters, while BAO and RSD data mostly reproduce the prior space. This is within our expectations, since BAOs only depend on the background quantities and therefore mainly provide additional constraints on the standard ΛCDM background parameters at late times. RSDs, on the other hand, probe the growth of structure and are thereby sensitive to modified gravity. However, RSD constraints are not competitive with those from cosmic shear and the CMB, yielding only a weak constraint on the braiding term. KiDS-Legacy and Planck provide a strong constraint on the effective Planck mass, for which we found a marginal mode and 68% HPDI of

Fig. 2. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ for KiDS-Legacy (blue), Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE (green), DESI DR2 BAO + eBOSS DR16 RSD (orange), and their combination (purple) in the parameter range constrained by the data. The list of top-hat priors for these parameters is provided in Table 1. Additionally, we show the posterior distribution of the derived value of the braiding term at z = 0, αB, 0. The black contour illustrates the prior volume. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. The dashed lines indicate the corresponding GR values.

$$\begin{array}{c}\hfill \mathrm{\Delta }{\widehat{M}}_{\ast }^2=0.{32}_{-0.21}^{+0.07}\end{array}$$(32)

Fig. 3. Samples of the posterior distribution of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ coloured by the value of the effective gravitational coupling, μ∞, in the small-scale limit, k → ∞, for the combination of KiDS-Legacy, Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE, DESI DR2 BAO, and eBOSS DR16 RSD datasets. The black contour illustrates the marginalised posterior distribution with the inner and outer contours corresponding to the 68% and 95% credible intervals, respectively. The dashed line shows a fit of the degeneracy for $\mathrm{\Delta }{\mu }_{\infty ,\mathrm{eff}}=\mathrm{\Delta }{\widehat{M}}_{\ast }^2{\widehat{D}}_{\mathrm{kin}}^{\gamma }$ with a best-fit value of γ = −1.220.

in the combined analysis of all probes. This excludes shifts towards lower values of the Planck mass, which is mainly driven by the cosmic shear constraint, and puts a strong upper limit on shifts towards higher values, but is still compatible with the ΛCDM assumption of $\mathrm{\Delta }{\widehat{M}}_{\ast }^2=0$ at 1.5σ. Additionally, we constrained the de-mixed kinetic term of the scalar field to be

$$\begin{array}{c}\hfill {\widehat{D}}_{\mathrm{kin}}=3.{74}_{-1.92}^{+0.69},\end{array}$$(33)

finding a 2σ preference for positive values of the braiding term, given by

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{B},0}=1.{44}_{-0.76}^{+0.39}.\end{array}$$(34)

We note that here, the upper boundary of the braiding term is not constrained by the data, but instead driven by the discontinuity at α_B = 2, as discussed in Sect. 2.1. While both Planck and KiDS data individually constrain modified gravity parameters, we find that the addition of KiDS data in the combined analysis leads to a reduction in the uncertainties of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ by 5% and 18%, respectively.

In the joint analysis of the KiDS-Legacy, Planck, DESI BAO, and eBOSS RSD datasets, we observe that the Horndeski model parameters are strongly degenerate, as can be inferred from the purple contour in Fig. 2. Here, we derive a correlation coefficient of r = 0.96 between ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, r = 0.97 between ${\widehat{D}}_{\mathrm{kin}}$ and α_B, 0, and r = 0.88 between $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and α_B, 0. This degeneracy is driven by the interplay between the model parameters, which can be inferred from the small scale limit of the effective gravitational coupling, μ_∞, given in Eq. (D4) in (CB24). At late times and for small ΔM_*², we can approximate μ_∞ by

$$\begin{array}{c}\hfill {\mu }_{\infty }=\frac{2D_{\mathrm{kin}}+{({\alpha }_{\mathrm{B}}+2\mathrm{\Delta }{M}_{\ast }^2)}^2}{2D_{\mathrm{kin}}(1+\mathrm{\Delta }{M}_{\ast }^2)},\end{array}$$(35)

where we set c_s² = 1 and approximated

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{M}}=3{\mathrm{\Omega }}_{\mathrm{m}}(a)\mathrm{\Delta }{M}_{\ast }^2(1+\mathrm{\Delta }{M}_{\ast }^2),\end{array}$$(36)

which follows from the specific functional form of ΔM_*². At late times and for a typical value of Ω_m ≈ 0.3, Eq. (36) yields α_M ≈ ΔM_*². From Eq. (35) we infer a degeneracy between D_kin, ΔM_*², and α_B. Here, α_B is determined integrating Eq. (27) in CB24, where the source term is controlled by D_kin and ΔM_*². We further investigate this degeneracy by in Fig. 3, which shows samples of the posterior coloured by the value of μ_∞. We find that for this particular parametrisation, Eq. (35) yields approximately constant values of μ_∞ along the degeneracy direction of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, which indicates a similar deviation from GR. Additionally, the values of Σ_L, ∞ are distributed similarly along the degeneracy direction, which is also reflected in the marginalised posterior distribution of μ_∞ and Σ_L,∞, discussed in Sect. 4.3. To describe the degeneracy between ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, we define

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mu }_{\infty ,\mathrm{eff}}=\mathrm{\Delta }{\widehat{M}}_{\ast }^2{\widehat{D}}_{\mathrm{kin}}^{\gamma },\end{array}$$(37)

which we fit to the posterior samples shown in Fig. 3, yielding

$$\begin{array}{c}\hfill \gamma =-1.220\pm 0.002,\end{array}$$(38)

with a marginal constraint on Δμ_eff of

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mu }_{\infty ,\mathrm{eff}}=0.066\pm 0.023.\end{array}$$(39)

Here, Δμ_∞, eff parametrises the deviation of μ_∞ from the GR value, which we find to be significant at approximately 2.9σ. However, this deviation from GR only leads to a small correction to the model predictions and a small improvement in the best-fit χ², as can be inferred from Table 2. Therefore, this model is found to not be statistically preferred over GR.

In Fig. 4, we show our constraints on the Horndeski parameter space when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ while fixing the effective Planck mass to the fiducial ΛCDM value, $\mathrm{\Delta }{\widehat{M}}_{\ast }^2=0$. Here, we observe that the constraints on modified gravity are almost entirely driven by the cosmic shear data, while the CMB dataset on its own does not yield significant constraints on the modified gravity parameters. However, as we show in Appendix B, the marginalisation over the lensing anomaly parameter A_L makes a significant impact on the inferred modified gravity parameter constraints for this particular parametrisation. We note that for this parametrisation the prior excludes negative values of α_B, 0.

Fig. 4. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ for KiDS-Legacy (blue), Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE (green), DESI DR2 BAO + eBOSS DR16 RSD (orange), and their combination (purple). Additionally, we show the posterior distribution of the braiding term at z = 0, αB, 0, and the combination ${\widehat{c}}_{\mathrm{sN}}^2={\widehat{D}}_{\mathrm{kin}}{\widehat{c}}_{\mathrm{s}}^2$, which is derived from the sampling parameters. The black contour illustrates the prior volume. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. In the GR limit, these quantities converge to zero, as the Horndeski scalar field does not exist within GR.

In this parametrisation, we mainly constrain the combination ${\widehat{c}}_{\mathrm{sN}}^2={\widehat{D}}_{\mathrm{kin}}{\widehat{c}}_{\mathrm{s}}^2$. This can be understood by inspecting the effective Newtonian coupling and the gravitational slip in the quasi-static approximation, given Eqs. (12) and (13). For ΔM_*² = 0, the small scale limits, μ_∞ and μ_Z, ∞, given in Appendix D in CB24, simplify to

$$\begin{array}{c}\hfill {\mu }_{\infty }={\mu }_{\mathrm{Z},\infty }=1+\frac{{\alpha }_{\mathrm{B}}^2}{2c_{\mathrm{sN}}^2},\end{array}$$(40)

yielding γ_QSA = 1. Thus, the modification to the lensing potential (i.e. the quantity our data are mostly sensitive to) becomes

$$\begin{array}{c}\hfill \mathrm{\Sigma }(k,a)=\frac{{\mu }_{\mathrm{p}}+{(k/aH)}^2{c}_{\mathrm{sN}}^2[1+{\alpha }_{\mathrm{B}}^2/(2c_{\mathrm{sN}}^2)]}{{\mu }_{\mathrm{p}}+{(k/aH)}^2{c}_{\mathrm{sN}}^2},\end{array}$$(41)

which only depends on c_sN² and the braiding term, α_B, that is determined through the inferred initial condition, α_B, 0. In practice, α_B is integrated from Eq. (27) in CB24, from which we can see that the source term for this model is primarily controlled by c_sN². This is reflected by the strong degeneracy between ${\widehat{c}}_{\mathrm{sN}}^2$ and α_B0 in Fig. 4. We additionally display the derived posterior distribution for ${\widehat{c}}_{\mathrm{sN}}^2$ in Fig. 4 and infer a marginal mode and HPDI of

$$\begin{array}{c}\hfill {\widehat{c}}_{\mathrm{sN}}^2=0.{34}_{-0.28}^{+0.28},\end{array}$$(42)

which is consistent with the ΛCDM expectation of ${\widehat{c}}_{\mathrm{sN}}^2=0$ at 1.2σ. Here, we find that the addition of KiDS-Legacy data in the combined analysis significantly reduces the uncertainty in ${\widehat{c}}_{\mathrm{sN}}^2$ by 67%. This highlights that, in this particular parametrisation, the constraints on modified gravity are largely driven by cosmic shear.

4.3. Constraints on derived modified gravity functions

As outlined in Sect. 2.1, the parametrisation of Horndeski gravity in terms of the α functions, defined in Eq. (6), is a common choice in the literature. In this section, we provide a comparison between the inferred α functions in the stable parameter basis and a direct sampling of the α functions. However, we note that the time evolution of the derived α differs from the time evolution of the parameters in Eq. (6), which are tied to the evolution of the dark energy density parameter. Therefore, we computed the evolution of α_B and α_M as a function of the scale factor and display the resulting constraints in Fig. 5 for the three modified gravity analysis setups. We note that when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, we assume a fixed value of the Planck mass ΔM_*² = 0, resulting in a fixed α_M = 0 in this particular setup. In Fig. C.1, we display samples of the derived α_B and α_M as a function of the scale factor, generated from the prior, in comparison to the posterior. The left panel of Fig. 6 shows a comparison between the marginalised 2D posterior distribution for the derived α_B and α_M at redshifts z = 0.0, z = 0.5, and z = 1.0 for the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$.

Fig. 5. Constraints on αB (left) and αM (right) as a function of scale factor for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour results from the analysis with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the orange contour corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The green contour shows the derived constraint on αB when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ for a fixed Planck mass and therefore imposing αM = 0. The solid line shows the mean of the posterior and the shaded region illustrates the 68% credible interval. The prior distributions of αB and αM for each model are displayed in Fig. C.1.

We find that the stable parameter basis in terms of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ prefers a stronger braiding term at late times compared to the common direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$. On the other hand, we find an overall lower run rate of the effective Planck mass, which we found to be close to the GR value, as indicated by the $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ constraint given in Eq. (32). Fixing ΔM_*² = 0, which implicitly assumes a zero run rate of the Planck mass, and sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ yields a α_B constraint that is compatible with the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$. Thus, we conclude that the inferred constraints on α_B and α_M are strongly dependent on the choice of model parameters in the stable parametrisation as already shown by, for instance, Noller & Nicola (2019) for the α-parametrisation. This is also reflected in the 2D posterior distribution, shown in the left panel of Fig. 6, which highlights different degeneracy directions between α_B and α_M at z = 0, which gradually align for increasing redshift. However, as can be inferred from Table 2, the data show no significant preference for a particular modified gravity model. We note that for the ${\widehat{\alpha }}_{\mathrm{B}}$-${\widehat{\alpha }}_{\mathrm{M}}$ parameter basis, the quantity (1 + z)α_i, depicted in Fig. 5, converges to a constant at low redshifts. This can be understood by expanding (1 + z)α_i for z → 0, yielding

Fig. 6. Posterior distribution of the derived αB and αM (left) and the effective Newtonian coupling μ∞ and the lensing modification ΣL, ∞ (right) at three distinct redshifts. Each posterior is derived from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contours are inferred with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the orange contours corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The black dashed lines indicate the GR limit set by αB = αM = 0 and μ∞ = ΣL, ∞ = 1. Each contour shows the 95% credible interval.

$$\begin{array}{c}\hfill (1+z)\frac{{\alpha }_i(z)}{{\alpha }_{i,0}}\approx 1+(1-3{\mathrm{\Omega }}_{\mathrm{m}})z,\end{array}$$(43)

which for a typical value of Ω_m ≈ 0.3 approaches an approximately constant value.

Additionally, we computed the effective Newtonian coupling and the lensing modification, as defined in Eqs. (15) and (17), for the three modified gravity parametrisations in the small-scale limit, k → ∞. The corresponding posterior distributions as a function of redshift are illustrated in Fig. 7, where the functional form is determined by the particular parametrisation of modified gravity. In Fig. C.2, we display samples of the derived μ_∞ and Σ_L, ∞ as a function of the scale factor, generated from the prior, in comparison to the posterior distributions. For the effective Newtonian coupling, we find a preference for μ > 1 at low redshift, indicating a deviation from GR, which predicts μ = 1. However, for the lensing modification, which is the quantity that is mostly constrained by the weak lensing data, we find that both the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and parametrisation via the α functions yield a posterior that is consistent with the GR value of Σ = 1. Therefore, the modified gravity models do not yield a significantly better fit to the data, resulting in the modified gravity model statistically not being preferred over ΛCDM, listed in Table 2. Additionally, we display the corresponding 2D posterior distributions for μ_∞ and Σ_L,∞ at three different redshifts in the right panel of Fig. 6. Here, we observe a degeneracy between both parameters in the stable basis, which is not present in the parametrisation in terms of the α functions. However, we find the constraints on the deviation from GR to be consistent between both parametrisations. We note that, as discussed in Sect. 4.2, we fix ΔM_*² = 0 when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ and therefore we impose μ = Σ for this particular parametrisation. Nevertheless, we do not find a significant deviation from ΛCDM for this parametrisation. We conclude that for the three parametrisations of modified gravity considered in this work, we find the data to be highly consistent with ΛCDM, preferring only small corrections to the ΛCDM prediction, as indicated by the constraints on the lensing modification, displayed in Figs. 6 and 7.

Fig. 7. Constraints on the effective Newtonian coupling μ∞ (left) and the lensing modification ΣL, ∞ (right) as a function of scale factor in the small-scale limit, k → ∞, inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour results from the analysis with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the green contour is inferred by sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, while the orange contour corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The dashed lines show the GR values of μ = ΣL = 1. The solid line shows the mean of the posterior and the shaded region illustrates the 68% credible interval.

4.4. Modified dark energy background evolution

In the previous sections, we have derived constraints on modified gravity under the assumption of a ΛCDM background. Here, we explore the possibility of further extending the background cosmological parameter space. In particular, we consider a dark energy density that is related to the pressure via a redshift-dependent equation of state, w(z), in the CPL parametrisation. This analysis is complementary to Reischke et al. (2025a), who derived constraints on dynamical dark energy from KiDS-Legacy in combination with additional low-redshift and CMB probes for a w₀w_aCDM cosmology. In Fig. 8, we show the marginalised posterior distribution for w₀ and w_a inferred from the combination of KiDS-Legacy, DESI, eBOSS, and Planck datasets for a modified gravity model parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ (blue contour) in comparison to a GR w₀w_aCDM analysis (orange contour).

Fig. 8. Constraints on w0 and wa for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour is inferred under the assumption of a modified gravity model parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ with a w0waCDM background cosmology. The black and orange contours show constraints for a w0waCDM cosmological model with fixed AL = 1 and AL as sampling parameter, respectively. The ΛCDM values are indicated by the dashed lines. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

We find the following constraints on the dark energy equation of state parameters:

$$\begin{array}{cc}\hfill w_0& =-0.{73}_{-0.13}^{+0.16},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill w_a& =-0.{75}_{-0.45}^{+0.36}.\hfill \end{array}$$(45)

We note that recent combined analyses of low-redshift and high-redshift probes have found evidence for dynamical dark energy (see for example Abdul Karim et al. 2025). Here, we find the equation of state parameter to be in better agreement with the ΛCDM value at the 2σ level. However, this shift can only be partly attributed to the modified gravity model, as can be inferred from a comparison of the blue and orange contours in Fig. 8. We observe a 0.5σ shift when comparing the GR w₀w_aCDM contour to the modified gravity contour with w₀w_aCDM background. Instead, the shift towards ΛCDM is mainly driven by the lensing anomaly parameter A_L, which we marginalised over in our analysis. For a fixed A_L = 1, we recover the preference for dynamical dark energy as indicated by the black contour, which is consistent with the findings of Abdul Karim et al. (2025) and Reischke et al. (2025a).

Figure 9 compares the marginalised posterior distribution of the modified gravity parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ for a w₀w_aCDM background (blue contour) to a ΛCDM background (green contour). We observe that opening up the parameter space of the background expansion does not make a significant impact on the inferred constraints on modified gravity. This retroactively confirms that for the cosmological probes and cosmological models considered in this work, imposing a ΛCDM background cosmology is a very good assumption. This is in agreement with Shah et al. (2026), who found that background and perturbations can be separated for most dark energy parametrisations. We constrain the structure growth parameter to be

$$\begin{array}{c}\hfill S_8=0.{834}_{-0.013}^{+0.015},\end{array}$$(46)

Fig. 9. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue and green contours assume a w0waCDM and ΛCDM background cosmology, respectively. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

which corresponds to a 1.76σ shift towards larger S₈ compared to the analysis with a ΛCDM background. A similar increase in S₈ was observed in Reischke et al. (2025a), who attributed this to an increase in Ω_m while σ₈ was found to be stable. Additionally, the w₀w_aCDM background yields a marginally better fit with Δχ² = −3.87 and is preferred over a ΛCDM background at 1.57σ in a model comparison via a Bayesian suspiciousness test. Thus, in agreement with Reischke et al. (2025a), we conclude that the preference for this extended background model is weak.

We note that since the modified gravity model considered in this work features a single field both affecting the growth of structure and driving the accelerated expansion of the Universe, the selection of a particular background imposes a mathematical stability prior on the modified gravity parameters. This is of particular importance in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, for which we found the ΛCDM background to only allow for positive values of the braiding term at all times. For a w₀w_aCDM background, on the other hand, we find the prior to be relaxed, allowing for a negative α_B. Furthermore, we find models with negative braiding to preferentially be located in the non-phantom region, defined by w_a > −1 − w₀, which, however, is disfavoured by the BAO and CMB data. Therefore, our parameter constraints under the assumption of a w₀w_a background cosmology are consistent with our fiducial constraints, inferred with a ΛCDM background. This conclusion holds for all three parametrisations of modified gravity considered in this work. However, by linear sampling in ${\widehat{D}}_{\mathrm{kin}}$ we exponentially down-weight large negative values for log₁₀α_B, 0. For a more exhaustive study of a dynamical dark energy background in future work, we therefore recommend a logarithmic sampling in ${\widehat{D}}_{\mathrm{kin}}$, which can in principle expose differences between the two backgrounds.

5. Conclusions

We present an analysis of the final cosmic shear dataset of the Kilo-Degree Survey (KiDS-Legacy) in a modified gravity framework. Given the consistency between KiDS-Legacy and other cosmological probes, we derived joint constraints with external data from BAOs, RSDs, and CMB anisotropies on the Horndeski class of modified gravity models, employing a parameter basis which satisfies stability criteria by construction. Furthermore, we investigated an extension to the assumed background cosmological model by adopting a dynamical dark energy model within the CPL parametrisation.

Our results demonstrate that cosmic shear can place significant constraints on the Horndeski parameter space, yielding limits that are competitive with, and can even surpass, those from the CMB. We found that modified gravity provides a slightly better fit to the combination of KiDS-Legacy, DESI DR2, eBOSS DR16, and Planck-Legacy data than GR for both ΛCDM and w₀w_aCDM cosmological models. However, Horndeski gravity is only preferred over ΛCDM at the 1.4σ level, indicating a high level of consistency of the data with the standard ΛCDM model. The inferred S₈ constraints are robust with respect to changes in the cosmological model, showing that cosmic shear data from KiDS-Legacy remains consistent with other probes when opening up the cosmological parameter space. Additionally, we observe that our constraints on modified gravity are mainly driven by the cosmic shear and CMB data, while BAO and RSD allow for the cosmological background model to be fixed at late times. Our results highlight that cosmic shear is a main driver of the modified gravity constraints. This is particularly evident in the parametrisation in terms of D_kin and c_s², where the addition of KiDS-Legacy data in the combined analysis reduces the uncertainty in the combination ${\widehat{c}}_{\mathrm{sN}}^2={\widehat{D}}_{\mathrm{kin}}{\widehat{c}}_{\mathrm{s}}^2$ by 67%.

At late times, our data prefer a modified gravity model in the stable basis of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, which deviates from GR at 2σ when considering the derived braiding term at the present time, α_B, 0. Computing the corresponding effective Newtonian coupling, we found a deviation from its GR value at 2.9σ. However, the overall correction to the GR prediction remains small, which is reflected in the low statistical preference over the GR model. Additionally, in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, the data mostly constrain the combination ${\widehat{c}}_{\mathrm{sN}}^2$, which we found to be consistent with GR at 1.2σ. Therefore, our best-fit Horndeski parameters are consistent with GR and indicate the level of deviation from GR that would remain compatible with the adopted data. Moreover, we compared our constraints on the braiding term, α_B, and the run rate of the effective Planck mass, α_M, to those inferred with a model featuring a direct sampling of the corresponding amplitudes, ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, with α_B(a) and α_M(a) proportional to the dark energy density, which is a commonly adopted parametrisation of Horndeski gravity in the literature. Here, we found the derived constraints on α_B(a) and α_M(a) to be strongly dependent on the choice of model parameters, which is consistent with earlier works on Horndeski gravity.

By allowing for an evolution of dark energy in our background cosmological model, we found our constraints on the modified gravity parameter space to be stable, which we attribute to the adopted data independently probing the evolution of the background and the growth of structure. When marginalising over the Planck lensing anomaly parameter, the data do not show a significant preference for a dynamical dark energy background over the standard ΛCDM background cosmology. We conclude that the final cosmic shear data from KiDS, together with a variety of other cosmological probes, are broadly consistent with structure formation governed by GR and CDM, while also remaining compatible with both Horndeski-class modified gravity and dynamical dark energy.

Acknowledgments

BS, MC, CH, and 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. RR and AD are supported by an ERC Consolidator Grant (No. 770935). BJ acknowledges support by the ERC-selected UKRI Frontier Research Grant EP/Y03015X/1 and by STFC Consolidated Grant ST/V000780/1. AL acknowledges support from the Swedish National Space Agency (Rymdstyrelsen) under Career Grant Project Dnr 2024-00171. AHW is supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR), under project 50QE2305, made possible by the Bundesministerium für Wirtschaft und Klimaschutz, and acknowledges funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source to Signal”. MA acknowledges support from 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 is supported by the Polish National Science Center through grants no. 2020/38/E/ST9/00395 and 2020/39/B/ST9/03494. CG is funded by the MICINN project PID2022-141079NB-C32. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. BG acknowledges support from the UKRI Stephen Hawking Fellowship (grant reference EP/Y017137/1). CH acknowledges support from the UK Science and Technology Facilities Council (STFC) under grant ST/V000594/1. HHi is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491, an ERC Consolidator Grant (No. 770935), and the DLR project 50QE2305. 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. SSL has received funding from the programme “Netzwerke 2021”, an initiative of the Ministry of Culture and Science of the State of Northrhine Westphalia. LL is supported by the Austrian Science Fund (FWF) [ESP 357-N]. CM acknowledges support from the Beecroft Trust, the Spanish Ministry of Science under the grant number PID2021-128338NB-I00, and from the European Research Council under grant number 770935. LM 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”. LP acknowledges support from the DLR grant 50QE2302. TT acknowledges funding from the Swiss National Science Foundation under the Ambizione project PZ00P2_193352. MvWK acknowledges the support by UK STFC (grant no. ST/X001075/1), the UK Space Agency (grant no. ST/X001997/1), and InnovateUK (grant no. TS/Y014693/1). MY acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program with Grant agreement No. 101053992. YZ acknowledges the studentship from the UK Science and Technology Facilities Council (STFC). Kilo-Degree Survey: 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. Dark Energy Spectroscopic Instrument: This research used data obtained with the Dark Energy Spectroscopic Instrument (DESI). DESI construction and operations is managed by the Lawrence Berkeley National Laboratory. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: www.desi.lbl.gov/collaborating-institutions. The DESI collaboration is honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies. SDSS-IV: Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. Planck: Based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Software: The figures in this work were created with MATPLOTLIB (Hunter 2007) and CHAINCONSUMER (Hinton 2016), making use of the NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (McKinney 2010), COSMOSIS (Zuntz et al. 2015), NAUTILUS (Lange 2023), and MOCHI_CLASS (Cataneo & Bellini 2024) software packages. Author Contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (BS,RR,MG,MC,BJ,AL,ASM), 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 since its inception.

References

Appendix A: Full cosmological parameter constraints

In this appendix, we provide additional constraints on the full set of parameters for the three parametrisations of modified gravity in comparison to a ΛCDM analysis. In Fig. A.1, we display the marginalised posterior distribution for the cosmological and nuisance sampling parameters inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. We note that we do not show the redshift and IA nuisance parameters, since their posterior distributions are dominated by the Gaussian prior listed in Table 1. Table A.1 lists the marginal mode and the 68% HPDI of all sampled and derived parameters. Additionally, we list the marginal modified gravity parameter constraints inferred from the individual datasets, corresponding to the marginalised posterior distributions presented in Figs. 2 and 4, in Table A.2.

Fig. A.1. Marginalised posterior distributions of the cosmological and nuisance sampling parameters inferred from the combination of KiDS-Legacy with DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The black contour shows constraints from a ΛCDM analysis. The blue, green, and orange contours illustrate the resulting posterior from modified gravity analyses using the stable parameter basis of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the stable parameter basis of ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, and the αB-αM parametrisation, respectively. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. The list of top-hat priors for these parameters is provided in Table 1.

Appendix B: Impact of the CMB lensing anomaly on modified gravity constraints

As discussed in Sect. 2.4, we marginalised over the lensing anomaly parameter A_L when analysing Planck data in our modified gravity model. In this appendix, we test the impact of fixing the value of A_L to unity on the inferred constraints on modified gravity. We conducted the analysis for the combination of all datasets considered in this work and display the resulting constraints on S₈ and Ω_m, the best-fit χ², the PTE, the difference in χ² between ΛCDM and modified gravity models, and the N_σ preference level for modified gravity in Table B.1. Comparing these results to Table 2, we find the cosmological parameter constraints to be consistent with our fiducial constraints, although we observe a tendency of S₈ shifting towards larger values by less than 1σ. Moreover, we find that fixing A_L overall provides a worse fit to the data in all models. These results are consistent with Planck Collaboration VI (2020), who reported that marginalising over A_L provides a better fit to Planck data under the assumption of a ΛCDM model and found a similar variation in the cosmological parameters. When changing from a ΛCDM model to a modified gravity model, we find that the improvement in best-fit χ² is larger for a fixed A_L, indicating that the modified gravity model can, at least to some extent, correct for the apparent systematics in the Planck data. This increases the preference for the modified gravity model, as can be inferred from the rightmost column in Table B.1. However, with values of up to N_σ = 1.84 we do not find a statistically significant preference for the modified gravity model in this analysis setup.

Fig. B.1. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ (left) and ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ (right) for a fixed AL = 1. Additionally, we show the posterior distribution of the derived value of the braiding term at z = 0, αB, 0. The blue contour shows constraints from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets, while the green contour illustrates constraints from an analysis with Planck data only. For comparison, the black contour shows our fiducial constraints with free AL for the combination of all probes. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

In Fig. B.1, we show the constraints on the modified parameters inferred for the two stable parametrisations with A_L = 1 from the combination of all probes (blue contours) and from Planck only. For comparison, we display the fiducial constraints for a free A_L for the combination of all probes (black contours). In the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, we find the Planck contours to prefer a similar region in parameter space as the ones shown in Fig. 2. Thus, the lensing anomaly parameter does not significantly alter the modified gravity parameter constraints for this particular parametrisation. In the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, however, we observe a significant shift in the contour inferred with Planck data. This can be seen most prominently in the posterior distribution of the derived parameter ${\widehat{c}}_{\mathrm{sN}}^2$, which directly impacts the modification to the lensing potential. Therefore, fixing A_L = 1 can directly be corrected by increasing ${\widehat{c}}_{\mathrm{sN}}^2$ because of the degeneracy between the two parameters. Additionally, the Planck data on their own reveal a strong preference for a non-zero braiding term in this case. When combining the CMB and cosmic shear data, however, we added more data constraining the lensing potential, which tend to prefer lower values of the lensing potential. This results in the peak of the ${\widehat{c}}_{\mathrm{sN}}^2$ posterior shifting towards zero. Nevertheless, we find that not accounting for the lensing anomaly can lead to a biased constraint on modified gravity parameters, and therefore we conclude that marginalising over A_L plays an important role in the analysis presented in this work.

Appendix C: Prior distribution of derived modified gravity functions

In this appendix, we present the prior distribution of the derived braiding parameter, α_B, the derived run rate of the effective Planck mass, α_M, the effective Newtonian coupling μ_∞ and the lensing modification Σ_L, ∞. For each parametrisation of modified gravity considered in this work, we generated 500 samples from the prior, listed in Table 1. For each sample, we computed the corresponding α_B(a), α_M(a), μ_∞(a) and Σ_L, ∞(a) using MOCHI_CLASS. In Fig. C.1 we present the prior samples of α_B(a) and α_M(a) in comparison to the posterior distributions inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets, shown in Fig. 5. Here, the three rows show the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, and the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, respectively. We note that when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, we assume a fixed value of the Planck mass $\mathrm{\Delta }{\widehat{M}}_{\ast }^2=0$, resulting in a fixed α_M = 0, and therefore we only infer the prior distribution of α_B(a) for this particular setup. In Fig. C.2 we present the prior samples of μ_∞(a) and Σ_L, ∞(a) in comparison to the posterior distributions inferred from the same combination of datasets. We note that here, the fixed value of the Planck mass in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ imposes μ_∞ = Σ_L, ∞.

Fig. C.1. Prior distribution of the derived α functions generated by drawing 500 samples from the prior, given in Table 1. Each line illustrates the evolution of the derived αB and αM as a function of scale factor for one set of parameters generated from the prior. Additionally the coloured lines and shaded regions show the corresponding posterior distribution inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The three rows show the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, and the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, respectively.

Fig. C.2. Prior distribution of the effective Newtonian coupling, μ∞, and the lensing modification, ΣL, ∞, generated by drawing 500 samples from the prior, given in Table 1. Each line illustrates the evolution of μ∞ and ΣL, ∞ as a function of scale factor for one set of parameters generated from the prior. Additionally the coloured lines and shaded regions show the corresponding posterior distribution inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The three rows show the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, and the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, respectively.

All Tables

All Figures

Fig. 1. Posterior distribution in the S8-Ωm plane for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contours illustrate the posterior inferred in a Horndeski model parametrised by the stable basis parameters ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, and the green contour shows constraints from sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$. The orange contour shows the posterior in a Horndeski model using the ${\widehat{\alpha }}_{\mathrm{B}}$-${\widehat{\alpha }}_{\mathrm{M}}$ parametrisation. The black contour corresponds to the fiducial ΛCDM analysis. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

In the text

Fig. 2. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ for KiDS-Legacy (blue), Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE (green), DESI DR2 BAO + eBOSS DR16 RSD (orange), and their combination (purple) in the parameter range constrained by the data. The list of top-hat priors for these parameters is provided in Table 1. Additionally, we show the posterior distribution of the derived value of the braiding term at z = 0, αB, 0. The black contour illustrates the prior volume. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. The dashed lines indicate the corresponding GR values.

In the text

Fig. 3. Samples of the posterior distribution of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ coloured by the value of the effective gravitational coupling, μ∞, in the small-scale limit, k → ∞, for the combination of KiDS-Legacy, Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE, DESI DR2 BAO, and eBOSS DR16 RSD datasets. The black contour illustrates the marginalised posterior distribution with the inner and outer contours corresponding to the 68% and 95% credible intervals, respectively. The dashed line shows a fit of the degeneracy for $\mathrm{\Delta }{\mu }_{\infty ,\mathrm{eff}}=\mathrm{\Delta }{\widehat{M}}_{\ast }^2{\widehat{D}}_{\mathrm{kin}}^{\gamma }$ with a best-fit value of γ = −1.220.

In the text

Fig. 4. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ for KiDS-Legacy (blue), Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE (green), DESI DR2 BAO + eBOSS DR16 RSD (orange), and their combination (purple). Additionally, we show the posterior distribution of the braiding term at z = 0, αB, 0, and the combination ${\widehat{c}}_{\mathrm{sN}}^2={\widehat{D}}_{\mathrm{kin}}{\widehat{c}}_{\mathrm{s}}^2$, which is derived from the sampling parameters. The black contour illustrates the prior volume. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. In the GR limit, these quantities converge to zero, as the Horndeski scalar field does not exist within GR.

In the text

Fig. 5. Constraints on αB (left) and αM (right) as a function of scale factor for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour results from the analysis with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the orange contour corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The green contour shows the derived constraint on αB when sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ for a fixed Planck mass and therefore imposing αM = 0. The solid line shows the mean of the posterior and the shaded region illustrates the 68% credible interval. The prior distributions of αB and αM for each model are displayed in Fig. C.1.

In the text

Fig. 6. Posterior distribution of the derived αB and αM (left) and the effective Newtonian coupling μ∞ and the lensing modification ΣL, ∞ (right) at three distinct redshifts. Each posterior is derived from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contours are inferred with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the orange contours corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The black dashed lines indicate the GR limit set by αB = αM = 0 and μ∞ = ΣL, ∞ = 1. Each contour shows the 95% credible interval.

In the text

Fig. 7. Constraints on the effective Newtonian coupling μ∞ (left) and the lensing modification ΣL, ∞ (right) as a function of scale factor in the small-scale limit, k → ∞, inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour results from the analysis with Horndeski gravity modelled by the stable parameter basis ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ and the green contour is inferred by sampling ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, while the orange contour corresponds to a direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$ with αB(a) and αM(a) proportional to the dark energy density, following Eq. (6). The dashed lines show the GR values of μ = ΣL = 1. The solid line shows the mean of the posterior and the shaded region illustrates the 68% credible interval.

In the text

Fig. 8. Constraints on w0 and wa for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue contour is inferred under the assumption of a modified gravity model parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ with a w0waCDM background cosmology. The black and orange contours show constraints for a w0waCDM cosmological model with fixed AL = 1 and AL as sampling parameter, respectively. The ΛCDM values are indicated by the dashed lines. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

In the text

Fig. 9. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ for the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The blue and green contours assume a w0waCDM and ΛCDM background cosmology, respectively. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

In the text

Fig. A.1. Marginalised posterior distributions of the cosmological and nuisance sampling parameters inferred from the combination of KiDS-Legacy with DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The black contour shows constraints from a ΛCDM analysis. The blue, green, and orange contours illustrate the resulting posterior from modified gravity analyses using the stable parameter basis of ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the stable parameter basis of ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, and the αB-αM parametrisation, respectively. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively. The list of top-hat priors for these parameters is provided in Table 1.

In the text

Fig. B.1. Posterior distribution of Horndeski parameters in the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$ (left) and ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$ (right) for a fixed AL = 1. Additionally, we show the posterior distribution of the derived value of the braiding term at z = 0, αB, 0. The blue contour shows constraints from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets, while the green contour illustrates constraints from an analysis with Planck data only. For comparison, the black contour shows our fiducial constraints with free AL for the combination of all probes. The inner and outer contours of the marginalised posteriors correspond to the 68% and 95% credible intervals, respectively.

In the text

Fig. C.1. Prior distribution of the derived α functions generated by drawing 500 samples from the prior, given in Table 1. Each line illustrates the evolution of the derived αB and αM as a function of scale factor for one set of parameters generated from the prior. Additionally the coloured lines and shaded regions show the corresponding posterior distribution inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The three rows show the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, and the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, respectively.

In the text

Fig. C.2. Prior distribution of the effective Newtonian coupling, μ∞, and the lensing modification, ΣL, ∞, generated by drawing 500 samples from the prior, given in Table 1. Each line illustrates the evolution of μ∞ and ΣL, ∞ as a function of scale factor for one set of parameters generated from the prior. Additionally the coloured lines and shaded regions show the corresponding posterior distribution inferred from the combination of KiDS-Legacy, DESI DR2 BAO, eBOSS DR16 RSD, and Planck 2018 TTTEEE, low-ℓ TT, and low-ℓ EE datasets. The three rows show the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and $\mathrm{\Delta }{\widehat{M}}_{\ast }^2$, the direct sampling of ${\widehat{\alpha }}_{\mathrm{B}}$ and ${\widehat{\alpha }}_{\mathrm{M}}$, and the stable basis parametrised by ${\widehat{D}}_{\mathrm{kin}}$ and ${\widehat{c}}_{\mathrm{s}}^2$, respectively.

In the text