Extended Dark Energy analysis using DESI DR2 BAO measurements

(DESI Collaboration) {sup}29`Institut d’Estudis Espacials de Catalunya (IEEC), c/ Esteve Terradas 1, Edifici RDIT, Campus PMT-UPC, 08860 Castelldefels, Spain

{sup}30`Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain

{sup}31`The Ohio State University, Columbus, 43210 OH, USA

{sup}32`Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

{sup}33`Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

{sup}34`Instituto de Física, Universidad Nacional Autónoma de México, Circuito de la Investigación Científica, Ciudad Universitaria, Cd. de México C. P. 04510, México

{sup}35`CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

{sup}36`NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA

{sup}37`Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

{sup}38`Department of Physics & Astronomy and Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara Street, Pittsburgh, PA 15260, USA

{sup}39`University of Chinese Academy of Sciences, Nanjing 211135, People’s Republic of China.

{sup}40`University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, USA

{sup}41`Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia

{sup}42`Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia

{sup}43`Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

{sup}44`Scientific Computing Core, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

{sup}45`Departament de Física Quàntica i Astrofísica, Universitat de Barcelona, Martí i Franquès 1, E08028 Barcelona, Spain

{sup}46`Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), c. Martí i Franquès, 1, 08028 Barcelona, Spain.

{sup}47`Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA

{sup}48`Institut d’Astrophysique de Paris. 98 bis boulevard Arago. 75014 Paris, France

{sup}49`Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

{sup}50`Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

{sup}51`School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

{sup}52`Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA

{sup}53`Department of Physics and Astronomy, University of California, Irvine, 92697, USA

{sup}54`Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France

{sup}55`Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, 1156 High Street, Santa Cruz, CA 95064, USA

{sup}56`Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95065, USA

{sup}57`Department de Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain

{sup}58`Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA

{sup}59`Laboratoire de Physique Subatomique et de Cosmologie, 53 Avenue des Martyrs, 38000 Grenoble, France

{sup}60`Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

{sup}61`Department of Physics and Astronomy, Siena College, 515 Loudon Road, Loudonville, NY 12211, USA

{sup}62`Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

{sup}63`Department of Physics & Astronomy, University of Wyoming, 1000 E. University, Dept. 3905, Laramie, WY 82071, USA

{sup}64`Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada

{sup}65`Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada

{sup}66`Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada

{sup}67`Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA

{sup}68`Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

{sup}69`Max Planck Institute for Extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany

{sup}70`Department de Física, EEBE, Universitat Politècnica de Catalunya, c/Eduard Maristany 10, 08930 Barcelona, Spain

{sup}71`Université Clermont-Auvergne, CNRS, LPCA, 63000 Clermont-Ferrand, France

{sup}72`Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Republic of Korea

{sup}73`Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia

{sup}74`Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, USA

{sup}75`Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia

{sup}76`Department of Physics & Astronomy, Ohio University, 139 University Terrace, Athens, OH 45701, USA

{sup}77`Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany

{sup}78`University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81677 München, Germany

Introduction

The Λ Cold Dark Matter (ΛCDM) model has withstood the test of time as the standard framework of modern cosmology, and it provides a robust foundation for understanding the Universe. It describes a spatially flat universe that is homogeneous and isotropic on large scales, governed by Einstein’s general relativity [1]. The model incorporates two key components: about 70% is in dark energy that is described by the vacuum-energy contribution (corresponding to the cosmological constant Λ in the equations), while another 30% is in pressureless matter that is made up of a combination of cold dark matter (CDM) and baryons. Despite its elegant simplicity, ΛCDM has successfully explained a broad range of cosmological observations [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. On the whole, measurements made over the past several decades have largely confirmed this paradigm and, in particular, cemented dark energy [15, 16, 17] as the essential component of concordance model to explain the observed accelerated expansion of the Universe [2, 3].

While the cosmological constant has been a cornerstone of the standard model of cosmology, various dark energy models with an evolving equation of state have been proposed as alternatives [18, 19, 20, 21, 22, 23, 24, 25, 26]. Specifically, we are motivated to study these time-evolving alternatives by the recent cosmological results from the Dark Energy Spectroscopic Instrument (DESI) [27, 28]. DESI is able to measure multiple spectra simultaneously by means of its 5,000 fibers [29] and a robotic plane assembly [30] across the field of view given its 3.2° diameter prime focus corrector [31]. This is complemented by a high-performance spectroscopic data processing pipeline [32] and a streamlined operations plan [33]. DESI is designed to help better understand the nature of dark energy [34] and its successful survey validation [35] based on early data [36] showed that it meets the expected requirements of a Stage-IV survey. In particular, its Data Release 1 (DR1) [37] has already provided new insights into the behavior of dark energy. DESI DR1 measured the baryon acoustic oscillations (BAO) signature in the clustering of galaxies and quasars [38], as well as the Lyman-$\alpha$ forest [39]. The combined constraints from DESI DR1 BAO and external data [40], followed up with a similar analysis that combines the BAO with the full clustering information from DESI galaxies and other tracers [41, 42], as well as the supporting DESI DR1 papers that considered alternative descriptions of the dark energy sector [43, 44], all showed tantalizing hints of the departures from the cosmological constant dark energy model. Cosmological hints in the dark energy sector are currently a source of debate, and it is of high priority to explore them with more data. In this work, we make use of the BAO measurements from the second data release (DR2 [45, 46, 47]) from DESI to explore the possibility of an evolving, dynamical dark energy, and evaluate whether existing observations support such a paradigm shift. This paper is part of a set of supporting papers that aim to extend the cosmological analysis presented in [47] (see [48] for the supporting paper focusing on neutrino constraints).

An essential ingredient, in a study confronting dark-energy models with data, is the physical description of dark energy (DE). In the standard concordance model, $\Lambda$CDM, it is described by its contribution to the stress-energy tensor, $\Lambda$ or, equivalently, by its energy density relative to critical, $\Omega_{\mathrm{DE}}$. A dynamical dark energy is enabled by allowing the equation of state, $w = P/(\rho c^2)$, to differ from its $\Lambda$CDM value of $-1$. There are many possible ways to achieve this, a large number of which have been introduced and tested in the literature [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62]. We can classify them as parametric and non-parametric approaches. Parametric approaches rely on predefined functional forms for quantities like $w(z)$ (where $z$ is redshift), while non-parametric methods seek to reconstruct these quantities directly from data without assuming predefined functional forms or specific cosmological models. Both methods have advantages and disadvantages. On the one hand, parametric models are mathematically simple and easy to interpret, but may lead to biased inferences if the assumed parametric form deviates substantially from reality. On the other hand, non-parametric methods offer greater flexibility and are less subject to model-dependent biases. However, these are harder to implement and require careful validation with simulations. For this reason, we perform initial tests using simulated (mock) datasets. While there is no substitute for comprehensive validation, these tests check the methodology’s implementation and mitigate potential biases that could affect the results. We remind readers that all the analyses in this paper rely on the assumption that the data used are reliable and free from unknown systematics.

The paper is structured as follows. In Section II, we introduce the datasets and general methodology used in the analyses, followed by Section III, where we summarize the current status of the DESI results from the $w_0w_a$ parametrization [47]. Various alternative dark energy parametrizations are explored in Section IV, before the implementation of non-parametric methods in Section V. Section VI provides an interpretation of the possible physical mechanisms behind deviations from $\Lambda$CDM. Finally, in Section VII, we present our conclusions.

Datasets and Methodology

In this section, we provide a brief cosmological background on distance measurements relevant to DESI, with an emphasis on dark energy. We start by introducing the relevant cosmological functions, before proceeding to describe the datasets used and the statistical tools employed in our analysis.

As shown in [40], the evidence for spatial curvature in the Universe is not significant. Therefore, we assume a flat universe for all the results presented in our analyses. The time-dependence of the dark energy density is enabled via the equation of state $w(z)$; the expansion rate reads

where $\Omega_{\mathrm{bc}} = \Omega_c + \Omega_b$, $H_0$ is the Hubble parameter today, and $\Omega_b$, $\Omega_c$, $\Omega_\gamma$, $\Omega_\nu$, and $\Omega_{\mathrm{DE}}$ are the present-time energy density parameters in baryons, cold dark matter, radiation, massive neutrinos and dark energy, respectively. The neutrino species contribute to the matter content of the Universe at the present day, since they behave as non-relativistic matter once they have passed through a transition redshift during the matter domination era (e.g., transitioning around a redshift $\sim100$ for a neutrino mass eigenstate with a mass of $0.06$ eV) [63, 64]. This detail will be important when defining our cosmic microwave background (CMB) compression scheme, since relativistic neutrinos do not contribute significantly to matter density at the time of recombination. We define $\Omega_m = \Omega_{\mathrm{bc}} + \Omega_\nu$ to denote the matter content that scales $(1+z)^3$ when neutrinos are non-relativistic.

For a dark energy component with an equation-of-state parameter $w(z) = P_{\mathrm{DE}}(z)/(\rho_{\mathrm{DE}}(z)c^2)$, the energy density $\rho_{\mathrm{DE}}$ normalized to its present value evolves as

For a constant value of $w(z)$, the dark energy density becomes proportional to $(1+z)^{3(1+w)}$, while for a model based a cosmological constant ($w=-1$), the right-hand side of (2) is unity. The conventional $w_0w_a$ parametrization for time-varying $w(z)$ is [50, 51]

with energy density following the expression:

BAO measure the comoving distance at the effective redshift of a given galaxy sample, in units of the sound horizon ($r_s$) at the drag epoch, labeled as $r_d \equiv r_s(z_d)$. The drag epoch corresponds to the release of baryons from the drag of CMB photons, occurring at a redshift ($z_d$). The scale $r_d$ is thus the distance that sound waves in the photon-baryon fluid were able to travel all the way from the big bang, slightly after the time of recombination, to the drag epoch, given by

where $c_s(z)$ is the speed of sound waves in the fluid, and $z_d \approx1060$ is the redshift at which photons and baryons decouple [7]. The BAO measurements are sensitive to the distance in the direction transverse to the line-of-sight, corresponding to the comoving angular diameter distance

BAO also measure the comoving distance along the line-of-sight, which is directly related to the expansion rate

However, as described in [47], some DESI BAO measurements are isotropic, as in the case of the BGS tracer. Hence, we also make use of the spherically-averaged distance $D_V(z)$ that quantifies the average of the distances measured along, and perpendicular to, the line-of-sight to the observer [65], and is given by

Since these measurements are relative to the sound horizon $r_d$, which sets the BAO scale, the directly constrained quantities are the ratios $D_M/r_d$, $D_H/r_d$, and $D_V/r_d$. With this, we can now define the primary dataset used for our searches of dynamical dark energy, based on the latest DESI data:

• Baryon acoustic oscillations (BAO): We use the BAO distance measurements from DESI DR2, as detailed in Table III in [47]. In particular, for the BGS tracer, we use measurements of $D_V/r_d$ providing compressed low redshift information from the range $0.1<z<0.4$. For the rest of DESI tracers, we use the BAO distance measurements of $D_M/r_d$ and $D_H/r_d$. Explicitly, we use two LRG bins in the ranges $0.4<z<0.6$ and $0.6<z<0.8$, a combined tracer measurement for LRG+ELG in the range $0.8<z<1.1$, a measurement spanning $1.1<z<1.6$ for the ELG tracer and the QSO in the range $0.8<z<2.1$. The systematics tests associated with the BAO measurements from galaxy and quasar clustering are presented in [66]. We also include the Ly$\alpha$ measurements in $1.8<z<4.2$, which provides our highest redshift data-point. This measurement is described in detail in [46] (see also [67] for validation tests and [68] for specific catalog details). We refer to this whole dataset, encompassing information from redshift $0.1$ to $4.2$, split into seven main samples, as “DESI”.

We now proceed to define the cosmological datasets that we use, in combination with DESI, to obtain constraints on cosmological parameters. The cosmological probes and specific external datasets used in our analysis are:

• Supernovae Ia (SNe Ia): We combine DESI data with either of the following three SNe Ia datasets, namely PantheonPlus, Union3, and DESY5. The PantheonPlus [69] dataset comprises 1550 spectroscopically-confirmed SNe Ia in the redshift range $0.001<z<2.26$. The Union3 compilation [70] has 2087 SNe Ia in the redshift range $0.01<z<2.26$, 1363 of which are common to PantheonPlus, though the analysis methodologies are

substantially different. Finally, the DESY5 dataset [71] is a sample of 1635 photometrically-classified SNe Ia with redshifts in the range $0.1 < z < 1.13$, complemented by 194 historical low-redshift SNe Ia (which are also present in the PantheonPlus sample) spanning $0.025 < z < 0.1$.

• Cosmic microwave background (CMB): We include temperature and polarization measurements of the CMB from the Planck satellite [72]. In particular, we use the high-$\ell$ TTTEEE likelihood (planck_NPIPE_highl_CamSpec.TTTEEE), together with low-$\ell$ TT (planck_2018_lowl.TT) and low-$\ell$ EE (planck_2018_lowl.EE) [73, 74], as implemented in Cobaya [75]. Additionally, we combine temperature and polarization anisotropies with CMB lensing measurements from the combination of NPIPE PR4 from Planck [76, 77] and the Atacama Cosmology Telescope (ACT) DR6 [78, 79].

• Compressed CMB: We use the Gaussian correlated prior over $\omega_b \equiv\Omega_b h^2$, $\omega_{bc} \equiv\Omega_b h^2 + \Omega_c h^2$ and $\theta_*$ as defined in [47]. Here, the angular acoustic scale $\theta_*$ adds extra geometrical information from the CMB, while $\omega_b$ and $\omega_{bc}$² serve to set the sound horizon $r_d$ and calibrate our BAO measurements. These CMB-based quantities capture most of the relevant information from the early CMB by marginalizing over contributions from late-time effects, such as the integrated Sachs–Wolfe (ISW) effect and CMB lensing, resulting in a robust CMB compression for testing late-time physics [80]. In particular, we use these compressed measurements as a conservative alternative for constraining dark energy at the background level, thereby allowing for negative $f_{\rm DE}(z)$, as in Sections IV B and V A. For brevity, we refer to these as $(\theta_*,\omega_b,\omega_{bc})_{\rm CMB}$.

In our analysis, we utilize Markov Chain Monte Carlo (MCMC) sampling to explore the parameter space using the Metropolis-Hastings algorithm [81, 82] as implemented in Cobaya [83, 84]. For the alternate parametrizations, non-parametric methods, and DE classes, we adopt priors similar to [40], with exact specifications presented in Table I, and have modified the Boltzmann solver camb [85, 86], incorporating a generalized equation of state for dark energy for the theoretical prediction of observables. We employ the parametrized post-Friedmann (PPF) framework [87, 88] to compute cosmological perturbations for the time-dependent equation of state $w(a)$, where $a$ is the scale factor, which permits transitions across the phantom divide at $w=-1$. Additionally, we use custom theory code in Cobaya for the analysis of $f_{\rm DE}(z)$ binning and crossing statistics. For quintessence

parametrization	parameter	default	prior
Baseline	$\omega_{\mathrm{cdm}}$	—	$\mathcal{U}[0.001,0.99]$
	$\omega_b$	—	$\mathcal{U}[0.005,0.1]$
	$100\theta_{\mathrm{MC}}$	—	$\mathcal{U}[0.5,10]$
	$\ln(10^{10}A_s)$	—	$\mathcal{U}[1.61,3.91]$
	$n_s$	—	$\mathcal{U}[0.8,1.2]$
	$\tau$	—	$\mathcal{U}[0.01,0.8]$
in absence of $\theta_{\mathrm{MC}}$	$H_0$	—	$\mathcal{U}[20,100]$
Alt. Parametrization	$w_0$	$-1$	$\mathcal{U}[-3,1]$
	$w_a$	$0$	$\mathcal{U}[-3,2]$
Crossing	$C_0$	$1$	$\mathcal{N}[1,1^2]$
	$C_i$	$0$	$\mathcal{N}[0,1^2]$
Binning	$w_i$	$-1$	$\mathcal{U}[-3,1]$
	$f_{\mathrm{DE},i}$	$1$	$\mathcal{U}[-5,5]$
Gaussian Processes	$\ell_f$	—	Eq. (C3)
	$w(z_l)$	$-1$	$\mathcal{N}[-1,1^2]$
Dark Energy Classes
Calib. Thawing	$w_0$	—	$\mathcal{U}[-3,1]$
Algebraic Thawing	$w_0$	—	$\mathcal{U}[-1,1]$
	$p$	—	$\mathcal{U}[0,30]$
Emergent	$\Delta$	—	$\mathcal{U}[-3,10]$
Mirage	$w_0$	—	$\mathcal{U}[-3,1]$

models, we have a modified version of the class [89, 90] integrated into our inference pipeline. We switched to the Refast option for recombination as it does not assume anything about the equation of state. We assume one massive and two massless neutrino species with a $\sum m_\nu= 0.06$ eV and $N_{\rm eff}=3.044$. For the SNe Ia likelihoods (PantheonPlus, Union3, and DESY5), we analytically marginalize over the absolute magnitude $M_B$. For clarity of presentation, we utilize Union3 in the figures as a conservative result, as it has larger uncertainties compared to the PantheonPlus and DESY5 datasets. Nevertheless, we will also discuss constraints derived from other supernova datasets wherever they are relevant to our analysis. Finally, $\Delta\chi^2$ is defined with respect to the $\Lambda$CDM best fit. For the calculation of the best fit points themselves, we start with the maximum a posteriori (MAP) points from the four chains produced during the MCMC sampling, and make use of the iminuit [91] minimizer. Thus, the quantity used in model comparisons is more precisely $\Delta\chi^2_{\rm MAP}\equiv-2\Delta\ln\mathcal{L}$, which is the difference in the log posterior values at the calculated maximum posterior points, scaled by $-2$. Since the posterior depends on the product of the likelihood and priors, we also take into account the ratios of different model priors to ensure that there is no additional penalty in the $\Delta\chi^2_{\rm MAP}$ from comparing two models.

Overview of the $w_0w_a$CDM results

We begin by summarizing the dark energy main findings of the DESI DR2 BAO key paper [47], assuming the $w_0w_a$ parametrization given in Eq. (3). As an example, the marginalized constraints in the $w_0w_a$ plane are shown in Figure 1 for the DESI+CMB+Union3 data combination, together with the constraints from DESI DR1 BAO with those obtained with DESI DR2 BAO, corresponding to one and three-years worth of data, respectively. The combined data favor the region $w_0 > -1$ and $w_a < 0$, away from a cosmological constant, implying that the equation of state was phantom-like ($w(z) < -1$) in the distant past and has since evolved to $w(z) > -1$ at present, as shown in the top panel of Figure 2. This preference was observed in previous DESI analyses [40, 41, 42, 43, 44] and persists even when allowing for variations in the spatial curvature ($\Omega_K$) [40], modified gravity [92], or modifications to the pre-recombination physics [93].

DESI DR2 BAO data show that the mean posterior distributions have shifted slightly toward the $\Lambda$CDM-expected values, while the reduced uncertainties have marginally increased the statistical significance of the deviations from $\Lambda$CDM to $2.8 - 4.2\sigma$ (with improvements in fit ranging from $-21.0 \leq\Delta\chi^2 \leq-10.7$), compared to $2.5 - 3.9\sigma$ from DR1 [40, 47]. Similar conclusions follow when using the weighted posterior average of log-likelihood, the Bayesian counterpart of $\Delta\chi^2$. For a detailed Bayesian model comparison, see Appendix A.

Interestingly, with the increased precision, the combined DESI+CMB data already suggest a $\simeq3\sigma$ deviation from $\Lambda$CDM, independent of any SNe Ia compilation, with similar conclusions drawn from the DESI+DESY3 (3$\times$2pt) combination, though exhibiting a lower tension; see Figure 14 in [47].

Physically, a phantom equation-of-state ($w(z) < -1$) translates into an energy density that increases with the expansion ($d\rho_{\mathrm{DE}}/da > 0$), before reaching its maximum (at $z_c \simeq0.45$, in our case), when the equation of state crosses the phantom line ($w(z_c) = -1$), and $\rho_{\mathrm{DE}}$ starts diluting again as the universe expands. The mean redshift at which this transition occurs in the $w_0w_a$ parametrization is indicated by a solid vertical line in Figure 2. We should note at this stage that the exact redshift at which this crossing happens depends on the dataset combination under consideration.

These results may naively suggest a “phantom crossing” [94] at high redshifts. From a theoretical perspective, this so-called phantom crossing is challenging to accommodate within standard scalar-field models of dark energy that are minimally coupled to gravity, as these are constrained to satisfy $-1 \leq w \leq1$. In particular, within general relativity, a single-field dark energy component with $w < -1$ would necessarily violate the null energy condition (NEC), given by $\rho c^2 + P \geq0$ [95]. If confirmed, the phantom crossing would have profound implications for fundamental physics, as it would indicate a significantly more complex dark sector than traditionally assumed. However, it is important to emphasize at this stage that the $w_0$-$w_a$ parametrization is particularly effective at capturing the impact of various, possibly more fundamental, dark energy models on observables such as distances and the expansion history within $\sim0.1\%$ accuracy [54] and may fail to accurately approximate the true behavior of $w(z)$ itself, potentially leading to a spurious indication of phantom crossing. Thus, restricting our analyses to models satisfying $w > -1$ might artificially bias our inference. For more discussion on phantom crossing, see Section VI.E.

Before extending our analysis beyond $w_0w_a$CDM, we introduce two key quantities that will be useful throughout this work. Figure 3 presents the $Om(z)$ diagnostic [96] and the deceleration parameter $q(z)$ for the $w_0w_a$CDM model, where

and the deceleration parameter is given by

These two functions constitute a sensitive probe of new physics, as they are only sensitive to the ‘shape’ of the (normalized) expansion history $h(z) \equiv H(z)/H_0$. Thus, they are unaffected by the degeneracies that may exist between the dark energy and matter densities at the background level [97, 98, 99]. Indeed, one can readily see from Eq. (9) that the quantity $Om(z)$ is strictly constant and equal to present matter density ($Om(z) \equiv\Omega_{\rm m}$) if dark energy is in the form of a cosmological constant. Thus, $Om(z)$ serves as a null-test of $\Lambda$CDM, and any significant deviation from a constant value would indicate dynamics in the dark energy density. The reconstructed $Om(z)$ in Figure 3 shows a clear ($> 2\sigma$) deviation from constancy in the range $0 \lesssim z \lesssim0.5$, where the black dashed line represents the best-fit $\Lambda$CDM value of $\Omega_{\rm m} = 0.302$. On the other hand, $q(z)$ tracks the logarithmic derivative of $h(z)$, rather than its shape, approaching a value of 0.5 during matter domination. The reconstructed $q(z)$ suggests that the Universe’s acceleration ($q < 0$) began earlier in cosmic history ($z_{\rm acc} \simeq0.8$) than predicted by $\Lambda$CDM ($z_{\rm acc} \simeq0.65$) with a slowing down of cosmic acceleration at recent times. These trends in $Om(z)$ and $q(z)$ were previously observed with DESI DR1 data and persist with slightly more statistical significance in DESI DR2.

Parameterizing Dark Energy

To more closely explore the possible dynamical nature of dark energy, we now turn to parameterizations of either the equation of state $w(z)$, or energy density $\rho_{\rm DE}(z)$. Since different parameterizations can lead to differences

in the inferred evolution of dark energy, it is crucial to explore multiple forms to assess the robustness of any detected deviation from a cosmological constant. We examine various two-parameter functional forms as alternatives to $w_0w_a$CDM. In addition, we increase the degrees of freedom available for $w(z)$, to probe the trends present in the data. While the parameterizations investigated here are not necessarily tied to a specific physical model, they cover distinct functional spaces, helping to ensure the results are not driven by the choice of parameterization.

Alternative Lg parameterizations

In this section, we explore four alternative parameterizations from the literature (see [100, 101] for the equivalent DR1 results) that, like $w_0w_a$CDM, introduce two additional parameters: the present-day equation of state $w_0$ and an evolution parameter $w_a$, but with different functional forms.

Figure 4 presents constraints on these alternative models as defined in Table II, i.e.: Barboza-Alcaniz (BA) [102, 103], exponential (EXP)⁴ [104, 105], logarithmic (LOG) [103] and Jassal-Bagla-Padmanabhan (JBP) [106], alongside the Chevallier-Polarski-Linder (CPL) in blue for comparison. The shaded bands, representing $1\sigma$ uncertainties, are derived from a combination of DESI, CMB, and Union3. With the exception of JBP, all models exhibit similar low-redshift behavior, including a phantom crossing near $z\sim0.5$. In Table II, we present the alternative functional forms of $w(a)$ and $\Delta\chi^2$ values relative to $\Lambda$CDM, showing that BA, CPL, LOG, and EXP provide statistically comparable fits to the data. The functional form of the JBP parametrization, which forces it to assume identical early- and late-time behavior, results in a slightly poorer fit. These findings confirm that constraints from CPL are broadly representative of the alternative $w_0w_a$ models considered, with no significant improvement observed for any alternative form. This suggests that current data lacks the sensitivity to distinguish between these parameterizations at $z>2$, a conclusion that remains unchanged across different SNe Ia datasets.

Param.	Functional Form	$\Delta\chi^2$
BA	$w_0+w_a\frac{1-a}{a^2+(1-a)^2}$	$-17.3$
EXP	$(w_0-w_a)+w_a\exp(1-a)$	$-17.5$
LOG	$w_0-w_a\ln a$	$-17.6$
JBP	$w_0+w_a a(1-a)$	$-13.6$
CPL	$w_0+w_a(1-a)$	$-17.4$

Crossing Statistics

Rather than exploring different redshift evolutions for $w(z)$, one can instead gauge the impact of introducing additional degrees of freedom in the DE characteristics. Following the methodology detailed in [43], we expand the equation of state of dark energy $w(z)$ in terms of Chebyshev polynomials (see also [107, 108, 109, 110]),

where $C_i$ are free coefficients, and $T_i(x)$⁵ are Chebyshev polynomials of the first kind, forming a complete basis for continuous functions in the large-$N$ limit, although $N\simeq3$ is generally sufficient to capture smooth functions. We note that $\Lambda$CDM is recovered for $C_0=1$ and $C_{i>0}=0$. Alternatively, one may want to work with the normalized dark energy density $f_{\mathrm{DE}}(z)$ instead

Expanding in $f_{\mathrm{DE}}(z)$ offers the advantage of allowing the (effective) energy density $\rho_{\mathrm{DE}}(z)$ to change sign, thereby encompassing a broader class of models [111, 112, 113, 114], including modified gravity scenarios [115, 116, 117, 118] and complex dark sector interactions that are difficult to capture with a parametrized $w(a)$. We note that the expansion in Eq. (12) has one less degree of freedom relative to that of Eq. (11), as all samples must satisfy $f_{\mathrm{DE}}(z=0)=1$.

The top panel of Figure 5 shows the reconstructed $w(z)$ for the DESI+$({\theta}_*,\omega_b,\omega_{bc})_{\mathrm{CMB}}$ combination, with (blue) and without (orange) SNe Ia data from the Union3 compilation. The bottom panels show the reconstructed $f_{\mathrm{DE}}(z)$. It is noteworthy that not only do the expansions in $w(z)$ and $f_{\mathrm{DE}}(z)$ yield similar behaviors independently of SNe Ia data, but they also agree with the main results obtained using the $w_0w_a$CDM parametrization [47], as shown in Figure 2. This consistency further strengthens the robustness of these results.

While these results perfectly align with the $w_0w_a$CDM results, the expansions in Eqs. (11) and (12) offer greater flexibility, enabling it to capture features in the evolution of dark energy beyond the linear parametrization given by Eq. (3). Despite this additional flexibility—and also confirmed by our independent analyses—the combined data favor a smooth evolution, well described by

$w_0w_a$CDM within the probed low-redshift range. The improvement in fit, quantified by $\Delta\chi^2$, is shown as a function of the number of free parameters in Figure 6. A two-parameter expansion in $f_{\mathrm{DE}}(z)/w(z)$ captures the main trends in the data, as already noted in [119]. Introducing additional degrees of freedom does not significantly improve the fit to the combined data and would be disfavored from a model comparison perspective, as the added complexity is not justified by the data. We note that due to the complications that can arise in the treatment of perturbations when allowing for $\rho_{\mathrm{DE}} < 0$, this part of the analysis is restricted to the “compressed” CMB information, denoted as $(\theta_*,\omega_b,\omega_{bc})_{\mathrm{CMB}}$, rather than the full Planck likelihood. We have verified that $(\theta_*,\omega_b,\omega_{bc})_{\mathrm{CMB}}$, as described in Section II, yields almost identical constraints as those using the primary CMB anisotropies.

While all the parametric models we tested above suggest a phantom crossing, the exact redshift at which this happens depends on the chosen parametrization, as seen from Figure 5 which suggests a slightly higher value for $z_c \simeq0.5$ than $w_0w_a$CDM. This variation—although not statistically significant—is expected due to the inherent limitations of parametric fitting, as each functional form has a restricted degree of flexibility.

Non-Parametric Methods

In contrast to the techniques explored in Section IV, non-parametric techniques focus on determining the true function of quantities such as $h(z)$, $f_{\mathrm{DE}}(z)$, and $w(z)$ from observational data, rather than merely estimating the parameters of a pre-specified form for $w(z)$. We are not interested in model comparison here per se, but rather the robustness of the observed trends in the data under different non-parametric reconstruction techniques.

We explore two techniques: binning and Gaussian process regression. We have also tested the cosmographic expansion up to $O((t-t_0)^5)$ where $t$ is the lookback time, and $t_0$ denotes the current epoch [120, 121]. However, we do not present those results here since they did not pass validation tests with the full DESI data and require a redshift cut-off to make unbiased inference.

Binning

Binning is a technique widely used in cosmology that allows for comparison of different redshift intervals, without the assumption of a specific functional form; see [49, 54, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131] for some examples. Here we focus on binning the equation of state of dark energy, $w(z)$ and the dark energy density, $f_{\mathrm{DE}}(z)$, permitting localized analyses of the behavior motivated by the data. The additional degrees of freedom introduced make it possible to probe for potential variations or trends in $w(z)$ across redshifts, which may help to identify deviations from the standard $\Lambda$CDM model.

In this section, we supplement the 3 uniform redshift bin scheme for the dark energy equation of state parameter from [47] (see Figure 12 there) in order to assess the impact of different choices for the implementation. The binned function takes the general form

where $w_i$ are the bin amplitude parameters, $N$ the number of bins and $s$ the smoothing scale⁶, which controls the sharpness of the transitions around the edges $z_i$ between bins. We assume no prior correlation between bins.

Several different additional schemes for $w(z)$ were tested, including logarithmic binning, binning aligned with the redshifts of the tracer types, and various uneven binning approaches across the constrained redshift interval. However, for clarity, we present results only for schemes with uniform redshift bins between $z=0$ and $z=2.1$, as results do not change qualitatively across the different binning schemes. We consider the combination of DESI, Union3 SNe Ia, and CMB. In the case of $f_{\mathrm{DE}}(z)$, the compressed CMB is used to avoid the computational complexity associated with correctly modifying the behavior of dark energy in a Boltzmann solver to account for $f_{\mathrm{DE}} < 0$, while also constraining the parameters exclusively with early CMB information.

Figure 7 (upper panel) shows the median values of $w_i$, with $1\sigma$ and $2\sigma$ error bars, positioned at the center of their respective bins’ redshift intervals. These intervals are shown in the same colors in between the panels. The highest redshift interval effectively extends to high redshifts, with the corresponding amplitude positioned at $z = 2.8$ merely for convenience. The constrained amplitudes for overlapping bins between different schemes are all within $\sim1\sigma$ of one another. Superimposing with the median (dashed grey line) and $1\sigma$, $2\sigma$ confidence levels of the corresponding $w_0w_a$CDM result, we see that the behavior recovered by each different scheme is in good general agreement, with median points on either side of the line $w(z) = -1$. The data provide the tightest constraints in the lowest redshift bin, where they prefer a $w(z)$ that is more than $3\sigma$ away from $\Lambda$CDM value of $-1$, whereas the higher redshift bin amplitudes remain, at most, within $2\sigma$ of $\Lambda$CDM. The question of an actual crossing is more subtle, since it would have to occur at the edge between two adjacent bins, meaning that it would depend non-trivially on the number of bins and their chosen centers, and not make for a very robust ‘measurement’.

To allow explicit exploration of the region of parameter space with negative $f_{\mathrm{DE}}(z)$, which is excluded when binning with the amplitude of $w(z)$, we also test additional binning schemes where the bin variables are instead associated with the amplitude of $f_{\mathrm{DE}}(z)$. Figure 7 (lower panel) shows the same effective behavior between individual binning schemes, with good agreement to the $w_0w_a$ curves, indicating a turn-over somewhere in the region of $0.5 < z < 1.0$ and $f_{\mathrm{DE}}(z) > 0$ at around $2\sigma$ for most of the bins. The uncertainties increase with redshift, becoming progressively less Gaussian, with longer tails extending towards more negative $w$ values. Lastly, we note that the amplitudes in adjacent bins exhibit mild correlations, weakening with increasing redshift.

To decorrelate the bins and get additional insights into the contributions from different redshift intervals, we also perform a principal component analysis. Principal component analysis (PCA) is effectively a transformation that provides a new basis in which the new coefficients $q_i$, corresponding to the bin amplitude parameters, are decorrelated. There are, in general, infinitely many such decorrelated bases, but only one that is orthogonal. We may obtain it simply by finding the eigenvector basis that diagonalizes $C^{-1}$, the inverse covariance matrix of the bin amplitude parameters $w_i$, calculated by MCMC sampling [52]. See Appendix B for additional details.

We divide the equation of state parameter into 10 uniform bins of fixed amplitude between $z = 0.1$ and $z = 2.1$, with two additional free parameters, one each for the amplitude on either side. The covariance matrix of the resulting bin parameters $w_i$, $i = 0,1,\ldots,11$ with DESI+CMB+Union3 is used to determine the eigenvector basis. The basis functions, or principal components, corresponding to the 4 largest eigenvalues are presented in Figure 8, along with the corresponding errors (obtained as square roots of inverse eigenvalues)

The largest principal component is well localized in $z$, peaking in the range $0.1 < z \leq0.3$, while the second-largest component is mostly positive and peaks in the interval $0.3 < z \leq0.5$. The remaining components show increasingly more pronounced oscillatory behavior, with

at least $20\times$ the uncertainty of the first bin, $\sigma_0$. Overall, the binning results from different schemes are in good general agreement. The crossing of phantom divide line by $w(z)$, and turn-over in $f_{\mathrm{DE}}(z)$ followed by a decreasing trend towards higher redshifts, found in the other analyses [40, 43, 132, 133] are consistent with these results. Even so, the approach has its limitations. While it is well-suited to testing deviation from a constant function, capturing more complicated behaviors requires additional degrees of freedom, which increases the level of uncertainty [52, 54, 124, 125, 128]. In particular, though the data seem to be consistent with a phantom crossing, it is difficult to draw strong conclusions about the specific redshift where a phantom crossing of $w(z)$, or turn-over in $f_{\mathrm{DE}}(z)$, might occur. The limitations present in this approach make it important to understand what kind of biases may be introduced by the implementation. In Appendix E, we perform some tests on simulated data in an attempt to address this.

Gaussian Process Regression

In this section, we discuss Gaussian processes, which can be thought of as a generalization of binning where the amplitudes at every redshift are sampled but are subject to some constraints (prior assumptions) on the form of the resulting functions. This allows for a complementary analysis with the possibility of improving the trade-off between flexibility and constraining power.

Gaussian Process (GP) regression [135] is a powerful, non-parametric statistical tool widely used in various fields, including cosmology [136, 137, 138], to reconstruct smooth functions from noisy data without assuming a specific functional form (see e.g., [62, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148] for a non-exhaustive list). For the purpose of this work, GP can be thought of as a way of sampling the space of continuous functions in a non-parametric manner. This allows data-driven reconstructions of the quantities of interest for dark energy, namely $w(z)$ or $f_{\mathrm{DE}}(z)$ with minimal assumptions [62, 136, 149, 150]. More specifically, at each point in parameter space, we draw a sample (realization) of $w$ from a multivariate Gaussian distribution, e.g., $w(z) \sim\mathrm{GP}(m(z) = -1, K)$ where $K$ is a given covariance function — known as kernel — encoding our prior assumptions about the smoothness of the reconstructed function. We further impose $w(z \geq z_{\max}) = -1$ to recover a standard ($\Lambda$CDM-like) expansion history at early times. We have chosen $z_{\max} = 10$, after checking that this choice does not significantly alter our conclusions. This is implemented in a modified version of the Boltzmann solver camb, with more details on GP given in Appendix C.

Figure 9 illustrates the reconstructed dark energy properties using GP with various datasets: DESI + CMB (left), DESI + Union3 (middle), and DESI + Union3 + CMB (right). The top row presents the reconstructed $w(z)$, indicating deviations from $\Lambda$CDM at low redshift and hints of $w(z)$ crossing into the phantom regime around $z \approx0.5$. The inclusion of CMB data (left and right panels) results in tighter constraints on $\Omega_m$, which strengthen the significance of deviations from $w = -1$, whereas the DESI + Union3 combination allows for lower values of $\Omega_m$ and a broader range of variations $w(z)$. The second row demonstrates a notable bump in the evolution of the dark energy density, while constraints derived without CMB allow for a wider variety of $w(z)$. The third row displays the normalized Hubble parameter $h(z)/h_{\Lambda\mathrm{CDM}}(z)$. The fourth row presents the $\Omega_m$ diagnostic, clearly showing the evolution as a function of $z$, indicating a deviation from $\Lambda$. Lastly, the final row depicts the reconstructed deceleration parameter $q(z)$, which slightly exceeds the expectations of the $\Lambda$CDM model, suggesting a slowdown in the acceleration rate.

In Figure 10, we present a comparison of results obtained from GP reconstruction utilizing the $w_0w_a$ parametrization derived from DESI, CMB, and Union3 data. The GP reconstruction, illustrated in blue, aligns very well with the $1\sigma$ posterior predictions of the $w_0w_a$CDM model. We would like to remind readers that the GP approach imposes a Gaussian prior distribution on $w(z)$, centered at the mean function which we explicitly choose to be $w = -1$, as represented by the black dotted line. This effectively places more prior weight on $\Lambda$ and any observed deviations from $w(z) = -1$ are largely driven by data. Finally, we would like to emphasize that although GP regression offers advantages over parametric methods, it is important to interpret the reconstructed $w(z)$ with caution. While flexible, the method may not fully capture certain behaviors of $w(z)$, as illustrated in Appendix E using simulated data. Nevertheless, GP remains a valuable tool for assessing the dynamical nature of dark energy in a non-parametric manner.

Implications for Dark Energy

The various methods explored in Sections IV and V provide a flexible way to test deviations from $\Lambda$CDM and ensure robust results without committing to a specific dark energy model. However, interpreting the deviations from a cosmological constant and understanding its implications for fundamental physics necessitates a deeper exploration of physically-motivated dark energy models. Rather than constraining specific models, we focus on different classes, characterized by their dynamics [151, 152, 153] and inspired by theoretical considerations.

Thawing dark energy

The first class of models we consider is known as thawing dark energy [151]. This class characterizes quintessence models [154, 155, 156, 157], in which a minimally coupled scalar field remains frozen at early times due to Hubble friction, effectively behaving like a cosmological constant with $w=-1$. Only when the scalar field’s mass becomes comparable to the Hubble rate, $m_\varphi\sim H$, does the field begin to evolve dynamically, causing its equation of state to ‘thaw’ away from $w=-1$ into the quintessence regime, $w>-1$. Note that there exists a second class of DE dynamics, referred to as the ‘freezing’ class, where the field evolves towards a de Sitter state ($w=-1$) in the asymptotic future. Such dynamics are characterized by $w_a>0$ and are not favored by observations. For a review on quintessence models, see e.g., [158, 159]

This behavior is typical of pseudo-Nambu-Goldstone-Boson (PNGB) quintessence models [160] and simple potentials such as $V\propto m^2\varphi^2$ and $V\propto\lambda\varphi^4$, both of which are ubiquitous in high-energy physics [24]. Interestingly, Ref. [54] demonstrated that the phase-space dynamics of these models can be well-approximated using the $w_0w_a$ parametrization. Many thawing potentials map onto a narrow region in the $w_0w_a$ plane, approximately following the relation

This ‘calibrated thawing’ relation provides a form that acts as a good approximation for the thawing dynamics. However, it also allows the equation of state to cross the phantom line ($w=-1$), which is unphysical for quintessence models [94, 161, 162]. This occurs because Eq. (14) is designed to approximate the expansion rate $H(z)$ and distance measures $D(z)$ at sub-percent precision — precisely the quantities probed by cosmological observations — but does not necessarily approximate $w(z)$ itself [54].

Nevertheless, it is possible to describe thawing dynamics while ensuring that $w>-1$ at all times. Following [163, 164] (see also [165]), the evolution of the thawing equation of state can be parameterized by the algebraic expression

where $p$ and $w_0$ are free parameters, and $b=0.5$ [164]. Notably, this formulation, referred to as ‘algebraic thawing’ is more general, where the case $p=1$ and Eq. (14) were found to yield nearly identical late-time constraints as shown in Appendix A of [44].

The reconstructed posterior distribution of $w(z)$ for the thawing class is shown in the top panel of Figure 11 for the DESI+CMB+Union3 data combination. This assumes the algebraic form, which enforces $w\geq-1$, and where we have marginalized over the parameter $p$. However, mild degeneracies with $w_0$ leave the posteriors for $p$ largely unconstrained. In particular, it is seen that large values of $p$ are allowed by the data, resulting in our posterior hitting the prior bound $p=30$, as shown in the bottom panel of Figure 11. However, we do not extend our analysis to larger values of $p$, as numerical complications can arise when dealing with DE models with very rapidly varying $w(a)$, particularly in the treatment of perturbations and CMB lensing.

The calibrated thawing relation Eq. (14) yields no significant improvement in fit, as seen from the $\Delta\chi^2$ values in Table III and from the posteriors in Figure 13 being consistent with $w_0=-1$, except for DESY5. The algebraic thawing parametrization in Eq. (15) can improve

DE classes	+PantheonPlus $\Delta\mathrm{DIC}$ ($\Delta\chi^2$)	+Union3 $\Delta\mathrm{DIC}$ ($\Delta\chi^2$)	+DESY5 $\Delta\mathrm{DIC}$ ($\Delta\chi^2$)
Thaw. (Cal.)	+0.4 (−1.6)	−0.6 (−2.5)	−5.8 (−7.1)
Thaw. (Alg.)	−1.0 (−2.9)	−4.6 (−6.9)	−10.1 (−13.2)
Emergent	+2.1 (−0.05)	+1.8 (−0.1)	+0.2 (−1.5)
Mirage	−9.1 (−10.5)	−13.8 (−16.2)	−18.7 (−20.7)
$w_0w_a$	−6.8 (−10.7)	−13.5 (−17.4)	−17.2 (−21.0)

the fit with respect to the standard model ($\Lambda$CDM), achieving a $\Delta\chi^2=-2.9$ for the DESI+CMB+Union3 data combination. Substituting the PantheonPlus data with Union3 or DESY5 yields a $\Delta\chi^2=-6.9$ and $-13.2$, respectively. However, this improvement in fit comes at the cost of including two additional degrees of freedom.

To illustrate that Eq. (15) correctly captures the phenomenology of thawing fields, and better quantify how the constraints would translate into constraints on the physical parameters of the theory, we consider one physically motivated, axion-like potential [160, 166, 167]

where $m_a$ is the mass of the boson particles related to the scalar field, and $f_a$ is regarded as the effective energy. Depending on the initial conditions, the axion cosine potential exhibits two distinct behaviors: the standard quadratic regime, the effective mass is positive ($m_{\rm eff}^2>0$), and the potential can be approximated by a quadratic form near its minimum, where the effective mass is defined as $m_{\rm eff}^2=\frac{\partial^2V}{\partial\varphi^2}$. Whereas, the hilltop regime [168] is, characterized by a negative effective mass $(m_{\mathrm{eff}}^2<0)$, when the field begins its evolution near the maximum of the potential (i.e., at $\varphi=0$) and rolls down toward the minimum at $\varphi=\pi f_a$. We refer the reader to Appendix D for more details on the model and its implementation in class.

In Figure 12, we report the marginalized posterior distribution for the equation of state parameter associated with the scalar field potential in Eq. (16), obtained using DESI, CMB, and three SNe compilations and obtain the following constraints for the physical mass: $\log_{10}(m_a c^2/\mathrm{eV})=-32.67^{+0.23}_{-0.25}$ (PantheonPlus), $-32.50^{+0.28}_{-0.30}$ (Union3) and $-32.63^{+0.26}_{-0.30}$ (DESY5) and effective energy scale: $\log_{10}(f_a/M_{\mathrm{Pl}})=-0.13^{+0.33}_{-0.29}$ (PantheonPlus), $-0.29^{+0.63}_{-0.35}$ (Union3) and $-0.09^{+0.66}_{-0.40}$ (DESY5). The constraints indicate that the field starts in the hilltop regime, with initial conditions of $\varphi_i/f_a\sim0.7-1.0$, rolls down the potential, and reaches the present value of $\varphi_0/f_a\sim1.1-1.4$, traversing approximately $\Delta\varphi\sim(0.2-0.4)M_{\mathrm{Pl}}$.

Emergent dark energy

The second family of DE models that we consider is the emergent class, where dark energy had a vanishing presence during most of cosmic history, and only ‘emerges’ in recent times. Following [169, 170], we parametrize the equation of state as

The parameter $\Delta$ determines the steepness of the transition in $w(z)$ and the transtion redshift parameter $z_t$ is determined by the equality $\rho_{\mathrm{DE}}(z_t)=\rho_m(z_t)$. The phenomenology we are trying to capture is that of abrupt changes in the equation of state, $w(z)$, driven by physical mechanisms, such as second-order phase transitions [171, 172, 173].

Despite hints of the sharp emergence of dark energy in recent times from non-parametric reconstructions, the DESI+CMB+SNe constraints on $\Delta$, as shown in the middle panel of Figure 13, indicate that such an emergent behavior is not statistically favored over $\Lambda$CDM, given the assumed $w(z)$. Note that while Eq. (17) can mimic the emergence of dark energy, it is limited by its inability to cross $w=-1$ or, equivalently, introduce a bump in $f_{\mathrm{DE}}(z)$; a feature that seems to be favored by the data. In principle, one can formulate an emergent dark energy model characterized by an effective equation of state that can cross $w(z)=-1$. Such behavior may be realized through the coupling of emergent dark energy with the dark matter sector [174, 175, 176].

Mirage dark energy

The last and more phenomenological class of models which we consider is that of mirage dark energy [177]. This refers to models in the $w_0w_a$ plane (see Figure 1) approximately living along the line

The mirage class is designed to describe a subset of dynamical dark energy models that preserve the distance to the surface of the last scattering as predicted by $\Lambda$CDM, a parameter tightly constrained by the CMB [177]. The name ‘mirage’ stems from the fact that these models would mimic $\Lambda$, yielding $\langle w\rangle\sim-1$ when fitting a constant $w$ to observations, as it could be seen in Table V [47] in DR2 and Table III in [40] for DR1 comparison. The mirage direction fully captures the DE phenomenology suggested by the data, with merely one degree of freedom $w_0$ that quantifies the strength of the mirage, with $w_0=-1$ corresponds to $\Lambda$CDM where the mirage is real. This mirage effect is also expected to persist in the growth of cosmic structures, provided that general relativity remains unmodified [177, 178, 179]. As noted in [44], by reducing the late-time dark energy density (i.e., increasing $\Omega_m$), one can make $w_0$ even less negative — and correspondingly $w_a$ more negative — enhancing the mirage effect. For comparison, DESI+CMB+Union3 prefers $\Omega_m\approx0.304$ in $\Lambda$CDM which increases to $\Omega_m\approx0.327$ in $w_0w_a$ which essentially lies along the mirage direction. From the data viewpoint, the mirage line in Figure 1 can be seen as the ‘principal component’, or ‘axis’ in the $w_0w_a$ plane carrying the most meaningful information, i.e., the eigenvector with the highest eigenvalue. Despite effectively reducing the dimensionality of the DE phenomenology, the exact physical mechanism for such rapid emergence of dark energy ($w(a)\ll-1$) remains unclear (see [180] for more discussion).

Model Comparison

In Figure 13, we show the constraints on the single additional parameters of the three DE classes. In the calibrated thawing (top panel), there is a mild deviation observed with the DESY5 dataset, but the other two SNe Ia datasets indicate overall consistency with ΛCDM. The emergent (middle panel) exhibits a similar trend, with constraints on $\Delta$ reflecting no significant departures from the standard cosmological framework ($\Delta= 0$). In contrast, the mirage class (bottom panel) demonstrates deviations from the value of $w_0 = -1$ across all three SNe Ia datasets.

We also make a quantitative comparison, examining the deviance information criterion (DIC)[181, 182], as defined below, along with the $\Delta\chi^2$. The former complements the $\Delta\chi^2$ by accounting for model complexity, which does not take into account the number of additional degrees of freedom in a particular model and could be made arbitrarily low if sufficient parameters were added. The DIC is defined as

where $p_D = \overline{D(\theta)} - D(\bar{\theta})$ is a penalty term, and $D(\bar{\theta}) = -2\ln\mathcal{L} + C$ is the ‘deviance’ of the likelihood, with constant $C$ vanishing in $p_D$. In practice, we use

and $p_D$ becomes effectively equivalent to the number of extra parameters in the limit of parameters that are well-constrained with respect to their prior. We consider the DE classes in this section and the $w_0w_a$ parametrization, for the data combinations DESI+CMB+SNe Ia. The comparisons are made between each model class and the ΛCDM model, and the key metric for each comparison being the DIC. The $\Delta\mathrm{DIC}$ (and $\Delta\chi^2$) values are reported in Table III, with a preference for the more complex model indicated by negative values, and for the simpler model, in this case ΛCDM, by positive values. For nested models, a decrease ($\Delta\mathrm{DIC} < 0$) of at least 2 is required for a not ‘insignificant’ improvement, and up to 5 constitutes a ‘positive’ preference over ΛCDM. A decrease of up to 10 is considered a ‘strong’ preference, and beyond this, the preference is ‘decisive’ [183]. However, for classes that are not nested within each other, there is no absolute scale for comparison and can only be quantitatively compared against ΛCDM.

The $w_0w_a$ parametrization achieves $-\Delta\chi^2 \sim10.7 - 21.0$ and $-\Delta\mathrm{DIC} \sim6.8 - 17.2$, indicating that it is strongly to decisively preferred over the standard ΛCDM. Comparatively, the calibrated thawing performs poorly. The algebraic thawing class improves the fit slightly more, and since the $\Delta\mathrm{DIC}$ is consistently larger than for calibrated thawing, the improvement in $\chi^2$ must be sufficient to reconcile the former’s second additional degree of freedom. The emergent dark energy class shows no significant improvement in fit, faring even worse than the thawing class. The final mirage class attains $-\Delta\chi^2 \sim10.5 - 20.7$, comparable to the $w_0w_a$ parametrization, as well as $-\Delta\mathrm{DIC} \sim9.1 - 18.7$. This is perhaps unsurprising, given how closely the mirage direction aligns the $w_0w_a$CDM constraints.

The advantage of the thawing and emergent classes is their connection to a physical interpretation, which is fairly straightforward, while for the mirage class less so. It is also important to note that the model comparison metrics used in this section serve both to quantify the data’s preference for each model—providing an absolute scale in the case of nested models—and to facilitate a relative ranking among non-nested models, such as the algebraic thawing and $w_0w_a$ parameterizations.

Is there evidence for phantom crossing?

Both the parametric and non-parametric methods discussed in Section IV and Section V indicate a possible crossing of the phantom divide line; however, as illustrated in Figure 16, this does not guarantee that the crossing is genuine. Specifically, the apparent crossing in the $w_0w_a$ parameterization may be spurious to match observables, raising the question of whether $w(a)$ truly crosses $-1$ or if the behavior is simply an artifact of the parameterization.

To address this, we analyze the behavior of thawing quintessence using the algebraic model described by Eq. (15), which enforces $w(z) \ge-1$. Although the algebraic thawing model, restricted by our prior $p < 30$, yields a better fit to the data than $\Lambda$CDM by $3 \lesssim-\Delta\chi^2 \lesssim13$ and $-1 < -\Delta\mathrm{DIC} < 10$, it is considerably less favored than the $w_0w_a$ model. This overall preference for $w_0w_a$CDM over algebraic-thawing model, however, cannot be straightforwardly converted into $p$-values or $n$-sigma levels because algebraic thawing and $w_0w_a$CDM models are not nested.

To achieve $-\Delta\chi^2$ comparable to $w_0w_a$ with thawing models would require an exceptionally fine-tuned potential $V(\phi)$ [184], precise initial field settings [185], or ‘data-informed’ prior choices [186], resulting in $w(z \gtrsim0.3)=-1$ followed by a rapid increase to $w(z \lesssim0.3)>-1$. This suggests that a sharp increase followed by a decrease in dark energy density may be a necessary feature, since models that do not cross the phantom divide tend to underperform compared to those that do. The specific behavior of $w(z)$ suggested by the data — a phantom crossing from $w(z \gtrsim0.5)<-1$ to $w(z \lesssim0.5)>-1$ — is not predicted by any of the simplest and most studied extensions of $\Lambda$CDM (see also the recent discussion in [180]).

While the phantom crossing can be seen as theoretically challenging due to stability issues, for example, within the framework of minimally-coupled scalar fields, obtaining such behavior is not difficult in extended theoretical frameworks. For example, phantom crossing can arise in models where dark energy possesses multiple internal degrees of freedom, such as multi-field scenarios [187, 188, 189, 190, 191, 192], non-standard vacuum models [171, 172], frameworks where dark energy interacts with dark matter [193, 194, 195, 196, 197, 198], and modified theories of gravity [199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209]. Because of the aforementioned multiple internal degrees of freedom in these models, the effective, observable equation of state $w(z)$ can cross the phantom divide even though the null-energy condition is not violated. Whether a compelling theoretical mechanism — one that does not require many extra degrees of freedom or exotic assumptions — can be constructed to cross the phantom divide in the way suggested by data remains an open question, although some models have recently been put forward as viable in this regard [206, 207].

Conclusions

This work presents constraints on dark energy from DR2 BAO, in combination with cosmic microwave background, and type Ia supernova data. We began our study by summarizing and expanding upon the $w_0w_a$CDM analysis presented in [47]. Under the assumption of a linearly evolving $w(a)=w_0+w_a(1-a)$, the latest DESI+CMB results indicate a $\simeq3\sigma$ deviation from $\Lambda$CDM. The data shows a clear preference for the $(w_0>-1,w_a<0)$ quadrant, and in particular $w_0+w_a<-1$, implying a past phantom-like equation of state transitioning to $w(z)>-1$ today. The reconstructed $Om(z)$ and deceleration parameter $q(z)$ also show clear deviations from $\Lambda$CDM, reinforcing the case for evolving dark energy.

To assess the robustness of our findings, we conducted a series of analyses: i) varying the redshift dependence of $w(z)$, by considering various parameterizations (Section IV A), and ii) studying the improvement in fit as more freedom is given to the dark energy characteristics (Section IV B). As in DR1, the results are rather stable under changes in the assumed form for $w(z)$, and the data does not seem to require more degrees of freedom in $w(z)$, beyond $w_0w_a$, as shown in Figure 6 (see also [43]).

Next, we implemented two non-parametric reconstruction techniques and applied them to the redshift-dependent equation of state $w(z)$ and dark energy density $f_{\mathrm{DE}}(z)$, in order to allow more flexibility than that available in the parametric methods. Overall, the conconstraints support the evolution indicated by the $w_0w_a$ parametrization, giving the tightest constraints at low redshifts, where they display a preference for a deviation with $w(z) > -1$, while suggesting a crossover to the phantom regime at higher redshift. The low redshift deviation is evidently independent of the chosen binning variable, although the constraints remain within $2\sigma$ of $\Lambda$CDM at higher redshifts. Gaussian process regression is better able to localize the redshift where the crossing should occur, around $z \sim0.5$.

In order to provide possible interpretations for the physical origin of the observed deviation, three model classes were considered, each endowed with a different dynamical behavior and motivated to various degrees by physical theory. The thawing and emergent classes are the least well-supported, indicating that the data might not favor dark energy evolution that arises from, respectively, either minimally-coupled scalar field models or emergent behavior in energy density. In contrast, the mirage class performs remarkably well, capturing DE phenomenology with just one additional degree of freedom, which warrants an inquiry into whether any underlying physics or systematic effects could explain this mirage.

In summary, irrespective of the parametric/non-parametric methods used, the evidence of deviation from $\Lambda$CDM is significant. Our findings suggest that the canonical $w_0w_a$ parametrization effectively captures the essence of dark energy evolution in our study. Decisive tests of dark energy and its possible deviations from the $\Lambda$CDM model will require a combination of complementary probes. The forthcoming DESI data releases, including constraints from redshift space distortions and peculiar velocities, will offer crucial insights into the nature of dark energy and gravity. The upcoming SNe measurements from the ZTF survey [210, 211], the Vera C. Rubin Observatory [212, 213] and Nancy Grace Roman Space Telescope [214] will extend the Hubble diagram probed by DESI to very low redshifts, improving constraints on $w_0$. Meanwhile, data from Euclid [215] and Rubin will serve as an important cross-check of DESI’s findings, helping to assess the impact of potential systematics. Finally, next-generation CMB experiments will further tighten constraints on early-universe parameters, breaking degeneracies with late-time observables. With these advancements, the next decade promises to determine if we are entering a new era in modern cosmology that necessitates a paradigm shift.

Data Availability

The data used in this analysis will be made public along the Data Release 2 (details in https://data.desi.lbl.gov/doc/releases/). The data points corresponding to the figures from this paper will be available in a Zenodo repository.

Acknowledgments

The authors thank Eric Linder for his valuable discussions and Robert Crittenden and Kazuya Koyama for their detailed comments. R.C. is funded by the Czech Ministry of Education, Youth and Sports (MEYS) and European Structural and Investment Funds (ESIF) under project number CZ.02.01.01/00/22.008/0004632. A.S. would like to acknowledge the support by National Research Foundation of Korea 2021M3F7A1082056. MI acknowledges that this material is based upon work supported in part by the Department of Energy, Office of Science, under Award Number DE-SC0022184, and also in part by the U.S. National Science Foundation under grant AST2327245. CGQ acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51554.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

This material is based upon work supported by the U.S. Department of Energy (DOE), 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 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 Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions.

The DESI Legacy Imaging Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall z-band Legacy Survey (MzLS). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory. Legacy Surveys also uses data products from the Near-Earth Object Wide-Field Infrared Survey Explorer (NEOWISE), a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Legacy Surveys was supported by: the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility; the U.S. National Science Foundation, Division of Astronomical Sciences; the National Astronomical Observatories of China, the Chinese Academy of Sciences and the Chinese National Natural Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. The complete acknowledgments can be found at https://www.legacysurvey.org/.

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.

The authors are 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.

References

Bayesian Model Comparison

Here we discuss the Bayesian model comparison and compute Bayes factor between the $w_0w_a$ and $\Lambda$ models, which is given by the corresponding ratio of their evidence $\mathcal{Z}$ under a given data set: $B_{w_0w_a\Lambda}=\mathcal{Z}_{w_0w_a}/\mathcal{Z}_{\Lambda}$. In practice, we compute Bayes factors using the nested sampler polychord [216], employing the same priors⁷ and the Boltzmann solver (CAMB) as in the posterior analysis discussed in the main text. The evidence is then estimated using the anesthetic package [217].

The logarithm of the evidence (log $\mathcal{Z}$) can be expressed as the contribution from two terms in the form [218]

where, $\mathcal{L}$ is the likelihood, $\mathcal{P}$ is the posterior, $\pi$ represents the prior distribution and $\langle\rangle_{\mathcal{P}}$ is the posterior weighted average. The posterior average of the log-likelihood $\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$ removes the prior-dependent Occam’s penalty, $\langle\ln(\mathcal{P}/\pi)\rangle_{\mathcal{P}}$, contribution from the log-evidence to provide a quantitative assessment of how well the model fits the data and can be considered the Bayesian equivalent of $\chi^2$.

In Table IV, we report the differences between the $w_0w_a$CDM and $\Lambda$CDM models, for values of $\ln B$ (= $\Delta\ln\mathcal{Z}$), $\Delta\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$, and $\Delta\langle\ln(\mathcal{P}/\pi)\rangle_{\mathcal{P}}$ for each data combination. The Bayesian evidence ratio indicates that the support for the $w_0w_a$ model increases with the latest DR2 release. On Jeffreys’ scale [219, 220], the DESY5 combination indicates almost a strong preference ($\leq5$) for the $w_0w_a$CDM model when combined with DR2 compared to a moderate preference ($\leq2.5$) with the DR1 dataset, while Union3 shows a moderate preference for $w_0w_a$ compared to the weak preference observed when combined with the DR1 dataset. In contrast, Pantheon-Plus provides inconclusive evidence, showing a preference for $\Lambda$, though this preference diminishes in the change from DR1 to DR2. In all three cases, we observe that the trends in $\Delta\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$ remain consistent with the frequentist $\Delta\chi^2$ results presented in Table III, with a clear improvement in the likelihood of the fit when switching from the DR1 to the DR2 dataset. We note that since the Bayesian evidence value includes likelihood contributions across a range of possible parameter values, it is less susceptible to random fluctuations in the data leading to a better fit at a specific point in parameter space by chance. Therefore, in $B$ or $\Delta\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$ comparisons should be less noisy and more robust than a $\Delta\chi^2$.

	+PantheonPlus	+Union3	+DESY5
DESI DR2 BAO + CMB
$\ln B_{w_0w_a\Lambda}$	$-0.66 \pm0.44$	$3.03 \pm0.44$	$4.88 \pm0.44$
$\Delta\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$	$4.53 \pm0.17$	$7.78 \pm0.16$	$9.42 \pm0.16$
$\Delta\langle\ln(P/\pi)\rangle_{\mathcal{P}}$	$5.18 \pm0.43$	$4.75 \pm0.43$	$4.54 \pm0.44$
DESI DR1 BAO + CMB
$\ln B_{w_0w_a\Lambda}$	$-1.54 \pm0.44$	$1.62 \pm0.44$	$2.35 \pm0.44$
$\Delta\langle\ln\mathcal{L}\rangle_{\mathcal{P}}$	$3.34 \pm0.17$	$6.03 \pm0.16$	$7.67 \pm0.16$
$\Delta\langle\ln(P/\pi)\rangle_{\mathcal{P}}$	$4.87 \pm0.43$	$4.41 \pm0.43$	$5.32 \pm0.43$

A larger Occam’s penalty (i.e., a greater compression from prior to posterior) is typically expected with more informative data. However, this is not the case for the DESY5 dataset, although it is a better fit for the data. Even so, all of the Occam penalty factors for DESI +

CMB + SNe Ia data are roughly consistent with each other, given the estimated error on these values.

Please note that in the [40] used the Boltzmann solver class with a different set of priors and CMB likelihood. For a better comparison, we recomputed evidence for DESI DR1 BAO using consistent priors and methodology as described in the Section II.

We remind readers that Bayesian model comparison relies on the prior, this is especially important in the case of testing phenomenological models like $w_0w_a$CDM, as the prior chosen for the extra parameters is only phenomenologically justified. Occam’s penalty can be adjusted by choosing different priors. Generally, a wider prior tends to favor the simpler model, in this case is the $\Lambda$CDM. However, the priors chosen in Tab. I are broad enough for $w_0w_a$, and the posterior of these parameters is well constrained inside this prior when the datasets are combined. A different choice of prior on these parameters could change the ranking, such that $\Lambda$CDM would be favored over $w_0w_a$CDM for all data combinations, but to do this in the most extreme case (DESI DR2 BAO + CMB + DESY5) the prior range would need to be expanded by a factor of more than 10 times for both $w$ parameters, leading to a greater than 100 times expansion in prior volume. Such a prior would not doubt be considered unphysical and unreasonable by most.

Details of Binning PCA

Principal component analysis (PCA) is a commonly-used strategy that leverages the number of bins, or additional free parameters, against the uncertainties in the constrained amplitudes, in an attempt to analyse how most efficiently to segment the data [52, 54, 123, 124, 128], though the interpretation is not necessarily straightforward, and it remains subject to various caveats [54].

In this appendix, we lay out the full mathematical expressions for the PCA performed on the binning results. We start by writing the equation of state in terms of the binning amplitude parameters and an initial basis

where $w_j$ are the bin amplitude parameters, defined as previously, and $e_j$ is simply the tanh smoothed top-hat function

which has value one in the interval $z_{j-1} \lesssim z \lesssim z_j$ and zero elsewhere, allowing us to recover the expression in Eq. (13).

If we diagonalize the covariance matrix of these coefficients to obtain

then the orthogonal matrix $O$ contains the eigenvectors of $C^{-1}$ and the diagonal matrix $\Lambda$ contains the corresponding eigenvalues [52, 54]. For convenience, we normalize $O$ to have a determinant of one.

To obtain the new coefficients, we use the rows of $O$ as weights on the originals

Since $O$ is orthogonal, i.e. $O^{-1}=O^{T}$, the coefficients and basis will transform in the same way. As such, the new basis functions (see Figure 8) may be obtained in a similar fashion using [52, 54]

The uncertainty in the new parameters is given by the inverse of the eigenvalues [52]

While the PCA is no doubt a useful approach, it must be noted that, in general, eigenvectors are formally not well-defined for the inverse covariance matrix, and that the set of eigenvectors found will themselves depend on the binning parametrization and variable. Finally, there is no clear a priori interpretation of the size of $\sigma[q_i]$ without making further assumptions about the form of the equation of state [54].

Details on Gaussian Process Regression

In this appendix, we discuss the details of the methodology and implementation of the Gaussian process regression results presented in § 2.1 B. We adopt a flexible and non-parametric approach to model the dark energy equation of state $w(z)$ utilizing a Gaussian Process with a squared exponential kernel given by

where $\sigma_f$ denotes the typical deviations of $w$ from the mean function, which we consider to be the $\Lambda$CDM value $w(z)=-1$. The parameter $\ell_f$ controls the correlation length of samples. We draw samples from a joint Gaussian distribution across a dense uniform grid ranging from $z=0$ to $z=10$ that smoothly go back to $w(z)\to-1$ at high redshifts to maintain numerical stability. We fix $\sigma_f=1$ for accommodating a broad spectrum of dark energy behaviors and impose a generalized inverse Gaussian prior on $\ell_f$:

where $K_p(\sqrt{b})$ denotes the modified Bessel function of the second kind. We select $p=3$ and $b=2$ to penalize low $\ell_f$ values, thus preventing excessive freedom in $w(z)$, while also constraining high $\ell_f$ values to avoid oversampling nearly linear functions. To improve the efficiency of our sampling procedure, we introduce a latent variable that directly samples $w(z_l)$ at $z_l$ (we fix $z_l=0.4$) from a Gaussian prior with width $\sigma_f$. These sampled values are subsequently employed to generate samples of $w(z)$ from a conditional distribution, facilitating significantly quicker convergence. We have assessed the impact of varying $z_l$ (over the range $z\in[0.2,1.5]$) and marginalizing $\sigma_f$ with a uniform prior ($\pi(\sigma_f)=\mathcal{U}[0,2]$), finding our results and conclusions to be reasonably robust.

Quintessence

We discuss the constraints on a quintessence model with a scalar field $\varphi$ minimally-coupled to gravity, with action given by (e.g. [221])

where $g$ is the determinant of the metric $g_{\mu\nu}$. For our purposes, the scalar field potential has a periodic dependence on the field $\varphi$ in the form,

where $m_a$ denotes the mass of the boson particles related to the scalar field, and $f_a$ is regarded as the effective energy scale of the theory. This model, also referred to as PNGB in [160], has since been a recurrent feature in cosmological analyses [54, 222, 223, 224, 225, 226, 227, 228, 229]. The Klein-Gordon equation derived from the action in Eq. (31), in a standard cosmological setting, is

with $H$ being the Hubble parameter. We then perform the following polar transformation on the field variables $(\varphi,\dot{\varphi})$ (see [230, 231] and references therein),

Here, $M_{\mathrm{Pl}}$ is the reduced Planck mass, $\Omega_\varphi$ is the density parameter of the scalar field, and $\theta$ is an angular variable directly related to the scalar field equation of state $w_\varphi$,

Equation (D5) directly shows that, for quintessence models, the equation of state only varies in the range $-1\leq w_\varphi\leq1$.

As a result, the Klein-Gordon equation of motion can be rewritten as a dynamical system in the new variables, namely,

where a prime denotes a derivative with respect to $N=\ln a$, and $w_{\mathrm{tot}}$ is the total equation of state of the matter budget of the universe. The boson mass $m_a$ appears implicitly in the definition of the new dynamical variable $y\equiv2m_a/H$, which directly measures the ratio of the boson mass to the Hubble parameter. Likewise, the effective energy scale $f_a$ appears implicitly in the new parameter $\alpha=3/(f_a/M_{\mathrm{Pl}})^2$ and adopt following prior $m_a\in\mathcal{U}[0.1,10]$ and $\alpha\in\mathcal{U}[10^{-6},100]$ to sample mass and effective scale.

To solve the dynamical system Eq. (D6), we follow the prescription in [231]. The initial value of $y$ is determined by choosing a value of the boson mass $m_a$, that is, $y_i=2m_a/H_i$, while the initial angular variable is given by $\theta_i=(1/5)y_i$, which is an attractor solution at early times. As we will assume that $m_a\sim H_0\ll H_i$, the initial value of the equation of state is $w_{\varphi_i}\simeq-1$, which means that the field $\varphi$ starts its evolution close to a slow-roll regime. Finally, once a value of $\alpha$ is chosen, the initial value $\Omega_{\varphi_i}$ is adjusted using a numerical shooting routine inside the Boltzmann solver class until the desired value of $\Omega_\varphi$ at the present time is obtained. Using the polar transformations (D4b), the physical parameters $(m_a, f_a)$, and the initial values $(H_i, \Omega_{\varphi i})$, one can calculate the initial values of the original field variables $(\varphi_i, \dot{\varphi}_i)$ corresponding to the attractor solution at early times.

Validation on Mocks

Non-parametric approaches offer significant advantages by minimizing assumptions about the underlying physical properties of dark energy, allowing for a more flexible and unbiased reconstruction of cosmic expansion and structure growth. However, their power comes with the challenge of ensuring robustness, as the lack of an explicit model can introduce degeneracies and reconstruction artifacts. Careful validation using mock datasets is essential to assess the reliability of these techniques and identify potential biases that may arise due to the methodology itself, ensuring that the inferred constraints are driven by true cosmological signals and not systematic effects.

In this appendix, we assess the robustness of the non-parametric methodologies used in the main text with a set of two mock datasets, each generated from a distinct dark energy model and comprising a combination of simulated DESI, PantheonPlus SNe Ia, and CMB data. The first synthetic dataset assumes a $\Lambda$CDM cosmology $(w_0=-1,w_a=0)$, with the cosmological parameters set to their bestfit Planck values.

The second one is a $w_0w_a$CDM realization, with $w_0=-0.75$ and $w_a=-0.85$. We refer to these as the $\Lambda$CDM mock and the $w_0w_a$CDM mock, respectively. While admittedly not an exhaustive sample of possible dark energy models, these different mocks span a range of physical behaviors that allow us to quantify the statistical uncertainties and possible biases introduced by the methods themselves.

Figure 14 shows the binned $w(z)$ reconstruction in the case of 3 uniform bins between $z = 0$ and $z = 2.1$. The median values with 68% and 95% CL for the two mock datasets are shown. The true $w(z)$ for both mocks, plotted with a black dashed line, lies mostly within the $1\sigma$ contours, sometimes straying into the $2\sigma$ range near the edges of each bin where, as expected, the three uniform bin scheme is not necessarily flexible enough to capture the precise behavior. The deviation from $\Lambda$CDM (shown in gray dashed line) is correctly detected at more than $2\sigma$ in the two lowest redshift bins in the case of the $w_0w_a$CDM Mock. Figure 15 shows that Gaussian process (GP) reconstructions of $w(z)$ for the two mock datasets. The figure includes the median values as well as 68% and 95% CL for both of the mock datasets. The true expected function lies well within the $1\sigma$ contours in each case.

In the interest of completeness, we also test an extreme case, falling in the thawing class of models, by preparing a third mock following Eq. (15) with $w_0 = -0.5$ and $p = 20$. As shown in Figure 16, this demonstrates the possibility of limitations in the non-parametric implementation used here to pick up some extreme behaviors. We simply note that our implementation does have some limitations, and further improvements, such as investigating the impact of priors on hyperparameters and the choice of kernel, are left to future work.

Aside from the extreme case, the two approaches investigated here are seen to recover the simulated mock data well, without any significant bias detected in the mocks tested. Though not shown here, the same tests are performed for direct $f_{\rm DE}(z)$ reconstruction, with comparable analogous results. These non-parametric implementations are then applied without modification to the real data for the actual analysis.

Comparison with DESI DR1

Figure 17 compares the results of a Chebyshev expansion of $w(z)$ and $f_{\rm DE}(z)$ using DR1 BAO vs DR2 BAO data, in combination with $(\theta_\ast,\omega_b,\omega_{bc})_{\rm CMB}$ and Union3 measurements. Figure 18 presents a similar comparison for the GP reconstruction of $w(z)$ using DESI in combination with CMB and Union3. It is seen that while the main trends in the DE remain unchanged between the two data releases, the uncertainties in the reconstructions have significantly decreased with DR2.