1 Introduction↩︎

Accelerated cosmic expansion driven by a cosmological constant (\(\Lambda\)) has been remarkably successful in explaining a plethora of cosmological observations [1]. This \(\Lambda\) corresponds to a constant dark energy density \(\rho_\Lambda=\Lambda/8\pi G\). However, recent baryon acoustic oscillations measurement by the DESI collaboration [2], [3], in combination with Type Ia supernova distances and CMB constraints, favor models where the effective dark energy density is decreasing at late times. Kinematically, there is a preference for decreasing the acceleration close to today, albeit, the statistical strength depends on the supernova sample and on the adopted modeling assumptions (e.g. [2], [4]–[6]).

Canonical quintessence, a scalar field rolling in a potential, provides one of the simplest controlled realizations of dynamical dark energy [7]–[10]. However, conventional slow-roll or thawing models do not fully capture the DESI-preferred evolution. Phenomenological fits often point toward phantom-divide crossing (e.g. [4]), but such behavior is difficult to realize in stable, canonical theories and may partly reflect the flexibility of the chosen parameterization rather than a robust physical signal [11], [12]. These considerations motivate a simple canonical mechanism that remains close to \(\Lambda\)CDM until very late times, but then naturally reduces the dark-energy density today.

In this Letter we propose a “simple" dark energy model to explain these recent observations: a very late-time oscillating quintessence field. The general form of the quintessence potential we have in mind is illustrated schematically in 1. The field remains Hubble-frozen, or slowly rolling in the flatter part of the potential for most of cosmic history. At late times, it enters a fast-rolling, oscillatory regime around the minimum.

When confronted with recent observational data, we find the above model improves the fit significantly relative to \(\Lambda\)CDM, while remaining competitive with phenomenological dark energy parameterizations with the same number of additional parameters. For the best-fit region, the quintessence field starts oscillating at a redshift \(z_{\rm osc}\approx 0.1\) and completes one oscillation before today. The improvement in fit has a simple origin. When the field is Hubble frozen for most of cosmic history, its energy density is nearly constant and the model remains close to \(\Lambda\)CDM. After the recent onset of oscillations, the dark energy density (and associated acceleration) begins to decrease naturally.

This type of late-oscillating quintessence was studied in Ref. [13], where along with the change in the expansion rate due to the field oscillations, many rapid field oscillations were also shown to lead to resonant amplification of quintessence field fluctuations (also see [14]). This led to a corresponding scale-dependent growth of gravitational potentials and an enhanced Integrated Sachs-Wolfe effect [13]. Similar sequences of slow evolution, rapid oscillation, and possible resonance occur in other cosmological epochs, including post-inflationary dynamics and oscillating early dark energy [15]–[18]. Related oscillatory dark energy ideas have also been explored in connection with the coincidence problem [19].

In the present paper, we focus on whether this canonical scalar field mechanism can produce the background dynamics preferred by observations. We also include the associated impact of quintessence perturbations, and find that they remain small in the best-fit region.

2 Model & Dynamics↩︎

We consider the late-time oscillating quintessence potential introduced in Ref. [13], specialized to \(\alpha=0\), \[V(\phi)=\mathcal{V}\frac{(\phi/M)^2}{1+(\phi/M)^2},\quad \mathcal{V}\equiv \frac{m^2 M^2}{2}.\] Here \(m\) sets the curvature of the potential \(V''(\phi)\approx m^2\) at its (quadratic) minimum, while \(M\) controls the transition between a shallow large-field plateau and the small field quadratic minimum. Note that \(\mathcal{V}\) sets the height of the potential at large field values, and \(m\) sets the frequency of small amplitude field oscillations around the minimum.

While we use the above explicit form of the potential in what follows, we want to emphasize that the results are more general. In Appendix C, we show that a broader class of scalar-field potentials with qualitative features similar to the one in Fig. 1, and with varying high energy physics motivations [20]–[22], work equally well for our purposes. These models include cases in which the potential has non-quadratic minima (\(\sim |\phi|^{2n>2}\)).

For addressing the hints of late-time evolving dark energy, the dynamics we would like are ones where \(\phi\) is essentially frozen during radiation and matter domination with \(V(\phi_{\rm ini})\sim m_{\rm pl}^2H_0^2\), and begins oscillating close to today with frequency \(m\gg H_0\). This can be achieved if our quintessence field is misaligned to \(\phi_{\rm ini}\gg M\) at the end of inflation and \(\mathcal{V}\sim m_{\rm pl}^2H_0^2\). For the field to start entering the minimum region (\(\phi/M\sim 1\)) close to today, we would need the initial displacement of the field \(\phi_{\rm ini}/M\sim \sqrt{m_{\rm pl}/M}\). Overall, we need the high energy theory to provide \(m\) and \(M\) such that \(\mathcal{V}=m^2M^2/2\sim m_{\rm pl}^2 H_0^2\), with \(m\gg H_0\), while the initial misalignment must satisfy \(\phi_{\rm ini}/M\sim \sqrt{m/H_0}\).

We create a collection of model-realizations by specifying \(a_{\rm osc}\) and \(M\), drawn from a uniform distribution over \(a_{\rm osc}\) and \(\log_{10}(M/m_{\rm pl})\).¹ Here, \(a_{\rm osc}\) is the scale factor at which \(\phi\) first crosses zero, approximately characterizing the onset of oscillations. For each \(M\) and \(a_{\rm osc}\), we adjust \(\phi_{\rm ini}\) and \(m/H_0\) simultaneously so that spatial flatness (correct energy density at early times) is maintained while the desired value of \(a_{\rm osc}\) is reproduced. Field realizations with small \(a_{\rm osc}\) undergo extremely high-frequency oscillations that make the numerical computation challenging and slow. Moreover, the approximation adopted below for the field perturbations may not be sufficiently accurate in this regime. To avoid these problems, we impose a lower bound \(a_{\rm osc}=0.7\).

In general the quintessence field can be written as \(\phi(a,\boldsymbol{x})=\phi(a)+\delta\phi(a,\boldsymbol{x})\). We modified CLASS [23] to solve the background and (linearized) perturbation evolution. As the homogeneous \(\phi(a)\) undergoes multiple oscillations, parametric resonance leads to a scale-dependent, exponential growth in quintessence perturbations \(\delta\phi(a,\boldsymbol{x})\). This growth leads to breakdown of the linearized evolution of perturbations when \(\delta\phi(a,\boldsymbol{x})\sim \phi(a)\). Once this point is reached, the exponential growth is curtailed and nonlinear evolution with significant mode-coupling commences. Following Ref. [13], we define \(a_{\rm nl}\), the scalefactor when nonlinear evolution commences, as the moment when the envelope of \(\max_k\sqrt{k^3P_{\delta\phi}(a,k)/(2\pi^2)}\) exceeds the envelope of the homogeneous amplitude \(|\phi(a)|\). Here \(k^3P_{\delta\phi}(a,k)\sim \langle \delta\phi^2(a,\boldsymbol{x})\rangle_{L\sim k^{-1}}\) characterizes the typical variance of quintessence perturbations on the scale \(k^{-1}\).

From the Poisson equation, we can see that during the resonant growth phase of \(\delta\phi\), the gravitational potential \(\Psi\) also grows rapidly. For \(a>a_{\rm nl}\) we freeze the gravitational potential according to \(\Psi(a>a_{\rm nl},\boldsymbol{x})=\Psi(a_{\rm nl},\boldsymbol{x})\). In practice this is implemented with two runs of CLASS: the first determines \(a_{\rm nl}\) and the second applies the freeze-out prescription. As we show below, the best-fit models do not in fact enter the resonant growth regime by today (i.e., \(a_{\rm nl}>1\)).

3 Data↩︎

The low-redshift data set consists of DESI DR2 BAO measurements [3], which constrain the standard-ruler combinations \(D_M(z)/r_d\), \(D_H(z)/r_d\), and \(D_V(z)/r_d\), where \(D_M\) is the transverse comoving distance, \(D_H \equiv c/H(z)\) is the Hubble distance, \(D_V\) is the volume-averaged distance, and \(r_d\) is the sound horizon at the baryon drag epoch. It also includes the Union3 compilation [24] of 2087 Type Ia supernovae, which constrains the distance moduli \(\mu(z)\). We incorporate early universe information in two ways. First, we use the Planck PR4 CamSpec TT,TE,EE spectra [25] with \(\ell<1000\), supplemented by the PR3 low-\(\ell\) TT and EE likelihoods [26]. The high-\(\ell\) spectra are excluded because they are sensitive to lensing effect, and standard prescriptions for nonlinear corrections to the matter power spectrum are not calibrated for quintessence. As a complementary approach, we also incorporate CMB information by replacing the spectra with the compressed PR4 likelihood of Ref. [27], a multivariate Gaussian constraint on \((\theta_*,\Omega_bh^2,\Omega_{bc}h^2)\), where \(\theta_*\) is the angular size of the sound horizon at photon decoupling, \(\Omega_bh^2\) is the physical baryon density, and \(\Omega_{bc}h^2\equiv(\Omega_b+\Omega_c)h^2\) is the physical density of baryons + CDM. This prior preserves the CMB calibration of the sound horizon and physical matter densities while minimizing dependence on late-time lensing and nonlinear structure formation.

We sample the quintessence parameters along with \(\{\Omega_bh^2, \Omega_{c}h^2, H_0, A_s, n_s, \tau_\text{reio}\}\)² using Cobaya [28] until the Gelman-Rubin convergence statistic satisfies \(R-1<0.01\), and determine best-fit points with the optimizer BOBYQA [29], [30]. We also used prospect [31] to compute the likelihood profile.

4 Results and discussion↩︎

We perform a frequentist likelihood profile analysis for DESI + Union3 + compressed CMB first, in which we fix \(a_{\rm osc}\) to selected values and optimize the remaining parameters to obtain the corresponding best-fit \(\chi^2\). The result is shown in 2. It yields a best-fit with 68% confidence interval of \(a_{\rm osc}=0.916^{+0.018}_{-0.062}\). This best-fit value of \(a_{\rm osc}<1\) implies that the field entered the oscillatory regime in the past. We also note a mild reduction in \(\chi^2\) for \(a_{\rm osc}\) slightly above unity, corresponding to a thawing field that is about to enter the oscillatory regime. If the onset of oscillations is pushed further into the future, however, the quintessence dynamics approach those of \(\Lambda\)CDM and are correspondingly less favored.

3 shows the corresponding Bayesian marginalized posterior distribution of the field parameters in the \((\log_{10}(M/m_{\rm pl}),\log_{10}(m/H_0))\) plane for the two data combinations. The allowed region follows a narrow diagonal, reflecting the fact that reproducing the observed dark energy scale on the plateau approximately fixes the product \(mM\) through \(\mathcal{V}=m^2M^2/2\), which should be of the same order as the critical density today. The compressed CMB best-fit point lies near \(\log_{10}(M/m_{\rm pl})\simeq -2.2\) and \(\log_{10}(m/H_0)\simeq 2.5\).³ The region around the best-fit point lies on the \(a_{\rm nl}>1\) side of the boundary shown in the figure, where the field has not yet entered the nonlinear regime by today. Nevertheless, regions with \(a_{\rm nl}<1\) are still allowed.

4 shows the late-time evolution of the background quantities at the best-fit point for DESI + Union3 + compressed CMB data. Each time the scalar field crosses zero, the kinetic energy contribution temporarily dominates, causing the effective equation-of-state parameter \(w_{\rm DE}\) to rise rapidly to \(1\) before returning to \(-1\). As a result, the dark energy density \(\rho_{\rm DE}\) decays in a step-like manner. Consequently, the expansion rate \(H(z)\) is relatively reduced at very low redshift for \(a\gtrsim a_{\rm osc}\). For low-redshift observations, DESI and Union3 prefer slightly larger distances at very low redshift (\(z\lesssim 0.3\)) than at somewhat higher redshift, see 5. The reduction of \(H(z)\) induced by the scalar field oscillation at very low redshift naturally produces this behavior.

For DESI + Union3 + CMB (\(\ell<1000\)), we find that our model improves the fit relative to \(\Lambda\)CDM by \(\Delta\chi^2=-9.32\), and \(\Delta\)DIC = \(-4.49\)⁴. For DESI + Union3 + compressed CMB, we find an improvement of \(\Delta\chi^2=-9.19\), and \(\Delta\)DIC = \(-3.31\).

We also test the Padé-\(w\) phenomenological dark energy parameterization on the same data sets, \[w_{\rm DE}(a)=\frac{2\epsilon_0}{2+\eta_0(a^{-3}-1)}-1,\] which also has two additional free parameters, \(\eta_0\) and \(\epsilon_0\), always satisfies \(w_{\rm DE}(a)>-1\), and allows a rapid recent growth of \(w_{\rm DE}\) (blue dot-dashed line in 4). For this parameterization, the improvements relative to \(\Lambda\)CDM for these two datasets are \(\Delta\chi^2=-7.74\) and \(\Delta\chi^2=-6.64\), respectively. Our concrete quintessence model performs better than Padé-\(w\) parameterization.

For the \(w_0w_a\) (Chevallier–Polarski–Linder, CPL) parameterization [35], [36], which allows phantom crossing, \[w_{\rm DE}(a)=w_0+w_a(1-a),\] the improvements relative to \(\Lambda\)CDM are \(\Delta\chi^2=-11.28\) and \(\Delta\chi^2=-18.33\), respectively. Our model performs similarly to the \(w_0w_a\) parameterization on the DESI + Union3 + CMB (\(\ell<1000\)) data set, but is less favored than \(w_0w_a\) for the DESI + Union3 + compressed CMB data set which has more information from high-\(\ell\) part of CMB. ⁵ This is mainly because phantom crossing allows the dark energy contribution to be reduced during the matter-dominated era. As noted in the Introduction, however, current observations may be pointing more broadly to rapidly evolving low-redshift dark energy, rather than uniquely to phantom-divide crossing.

Before moving on to our conclusions, we reiterate that for our best-fit model, along with \(m\sim 10^2H_0\) and \(M\sim 10^{-2}m_{\rm pl}\), we find \(\phi_{\rm ini}\sim 10M\sim 0.1m_{\rm pl}\). The field excursion is therefore sub-Planckian. A model building question concerns the initial field displacement ⁶. If the potential remains flat significantly beyond \(10M\), the initial field misalignment would appear to be fine-tuned. A periodic potential with shallow regions between neighboring minima and a period of order \(10M\) could ameliorate such tuning. Related discussions initial-misalignment tuning for cosine potentials can also be found in, for example, [32], [37].

5 Conclusions↩︎

We have shown that late-time oscillating quintessence provides a simple and physically motivated realization of the low-redshift dark energy dynamics suggested by current distance data. In this scenario the scalar field is Hubble frozen, or slowly rolling, for most of cosmic history, so the model remains close to \(\Lambda\)CDM until very late times. Once the field reaches the steeper region near the minimum of its potential, it begins to oscillate, causing the dark energy density to decrease in a step-like manner. This behavior naturally reduces the expansion rate at very low redshift and thereby produces the pattern of distances preferred by the DESI BAO and Type Ia supernova data.

Using DESI DR2 BAO, Union3 supernovae, and two complementary treatments of CMB information, we find a preference for the onset of oscillations at \(a_{\rm osc}\simeq 0.9\), corresponding to \(z\simeq 0.1\). The improvement relative to \(\Lambda\)CDM is \(\Delta\chi^2\simeq -9\). The same model also gives a better fit than the Padé-\(w\) phenomenological parameterization with the same number of additional parameters, while remaining competitive with the CPL form on the data set using CMB spectra with \(\ell<1000\). Unlike purely phenomenological descriptions, however, the oscillating quintessence scenario arises from a canonical scalar field with a stable minimum and does not require crossing the phantom divide. More broadly, scalar fields that end in an oscillatory phase arise in several cosmological epochs, from the end of inflation and pre-recombination early dark energy and dark matter to late-time dark energy. Despite the difference in cosmological epochs and energy scales, the similar phenomenology provides an intriguing unifying theme.

In this work, the preference for the model arises primarily from its impact on the background expansion history, especially at very low redshift; more precise future Type Ia supernova measurements can therefore provide stronger constraints. At the perturbation level, we have treated the nonlinear evolution of the scalar field approximately. This is adequate in the region near the best fit, since the field has not yet entered the nonlinear regime there. For models in which the oscillations begin earlier, however, a more accurate nonlinear treatment in the late universe may be required (e.g. [38]). Such improvements could help to further distinguish these scenarios using observables such as the ISW effect and lensing.

6 Acknowledgments↩︎

We thank Christopher Cook and Francis-Yan Cyr-Racine for early work on this project, and Andrew Long, Mudit Jain and Scott Dodelson for helpful discussions. M.A.is supported by a U.S.Department of Energy (DOE) Grant DE-DE-SC0010103. MA also acknowledges that part of this work was done at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452.

7 Appendices↩︎

8 A. Full posterior distribution↩︎

This appendix summarizes the priors and posterior constraints used in the main analysis. The priors adopted in our analysis are summarized in Table 1. The full marginalized constraints are given in Table 2, and the corresponding one- and two-dimensional posterior distributions are shown in Fig. 6.

The two CMB treatments lead to broadly consistent posteriors for both the cosmological parameters and the scalar-field parameters. In both cases the marginalized posterior for \(a_{\rm osc}\) peaks at the region with oscillation in the past \(a_{\rm osc} < 1\).

9 B. Alternative parameterization↩︎

In the main analysis we use a flat prior on \(a_{\rm osc}\), which directly targets the epoch at which the field first crosses the minimum and begins oscillating. As a robustness check, we also consider a parameterization in terms of the underlying field parameters. Specifically, we impose flat priors on \(\log_{10}(M/m_{\rm Pl})\in[-6,3]\) and \(\log_{10}(\Lambda/\sqrt{H_0m_{\rm Pl}})\in[0,0.3]\), where \(\Lambda^4 \equiv \mathcal{V}\). This range covers solutions whose oscillations begin earlier than the lower bound \(a_{\rm osc}=0.7\) used in the main analysis, as well as thawing solutions whose first zero crossing would occur far in the future. Although we modified CLASS to estimate \(a_{\rm osc}\) for future oscillation solutions, many points in this branch do not cross zero within the numerical integration range, implying \(a_{\rm osc}\gtrsim 100\) for these samples. In the corresponding chains, samples with \(a_{\rm osc}>1\) account for \(21.68\%\) of the posterior sample weight, so the posterior still predominantly favors models in which the onset of oscillations has already occurred.

The one-dimensional marginalized posterior for the \(a_{\rm osc}<1\) branch is shown in 7. The alternative parameterization gives a posterior shape similar to that obtained with the flat prior on \(a_{\rm osc}\). This comparison indicates that the preference for late-time oscillations is not an artifact of using \(a_{\rm osc}\) as a sampling parameter. We also emphasize that the likelihood profile in 2 is a frequentist result obtained by optimizing the remaining parameters at fixed \(a_{\rm osc}\), and is therefore independent of the Bayesian prior used to sample the posterior.

10 C. Alternative scalar-field potentials↩︎

In this appendix we consider alternative scalar-field potentials that can produce similar late-time oscillatory behavior. As a representative example, we consider the T-model potential [39], \[V(\phi)=\Lambda^4 \tanh^{2n}\left(\frac{|\phi|}{M}\right).\] This potential has the same qualitative structure as the model used in the main text: it approaches a plateau at large field values, \[V(|\phi|\gg M)\simeq \Lambda^4,\] while near the minimum it reduces to a power-law form, \[V(|\phi|\ll M)\simeq \Lambda^4 \left|\frac{\phi}{M}\right|^{2n}.\] Oscillatory dynamics in this class of potentials have also been discussed in the context of post-inflationary evolution [16].

We study the cases \(n=1\) and \(n=2\), corresponding respectively to quadratic and quartic behavior near the minimum. The field evolution in these models is numerically more demanding than in the fiducial potential, because the late-time oscillations are more rapid. We therefore run the MCMC chains until they reach \(R-1<0.05\) and search the best-fit points based on the MCMC results. For \(n=1\), the best-fit point has \(a_{\rm osc}=0.915\), while for \(n=2\) it has \(a_{\rm osc}=0.911\). Both values are below unity, indicating that the field begins oscillating in the past. Relative to the best-fit \(\Lambda\)CDM model for the same data set, the corresponding improvements are \(\Delta\chi^2=-8.2\) and \(\Delta\chi^2=-7.7\), respectively. These results indicate that the preference for recent oscillatory dynamics is not a special feature of the particular potential adopted in the main analysis.

11 D. Results with alternative SNe Ia dataset↩︎

As a further robustness check, we compare the marginalized posterior distribution of \(a_{\rm osc}\) after replacing Union3 with alternative Type Ia supernova data sets, using either Pantheon+ [40], [41] or the new DES-Dovekie likelihood [42]. In each case we keep the DESI BAO data and compressed CMB prior fixed, so that the comparison isolates the impact of the supernova sample.

The comparison is shown in 8. For DES-Dovekie, the best-fit point is \(a_{\rm osc}=0.961\), while for Pantheon+ it is \(a_{\rm osc}=0.996\). Both values are below unity, although to different degrees, indicating that the field still prefers to have crossed the minimum before today. The preference is weakest for Pantheon+, where the best-fit onset is very close to the present epoch. This is consistent with the broader background-level behavior of these supernova samples: even without invoking oscillations, Pantheon+ is relatively less favorable to dynamical dark energy than Union3 or DES-Dovekie.

References↩︎