Main

One of the most enduring puzzles in modern cosmology is the accelerating expansion of the Universe. The impetus for introducing another component beyond ordinary matter first arose in the 1980s when early cosmic microwave background (CMB) anisotropy calculations showed that a low-density open model would conflict with small-scale CMB observations. Subsequently, large-scale clustering analyses1 and baryon-fraction arguments2 made it clear that a Universe consisting solely of critical-density dark matter was untenable, prompting the emergence of a cosmological constant (Λ) as the leading explanation. Observations of type Ia supernovae then provided direct evidence for an accelerating cosmic expansion3,4, and by the turn of the millennium, the Λ cold dark matter (ΛCDM) framework was firmly established as the consensus model of cosmology. This paradigm has been continuously supported by increasingly precise CMB measurements5 and baryon acoustic oscillation (BAO) data6, solidifying its status as the standard picture of cosmic acceleration.

Despite this consensus, the fundamental cause of the cosmic acceleration remains elusive. A widely favoured interpretation involves postulating an effective dark energy (DE) component that dominates the current energy density of the Universe (see refs. 7,8,9 for reviews). Provided its speed of sound is not too small, DE can be treated as nearly homogeneous on sub-sound-horizon scales, leaving its main cosmological signatures determined by its energy density ρ and pressure P. The equation of state, w ≡ P/c2ρ, serves as a crucial diagnostic: w = −1 corresponds to the vacuum energy, whereas dynamical DE models allow w to vary with time, as exemplified by quintessence10, phantom11, quintom12, k essence13, Chaplygin gas14 and holographic DE15. Subclasses of these scenarios can be further differentiated by examining their trajectories in the phase space spanned by (w, dw/dln a), where a is the cosmic scale factor16.

Recent observational progress, including supernova measurements from the Dark Energy Survey17, as well as other large supernovae compilations, such as Pantheon+18,19 and Union20, high-precision CMB data from the Planck Collaboration21, and BAO results from the Sloan Digital Sky Survey eBOSS programme22, has largely refined our ideas about DE. These complementary probes constrain DE-related parameters with growing precision, thus enhancing our understanding of cosmic acceleration. Although combined analyses that integrate all available datasets provide the most stringent overall limits, it remains instructive to consider each probe in isolation. Such a modular approach not only underscores the unique strengths and systematics of each measurement but also furnishes critical cross-checks, reinforcing the robustness of DE constraints.

Data Release 1 (DR1)23 and Data Release 2 (DR2)24,25 of BAO observations by the Dark Energy Spectroscopic Instrument (DESI) place stringent limits on the expansion history of the background of the Universe, achieving a level of precision on a par with recent supernovae surveys. These include the Pantheon+ compilation of 1,550 spectroscopically confirmed supernovae18,19, the Union3 compilation of 2,087 supernovae20 and the 5-year Dark Energy Survey sample (DESY5)17.

In this paper we probe DE dynamics along two complementary routes. First, we fit the distance ladder with a shape-function formalism that makes no explicit assumptions about the equation of state. Second, we carry out a fully Bayesian reconstruction of w(z) in which a correlation prior derived from the Horndeski gravity26 yields an interpretable Bayes factor, recovering the standard (w0, wa) result27,28 with two effective degrees of freedom. Here, w0 is the present-day equation-of-state value, and wa describes how w evolves with redshift. Although these techniques differ from the Gaussian-process and (w0, wa) analyses presented in the companion DESI studies25,29, all approaches paint a coherent picture: DESI DR1 and DR2 BAO measurements and three different type Ia supernova samples indicate the same mild, oscillatory departure from a cosmological constant. The convergence of results across datasets and methodologies strengthens the statistical case for an evolving DE component.

A shape-function analysis of DE

The shape function S0(a) illustrates how DE grows relative to matter, and the function S1(a) reveals whether and how the DE pressure changes with time. Finally, S2(a)—often referred to as the ‘state-finder’ parameter30,31—encodes both w(a) and its derivative, making it particularly useful for distinguishing among DE models (see Methods for the definition and more details of the shape functions). For example, a ‘freezing’ quintessence scenario predicts S2 < −1 (ref. 16), whereas barotropic fluid models require S2 > −1 (ref. 32). If S2 crosses the −1 threshold during cosmic evolution, that may indicate that several DE components dominate at different epochs.

Figure 1 displays the reconstructed S functions derived from DESI BAO and three supernovae datasets, based on equations (2)–(5). All panels reveal a systematic deviation from the ΛCDM baseline, and this pattern of deviation is consistent across the four datasets, indicating a shared underlying trend.

Fig. 1: Reconstructed shape functions Si (i = 0, 1 or 2) derived from several cosmic-distance measurements, displayed as a function of redshift z.
figure 1

All five distance probes show a consistent, ΛCDM-deviating trend in the S functions, strongest in DESI DR2. Each column corresponds to a different dataset. From left to right: DESI BAO DR1, DESI BAO DR2, Pantheon+, Union3 and DESY5. From top to bottom, the S0(z), S1(z) and S2(z), as defined in equations (5), which capture the time evolution of the DE energy density, the DE pressure and the DE equation of state, respectively. The solid white lines indicate the best-fitting reconstruction, and the blue-shaded region denotes the 68% confidence-level uncertainty, derived according to equations (2)–(5). The green dashed-dotted curve represents the best-fitting shape function from the CPL parameterization, and the horizontal dashed line is the corresponding prediction of the ΛCDM model. See the text for further details.

Full size image

A key point is that the DESI DR2 sample uncovers a stronger dynamical DE signal than DESI DR1. In particular, S0 and S1 both lie systematically above the ΛCDM expectation in DR2, pointing to a more pronounced departure from a pure cosmological constant. Even more notably, S2 in DR2 exhibits a sharper crossing behaviour relative to DR1, indicating that the observed expansion may not be fully described by standard ΛCDM.

This amplified crossing, coupled with smaller uncertainties, bolsters the case for DE having extra internal degrees of freedom. Reinforcing this conclusion, the DR2-based reconstructions closely match those obtained from three independent supernovae samples. Their broad consistency indicates that these deviations probably reflect the real dynamical evolution of the DE; for this not to be the case, a number of different systematics would all need to conspire in the same direction. Although DR2 provides stronger statistical evidence overall, uncertainties remain large at z ≳ 0.5, limiting definitive conclusions at higher redshifts. Even so, the alignment between the DESI DR2 and supernovae data underscores the importance of forthcoming high-precision measurements for clarifying the fundamental nature and potential evolution of DE.

Evolution of the equation of state

To explore this further, we adopt the commonly used w0–wa parameterization27,28, in which the DE equation of state is written as w(a) = w0 + wa(1 − a). Here, w0 is the present-day equation-of-state value, and wa describes how w evolves with redshift. Notably, ΛCDM is recovered by setting w0 = −1 and wa = 0. Fitting this model to the data constrains the possible dynamical nature of DE in a simple yet effective manner. Moreover, the w0–wa approach serves as a bridge between the shape-function reconstructions and the fully non-parametric w(z) methods discussed later, offering an intermediate step towards more general DE modelling.

Figure 2 shows the inferred constraints on (w0, wa) for different data combinations (see Methods for the datasets and priors used). Figure 2a corresponds to DESI DR1, and the Fig. 2b presents DR2. Contour colours indicate various datasets, with inner and outer contours marking the 1σ and 2σ confidence intervals, respectively. We divide the parameter space into four DE models: quintessence (w > −1 at all epochs, thus w0 > −1 and w0 + wa > −1)10, full phantom (w < −1 at all epochs, thus w0 < −1 and w0 + wa < −1)11, quintom A (w > −1 in the past but w < −1 today, thus w0 < −1 and w0 + wa > −1)12 and quintom B (w < −1 in the past but w > −1 today, thus w0 > −1 and w0 + wa < −1)12.

Fig. 2: Various data combinations with 68% and 95% credible-interval constraints on the DE equation-of-state parameters (w0, wa).
figure 2

a, Results based on DESI DR1 measurements. b, Results based on DESI DR2 measurements. Different colours indicate distinct dataset combinations, as labelled. In all cases, we assume in addition a BBN constraint on Ωbh2 and a CMB constraint on θ*. The parameter space is subdivided according to four DE models: quintom A (w > −1 in the past and w < −1 today), quintessence (w > −1 at all epochs), phantom (w < −1 at all epochs) and quintom B (w < −1 in the past but w > −1 today). The intersection point of the two dashed lines corresponds to the ΛCDM limit (w0, wa) = (−1, 0).

Full size image

Using only DESI DR1 BAO measurements yields constraints that are broadly consistent with ΛCDM, albeit with mild (~1.5–2σ) hints of departure. Incorporating supernova data narrows these contours further, indicating a slight preference for dynamical DE. A similar pattern emerges from the DR2-only analysis, which achieves notably tighter constraints than DR1 and favours w0 > −1 with wa < 0, consistent with a quintom B-like scenario. Even so, ΛCDM remains viable at ~1.5σ. This is weaker than the 2.3σ rejection of ΛCDM in ref. 25 from DESI + CMB because our CMB prior uses only part of the information and also because we have chosen to inflate the uncertainties. When DR2 is combined with supernova data, the best-fitting values stay in the quintom B region, pushing the tension with ΛCDM above 2σ, showing somewhat stronger preference for departure from a simple cosmological constant (see Supplementary Tables 1 and 2 for more details).

To delve deeper into the dynamical nature of DE, we perform a non-parametric Bayesian reconstruction of w(z) using the correlation-prior framework originally proposed in refs. 33,34 and subsequently applied in refs. 35,36. In this methodology, w(z) is treated as a free function and a correlation prior enforces smoothness and mitigates flat directions in the likelihood. The correlation prior itself can be motivated by theoretical considerations37; here, we adopt the version derived in ref. 37 from Horndeski theory26 (see Methods for details of the reconstruction).

Figure 3 (white curves and blue bands) illustrates the best fit and 68% confidence-level constraints on w(z) obtained from various data combinations. The upper row corresponds to DESI DR1 BAO data, and the lower row to DR2. In each row, the leftmost panel reports BAO-only constraints, whereas the other panels incorporate other supernova datasets: Pantheon+, Union3 and DESY5. In all cases, w(z) ≠ −1 is clearly indicated. Throughout the paper, we define the significance (that is, the signal-to-noise ratio, SNR) of the reconstructed w(z) deviating from a specific model wmod:

$${{\rm{SNR}}}^{2}\equiv {({\bf{w}}-{{\bf{w}}}_{\mathrm{mod}})}^\mathrm{T}\,{{{C}}}_{{\rm{w}}}^{-1}\,({\bf{w}}-{{\bf{w}}}_{\mathrm{mod}}),$$
(1)

where Cw is the posterior covariance for the w bins. For the DR1-based reconstructions (Fig. 3a–d), the SNRs of w ≠ −1 are 2.6, 3.9, 3.9 and 4.2, respectively. Switching to the DR2 dataset yields comparable significance at 2.6, 3.7, 4.3 and 4.5, respectively. It is interesting that the preference for deviations from ΛCDM when the supernovae samples are included remains strong with the enhanced statistical power of DESI DR2.

Fig. 3: DE equation of state w(z) reconstructed from several datasets.
figure 3

a–d, DESI DR1 BAO results, both alone (a) and in combination with three independent type Ia supernova datasets: DR1 + Pantheon+ (b), DR1 + Union3 (c) and DR1 + DESY5 (d). The results are for two different approaches: the correlation-prior method (bottom-layered, dark-blue band) and the (w0, wa) parameterization (top-layered, green band). The solid and dashed white lines indicate the best-fitting reconstructions of w(z), and the shaded regions are 68% confidence-level uncertainties. The horizontal dashed line denotes the ΛCDM prediction (w = −1). e–h, Results for the DESI DR2 BAO sample. e, DR2. f, DR2 + Pantheon+. g, DR2 + Union3. h, DR2 + DESY5. The BBN prior and constraints on the CMB acoustic scale θ* are imposed.

Full size image

All dataset combinations exhibit a persistent pattern in w(z): w > −1 at z ≲ 0.2 and w < −1 at z ≈ 0.75, which is consistent with the results reported in companion DESI papers25,29. Including supernova data highlights this feature and reveals pronounced oscillations in w(z). These patterns remain stable and gain further significance when DR1 BAO is replaced by DR2. In fact, oscillatory or non-monotonic w(z) evolution arises naturally in several classes of DE models, including axion-like quintessence scenarios with periodic potentials38, interacting dark sector frameworks, where energy exchange between DE and dark matter can induce transient oscillations in w (refs. 39,40), and multi-field quintom constructions that allow w(z) to cross the −1 threshold and back12,41. Crucially, the oscillatory feature seen in our reconstruction is not an artefact of any single dataset: it appears consistently across different type Ia supernova samples and remains essentially unchanged when the BAO dataset is upgraded from DESI DR1 to DR2. This stability across independent datasets and updated observations reinforces the interpretation that the observed oscillations in w(z) are a robust physical signal, potentially pointing towards an underlying dynamical DE component consistent with the aforementioned model classes.

The green bands in Fig. 3, which are reconstructions using the Chavallier–Polarski–Linder (CPL) parameterization, track the transition across w = −1 well and are broadly compatible with the non-parametric results, reinforcing the robustness of these conclusions.

Although the significance of w ≠ −1 stems from the reduction in χ2 (χ2 = −ln L, where L is the likelihood function) relative to ΛCDM, fully assessing the fit also requires determining the effective degrees of freedom in the reconstructed w(z). This step is not trivial because the correlation prior couples the w(z) bins. To address this, we perform a principal component analysis (PCA)42 to decorrelate the prior and posterior inverse covariance matrices of w(z). Further technical details of the PCA procedure can be found in Supplementary Information.

Figure 4 shows the PCA results for the reconstructed DE equation of state w(z). Figure 4a–d plots the eigenvectors of the first few well-constrained modes, highlighting their redshift dependence. Figure 4e–h displays the corresponding eigenvalues λi for different dataset combinations. PCA identifies independent modes constrained by the data. In Fig. 4a (BAO only), the first principal component (PC1) is smooth and relatively featureless, indicating that BAO data alone primarily constrain a single dominant mode of w(z). By contrast, adding supernova data (Fig. 4b–d) produces further constrained modes; PC2 and PC3 become more oscillatory, especially at low redshift, demonstrating that supernova observations notably improve the constraints on the redshift evolution of w(z). Because the main feature of w(z) is effectively a linear combination of these three modes, further progress will come from shrinking the PC1 error bar through improved low-z supernova calibration, tightening PC2 and PC3 with denser distance tracers in the 0.3 < z < 0.8 window and extending observations to z ≳ 0.8 to assess whether other PCs are required.

Fig. 4: PCA of the DE equation of state w(z).
figure 4

a–d, Plots of the first three principal components ei(z) as functions of redshift. a, BAO. b, BAO + Pantheon+. c, BAO + Union3. d, BAO + DESY5. The different dataset combinations constrain distinct modes in w(z). e–h, The corresponding eigenvalues λi versus the mode index. Solid lines distinguish DESI DR1 BAO analyses, and dashed lines indicate DESI DR2. e, BAO. f, BAO + Pantheon+. g, BAO + Union3. h, BAO + DESY5. Revealed are the relative information that each dataset provides compared with the correlation-prior baseline (green curves). The BBN prior and constraints on the CMB acoustic scale θ* are imposed in all cases. Larger eigenvalues signify better-constrained modes.

Full size image

Figure 4e–h shows the eigenvalues λi, which quantify the statistical significance of each mode (larger values correspond to tighter constraints). In Fig. 4e (also BAO only), four eigenvalue notably exceeds that of the correlation prior (green curve), which means that four eigenmodes well survive the correlation prior. When supernova datasets are incorporated (Fig. 4f–h), more modes become well constrained. Comparing DESI DR1 (orange) and DR2 (blue) eigenvalues reveals that DR2 consistently achieves larger eigenvalues for the constrained modes, indicating stronger statistical power in the DR2 sample. This gain is most evident when supernova data are included, as the gap between the first few eigenvalues of DR1 and DR2 widens, showing that DR2 improves constraints across several modes of w(z). We further validate these findings by computing the effective number of degrees of freedom Neff numerically with covariant PCA43, as detailed in the Supplementary Information. The covariant PCA results align with the PCA results, confirming the robustness of our assessment of how many independent modes of w(z) are constrained by the data.

We can project w(z) onto the PC eigenmodes 1 + w(z) = ∑iγiei(z) where the corresponding eigenvalue λi of the covariance matrix represents the error of the amplitude γi, λi = 1/σ2(γi). Supplementary Table 5 lists the amplitudes of the first five PCs. For DR2 and for DR1 or DR2 plus supernovae, the amplitude of the first PC is positive, indicating that this mode describes a behaviour of w(z) that increases from −1 over time. The improvement in the fits, Δχ2 = −∑i[γi/σ(γi)]2, is dominated by PC1. On including higher-order PCs, several crossings of w = −1 occur, as we observe in the reconstructed w(z) shown in Fig. 3. For DR2 plus supernovae, PC4 gives the second largest contribution to Δχ2.

An advantage of the Bayesian reconstruction is that the dependence of the reconstructed w(z) on the correlation prior is explicit. To demonstrate this, we deliberately change the amplitude of the correlation prior A and impose a stronger prior to avoid possible overfitting. By tuning A, the effective number of degrees of freedom Neff can be reduced. When Neff = 2 the reconstructed w(z) is consistent with the results that we obtained based on the w0–wa parameterization, as shown in Supplementary Fig. 1, with related Δχ2 shown in Supplementary Table 6.

Bayesian evidence and evolution

We evaluate the evidence against ΛCDM analytically and numerically (see Methods for details), as illustrated in Fig. 5. Consistent with expectations, ln L grows while ln V (where V is volume) decreases monotonically with increasing Δ (or Neff), which reflects the number of modes in w(z) that remain after applying the correlation prior. The analytic and numerical evidence estimates typically agree, except in BAO-only scenarios, where the data are not sufficiently constraining to maintain Gaussianity in the posterior. Figure 5 plots the model evidence and related quantities versus Δ, with the effective degrees of freedom Neff noted on the top axis. The top panels show DESI DR1 results, and the bottom panels DESI DR2. We display the change in log-likelihood Δ ln L (blue squares), the analytic Δ ln E (where E is evidence, orange circles), the numerical Δ ln E (red points) and the change in log-prior volume Δ ln V (green diamonds).

Fig. 5: Likelihood, prior volume and evidence as functions of the effective degrees of freedom Neff.
figure 5

Bayesian evidence peaks at ‘moderate’ support for dynamical w(z) once supernovae are included. The strength of this prior is varied by adding different values of Δ to its diagonal. For comparison, we show both the analytically estimated ln E (calculated as ln L + ln V) and the numerically computed values obtained with PolyChord. Top: results derived using DESI DR1 BAO data. Bottom: results with DESI DR2. The BBN priors and constraints on the CMB acoustic scale θ* are imposed in all cases. The shaded regions highlight values of Δ for which the correlation prior is affected by the flat prior of the w bins. See the text for further discussion.

Full size image

The increase in both Δ ln E and the SNR for detecting w ≠ −1 with larger Neff further supports the idea that adding extra degrees of freedom favours a dynamical DE interpretation over a pure cosmological constant. When using BAO data alone, Δ ln E remains consistent with zero for all Neff, confirming that BAO measurements by themselves offer limited evidence against ΛCDM. By contrast, once supernova datasets are added, Δ ln E rises substantially. For instance, DR1 + DESY5 yields a peak Δ ln E = 4.1 ± 0.6 (classified as moderate on Jeffreys’ scale44,45,46) at Neff = 2, increasing to 5.2 ± 0.7 (also moderate) for DR2 + DESY5 at Neff = 3. A similar pattern emerges with BAO + Union3: the evidence systematically increases from DR1 to DR2, reaching 3.3 ± 0.7 (moderate) at Neff = 3 for DR2.

By contrast, BAO + Pantheon+ shows weaker evidence, changing only slightly from DR1 to DR2, with Δ ln E peaking at 1.5 ± 0.6 (weak) in DR1 and 1.4 ± 0.7 (weak) in DR2, consistent with what is seen under the CPL parameterization. Although the Pantheon+ dataset alone can show a mild preference for w ≠ −1 (Supplementary Fig. 2), combining it with BAO data partially cancels that signal. One possibility is a mild tension between the Pantheon+-inferred distance scale and the BAO-preferred expansion history, pushing the joint constraints closer to ΛCDM. Further investigation is needed to determine whether this is attributable to systematics.

To gauge the strongest possible deviation from w = −1 without overfitting the data, we look at w(z) CDM models with the largest departure from ΛCDM while still having a positive evidence, that is Δ ln E ≥ 0. Under this condition, the significance of w ≠ −1 for BAO + DESY5 reaches 4.1 (DR1) and 4.3 (DR2). For BAO + Union3, the corresponding values are 3.8 (DR1) and 3.9 (DR2). Moreover, BAO + Pantheon+ yields maximum SNRs of 3.2 (DR1) and 3.1 (DR2), reinforcing the previous finding that Pantheon+ offers comparatively weaker support for w ≠ −1 than either Union3 or DESY5.

Summary and conclusions

With the start of the DESI survey, we now have access to exceptionally precise BAO measurements, which are a robust tool for probing the nature of DE. In this work, we analysed DESI DR2 BAO observations in tandem with complementary cosmic-distance indicators—supernovae data and the CMB acoustic scale θ*—to test for deviations from the ΛCDM paradigm. By employing a shape-function approach47,48, we derived key observables sensitive to DE and found consistent evidence for an evolving equation of state w(z) across several independent datasets.

Building on these findings, we performed a non-parametric Bayesian reconstruction of w(z) by jointly analysing DESI BAO, supernovae and CMB data, employing a correlation prior rooted in Horndeski theory. Our results exhibit strong statistical support for dynamical DE, at a confidence level exceeding 3σ. Rigorous model selection further indicates that DESI BAO, supernovae and CMB observations together constrain roughly three independent degrees of freedom without overfitting. In particular, combining DESI BAO, CMB and the DESY5 supernova data detects DE evolution at SNR = 3.9, with a Bayesian evidence of 5.2 ± 0.7 in favour of w(z) CDM over ΛCDM when Neff = 3.

These results indicate that DE may not be properly described by a simple cosmological constant, as they favour a time-dependent equation of state w(z) instead. Although earlier work36 hinted at similar conclusions, this study provides some of the tightest constraints so far, thanks in large part to the enhanced statistical power of DESI DR2 BAO data. Furthermore, adding supernova datasets substantially boosts the sensitivity to potential DE evolution: all three independent supernova samples (Pantheon+, Union3 and DESY5) show broadly consistent indications of deviations from ΛCDM.

Companion DESI papers25,29 investigate the time evolution of DE using complementary methods and also report the detection of the dynamical DE. The level of significance is generally greater in those works because the full Planck measurements are used there, as opposed to the θ* information used in this work. Moreover, recent results from the Atacama Cosmology Telescope indicate that the time evolution of DE is still favoured when the Planck data are replaced with the measurements made by the Atacama Cosmology Telescope49.

Extensive cross-checks indicate that known systematics cannot mimic the oscillatory pattern in Fig. 3. Photometric offsets and a possible evolution of the colour versus luminosity relation in the Pantheon+ data could shift w(z) by about 0.03 at z ≲ 0.4 (ref. 19), and the 0.04 mag low- to high-z mismatch removes only a smooth trend in (w0, wa) fits when using the DESY5 supernovae sample50. For the BAO analysis, DESI DR1 analyses bound template- and reconstruction-induced BAO shifts to 0.1% (ref. 51). All these effects vary smoothly with redshift: they can tilt or raise the whole w(z) curve but cannot produce the alternating-sign wiggles that we recovered. The feature, therefore, remains most naturally interpreted as evidence for genuine DE dynamics.

In conclusion, our study finds tentative evidence for departures from ΛCDM using the DESI DR2 data. Although our results favour a dynamical DE interpretation, continued observational and theoretical efforts will be critical in establishing whether these deviations signify new physics or stem from unknown systematics. The complete DESI sample, with increased statistical depth and extended redshift coverage, is expected to yield even stricter constraints on DE evolution. Other surveys, such as Euclid52 and the prime focus spectrograph at the Subaru Telescope53, will deliver complementary BAO and large-scale structure measurements, enabling cross-validations of DESI results. Likewise, next-generation CMB projects will improve our knowledge of early-Universe physics, strengthening the overall DE picture.

Methods

Datasets and priors

DESI at Kitt Peak National Observatory in Arizona marks a great step forwards in stage IV large-scale structure surveys54,55. Designed with a 3.2° diameter prime focus corrector56, DESI deploys 5,000 fibres57 through a robotic focal plane assembly58 to simultaneously measure several spectra. These hardware capabilities are bolstered by a high-throughput spectroscopic data-processing pipeline59 and an efficient operations plan60. Since beginning operations in 2021, DESI has been gathering high-fidelity spectra for numerous target classes: the bright galaxy sample (0.1 < z < 0.4)61, luminous red galaxies (0.4 < z < 1.1)62, emission-line galaxies (1.1 < z < 1.6)63, quasars (0.8 < z < 2.1)64 and the Lyman α forest (1.77 < z < 4.16)65. Early data analyses and a survey validation55,66 have confirmed that DESI satisfies the performance benchmarks required for a stage IV survey and is advancing our grasp of the role of DE in the cosmic expansion54.

DESI DR1 spans observations from 14 May 2021 to 14 June 202223 and has already yielded new insights into DE by detecting the BAO signature in galaxy and quasar clustering51 as well as in the Lyman α forest67. These results fit seamlessly with external data68, paving the way for more comprehensive analyses that incorporate the full clustering information from the DESI galaxy sample and other tracers69,70. The subsequent DR2, covering observations through 9 April 202424,25, not only encompasses the entirety of DR1 but also extends the redshift range and statistical reach of the survey. This expanded dataset offers a critical test of dynamical DE: if DE truly evolves over cosmic time, the enriched DR2 sample could reveal a more pronounced signature. By comparing DR1 and DR2 results, one can assess both the internal consistency and the potential evidence for evolving DE. This paper contributes to a suite of companion analyses that build on the main cosmological results25, including a related supporting study on DE29 and neutrino constraints71.

Our analysis draws upon three categories of measurements: BAO observations from DESI DR1 and DR2 (refs. 51,67), supernovae samples from Pantheon+, Union3 and DESY5, and distance constraints from the CMB. Because DE chiefly affects the expansion history of the background of the Universe and modelling its perturbations requires extra assumptions, we focus here on background observables. Specifically, we incorporated the DESI DR1 and DR2 BAO measurements, supernovae luminosity distances and the Big Bang nucleosynthesis (BBN) prior on the physical baryon density, Ωbh2 = 0.02196 ± 0.00063 (ref. 72). We also included CMB distance information from Planck PR4 (ref. 73), which imposes 100θ* = 1.04098 ± 0.00042, where θ* ≡ r*/DM(z*) is the ratio of the comoving sound horizon at recombination r* to the transverse comoving distance DM(z*) (ref. 21,68). The quantity r* was computed using RECFAST74 within CAMB75. As in ref. 68, we inflated the error bar from ref. 73 (column PR4_12.6 TTTEEE in their Table 5) by 75% to account for the potential broadening of θ* constraints in extended cosmological models.

Our baseline dataset comprises DESI BAO measurements, the BBN prior on Ωbh2 and the CMB acoustic-scale constraint on θ*. We then supplemented this baseline with each of the three supernovae samples in turn: DR1(2) + Pantheon+, DR1(2) + Union3 and DR1(2) +;DESY5 for the DESY5 sample. Each combination was analysed independently to gauge its impact on w(z) constraints.

The shape functions of DE

Galaxy surveys measure cosmic distances using (DM/rd, DH/rd) or DV/rd at specific redshifts, where DM is the comoving distance, DH ≡ c/H with H being the Hubble expansion rate and rd is the sound horizon scale. The volume-averaged distance is defined by \({D}_{{\rm{V}}}\equiv {\left(z{D}_{{\rm{M}}}^{2}{D}_{{\rm{H}}}\right)}^{1/3}\), with z being the redshift. These observations constrain the shape of H(z), which is directly related to DE. To link DH and DM, we adopted the parameterization introduced in ref. 76, which provides an accurate yet general description of cosmic distances across a wide range of cosmologies:

$$\begin{aligned}\frac{{D}_{{\rm{M}}}(z)}{{D}_{{\rm{M}}}^{{\rm{f}}}(z)}\,R&={\alpha }_{0}\left(1+{\alpha }_{1}x+{\alpha }_{2}{x}^{2}+{\alpha }_{3}{x}^{3}+\ldots \right),\\ \frac{{D}_{{\rm{H}}}(z)}{{D}_{{\rm{H}}}^{{\rm{f}}}(z)}\,R&={\beta }_{0}\left(1+{\beta }_{1}x+{\beta }_{2}{x}^{2}+{\beta }_{3}{x}^{3}+\ldots \right),\end{aligned}$$
(2)

where the superscript f denotes quantities evaluated in a fiducial ΛCDM model with the current fractional matter density ΩM = 0.3153. The normalization factor \(R\equiv {r}_{{\rm{d}}}^{{\rm{f}}}/{r}_{{\rm{d}}}\) accounts for discrepancies in the sound horizon scale, and \(x\equiv {D}_{{\rm{M}}}^{{\rm{f}}}(z)/{D}_{{\rm{M}}}^{{\rm{f}}}({z}_{{\rm{p}}})-1\), where DM denotes the comoving distance and zp is a pivot redshift for the series expansion.

Because DH and DM are interdependent, the coefficients βi follow from αi:

$${\beta }_{0}={\alpha }_{0}\left(1+{\alpha }_{1}\right),\quad {\beta }_{i}=(i+1)\,\frac{{\alpha }_{i}+{\alpha }_{i+1}}{1+{\alpha }_{1}}\quad (i > 0).$$
(3)

To avoid overfitting, the series expansion was truncated once the uncertainty in any αi exceeded unity. For DESI DR1 and DR2 BAO data, a cutoff at i = 3 (yielding four α parameters) proved adequate.

Given the values of α, one can construct the following function47,48:

$$S(a)\equiv{}A{H}^{2}(a){a}^{3}=BX(a){a}^{3}+C,$$
(4)

where A, B and C are constants, and X(a) ≡ ρDE(a)/ρDE(1) is the normalized DE density. As S(a) has the same functional form as ρDE(a)—up to a shift and overall rescaling—this ‘shape function’ effectively encodes the potential dynamical properties of DE.

Specifically, we defined three dimensionless functions that offer clear diagnostics for DE evolution:

$$\begin{aligned}{S}_{0}(a)&\equiv {a}^{3}-\frac{3\left[S(a)-S(1)\right]}{{S}^{{\prime} }(1)}={a}^{3}+\frac{X(a){a}^{3}-1}{w(1)}\mathop{\to}\limits^{\varLambda }1,\\ {S}_{1}(a)&\equiv \frac{1}{{a}^{3}}\,\frac{{S}^{{\prime} }(a)}{{S}^{{\prime} }(1)}=\frac{{P}_{{\rm{DE}}}(a)}{{P}_{{\rm{DE}}}(1)}\mathop{\to}\limits^{\varLambda }1,\\ {S}_{2}(a)&\equiv -\frac{{S}^{{\prime\prime} }(a)}{3{S}^{{\prime} }(a)}=w(a)-\frac{{w}^{{\prime} }(a)}{3w(a)}\mathop{\to}\limits^{\varLambda }-1,\end{aligned}$$
(5)

where derivatives are taken with respect to ln a. Here, PDE(a) and w(a) denote the pressure and equation of state of DE, respectively, and the arrows indicate their expected values under ΛCDM.

A non-parametric Bayesian reconstruction of DE

Within the w(z) CDM framework, w(z) is discretized into N = 29 piecewise-constant segments, w ≡ {w1, w2, …, wN}, uniformly spaced by a scale factor for z < 2.5. Additionally, an extra ‘wide’ bin wwide = −1 was assigned for z > 2.5 up to recombination. Allowing wwide to vary left it largely unconstrained by our data, and it showed minimal correlations with the other w bins. We, therefore, fixed wwide = −1 for simplicity. Alongside these wi, we varied several cosmological parameters: ΩM, the current fractional matter density; Ωbh2, the physical baryon density; H0, the Hubble constant; and \({{\mathcal{M}}}_{{\rm{b}}}\), the absolute supernovae magnitude (used only when supernova data were included).

We placed sufficiently broad flat priors on all parameters. Specifically, H0 ∈ [20, 100], ΩM ∈ [0.01, 0.99], w0 ∈ [−3, 1], wa ∈ [−3, 2] and wi ∈ [−6, 4] (for i = 1, …, 29). We explored the resulting parameter space using the Markov chain Monte Carlo (MCMC) algorithm implemented in Cobaya77. At each step, we computed a total χ2 that combines contributions from both the data and the correlation prior. The prior covariance CΠ was derived from Horndeski-based simulations37 spanning a broad parameter space, as investigated in ref. 37. Specifically, the prior is defined as

$$\begin{array}{lll}C(a,{a}^{\prime})&\equiv& \langle [w(a)-{w}_{{\rm{fid}}}(a)][w(a)-{w}_{{\rm{fid}}}({a}^{\prime})]\rangle \\ &=&\sqrt{C(a)C({a}^{\prime})}R(a,{a}^{\prime}),\end{array}$$
(6)

where C(a) ≡ C(a, a) and R(a, a′) is the normalized correlation function, which thus, equals unity for a = a′. The functional form of the correlation prior used in this work is taken from37:

$$\begin{aligned}C(a)&=0.05+0.8{a}^{2},\\ R(a,{a}^{{\prime} })&=\exp \left[-{\left(\left\vert \ln a-\ln {a}^{{\prime} }\right\vert /0.3\right)}^{1.2}\right].\end{aligned}$$
(7)

The total χ2 minimized in the MCMC process has two pieces:

$$\begin{aligned}{\chi }^{2}&={\chi }_{{\rm{data}}}^{2}+A{\chi }_{{\rm{prior}}}^{2},\\ {\chi }_{{\rm{prior}}}^{2}&={{\bf{D}}}_{{\bf{w}}}^\mathrm{T}\,{{{C}}}_{\Pi }^{-1}\,{{\bf{D}}}_{{\bf{w}}}={{\bf{w}}}^\mathrm{T}\,{\tilde{{{C}}}}_{\Pi }^{-1}\,{\bf{w}},\end{aligned}$$
(8)

where Dw = w − wfid = (I − S)w and \({\tilde{{{C}}}}_{\Pi }^{-1}\), the inverse modified covariance matrix of the prior, is defined as \({\tilde{{{C}}}}_{\Pi }^{-1}\equiv {\left({{I}}-{{S}}\right)}^\mathrm{T}\,{{{C}}}_{\Pi }^{-1}\left({{I}}-{{S}}\right)\). Matrices I and S are the identity matrix and the transformation matrix, respectively, and the fiducial model wfid = Sw was calculated by taking a local average of five neighbouring w bins through the transformation matrix S (refs. 34,35,36,78). The factor A balances the strength between data and prior, and we set A = 1 as a default to obtain the main result in Fig. 3.

Before applying it to actual observations, we validated our pipeline using mock datasets generated for four input DE models with the same data covariance matrix for DESI and supernova measurements. Further technical details are provided in Supplementary Information.

Calculating the Bayesian evidence

For a rigorous test of whether the w(z) CDM model is preferred over ΛCDM, we evaluated the Bayesian evidence using both PolyChord79,80 within Cobaya77 and an analytic approximation. Under the assumption that both prior and posterior distributions are Gaussian, the evidence E can be computed as35 ln E = ln V + ln L, where \(\ln V=(1/2)\left(\operatorname{ln}\operatorname{det}{{{C}}}_{{\rm{post}}}-\operatorname{ln}\operatorname{det}{\tilde{{{C}}}}_{{\rm{prior}}}\right)\). Here, \({\tilde{{{C}}}}_{{\rm{prior}}}\) and Cpost are the modified correlation-prior covariance and the posterior covariance, respectively, and ln L is the maximum log-likelihood from the data and prior. Because our reconstruction uses a running-average method to specify the fiducial w(z), \({\tilde{{{C}}}}_{{\rm{prior}}}\) has zero eigenvalues that render its volume ill-defined34,35 (Supplementary Information).

To address this, we added a diagonal term Δ−2 to the prior correlation matrix, interpolating between ΛCDM (Δ = 0) and a w(z) CDM model in which wfid is unconstrained by the correlation prior. The prior volume is controlled by Δ and, for Δ < 1, it became independent of the flat prior imposed on w in the MCMC analysis. The parameter Δ changed the effective number of degrees of freedom Neff by controlling the strength of the correlation prior.