Assessing observational constraints on dark energy

Observational constraints on time-varying dark energy ({\it e.g.}, quintessence) are commonly presented on a $w_0$-$w_a$ plot that assumes the equation of state of dark energy strictly satisfies $w(z)= w_0+ w_a z/(1+z)$ as a function of the redshift $z$. Recent observations favor a sector of the $w_0$-$w_a$ plane in which $w_0>-1$ and $w_0+w_a<-1$, suggesting that the equation of state underwent a transition from violating the null energy condition (NEC) at large $z$ to obeying it at small $z$. In this paper, we demonstrate that this impression is misleading by showing that simple quintessence models satisfying the NEC for all $z$ predict an observational preference for the same sector. We also find that quintessence models that best fit observational data can predict a value for the dark energy equation of state at present that is significantly different from the best-fit value of $w_0$ obtained assuming the parameterization above. In addition, the analysis reveals an approximate degeneracy of the $w_0$-$w_a$ parameterization that explains the eccentricity and orientation of the likelihood contours presented in recent observational studies.

One must be careful when portraying observational results in this way, because cosmological observables depend directly on the Hubble parameter H(z) and its integrals but only indirectly on the dark energy equation of state w(z).In particular, the functional form of w(z) can differ significantly from Eq. (1) in different models of quintessence, while still predicting very similar cosmological observables [12,13].In order to compare the predictions of different quintessence models to one another and to observations using a w 0 -w a plot, one can assign to each model the ordered pair (w 0 , w a ) that predicts an H(z) most similar to that of the quintessence model.This is similar to the approach suggested by de Putter and Linder [13], who showed that even just by assigning w a in this way, it is possible to match H(z) to ∼ 0.1% or better for some models of quintessence.
In this paper, we use this mapping technique to predict where observational preferences should fall on the w 0 -w a plane for a range of quintessence models driven by a scalar field φ with canonical kinetic energy density and potential V (φ).The models we consider, sometimes referred to as "thawing dark energy" [14], have the property that w(z) → −1 at early times (large z) because Hubble friction during the matter-and radiation-dominated eras is large enough to freeze the scalar field.At late times (small z), the Hubble friction becomes negligible, the scalar field accelerates down the potential, and w(z) increases as z decreases.We show that, when matching H(z)'s, these models are mapped onto the same sector of the w 0 -w a plot as observations currently prefer, even though they do not violate the NEC.As a result, observations favoring a region with w 0 + w a < −1 do not imply that dark energy must be NEC-violating at any redshift.
We detail our mapping protocol in Sec. 2 and discuss the uncertainties in this procedure (due partially to an approximate degeneracy in the w 0 -w a plane) in Sec. 3. The quintessence models we will test are presented in Sec. 4, and the results of their mapping onto the w 0 -w a plane are given in Sec. 5. Finally, in Sec.6, we summarize our findings and discuss the implications for interpreting observational likelihood contours on a w 0 -w a plane.In particular, we point out how the observations bear on issues such as NEC violation, consistency with supergravity and string theory, and on whether accelerated expansion continues forever or terminates and transitions to contraction.

Methods
As noted above, our procedure relies on determining which combination (w 0 , w a ) is the "best fit" to a given quintessence model of dark energy.Rather than choosing the fit that best mimics the quintessence field's equation of state w Q (z), we choose the fit that best matches the evolution of the Hubble parameter H Q (z) in the quintessence model.This allows the fit to most nearly reproduce key late-time cosmological observables, such as the Hubble distance the angular diameter distance and the luminosity distance Notably, when optimizing to match H(z), the integrated quantities D M (z) and D L (z) are found to match with similar accuracy, as we will show for a sample model in Sec. 5. Quantitatively, we define the "best fit" combination (w 0 , w a ) as one that minimizes the error Here, H Q (z) is the evolution of the Hubble parameter for a given quintessence model, while H fit (z) is the evolution for an artificial (w 0 , w a ) model.The maximization in Eq. ( 5) is performed over the finite interval z < 4 to roughly match the redshifts probed by BAO and SNe Ia measurements.These are also the redshifts during which the time-variation of dark energy is most relevant; at higher redshifts, the universe is strongly matter-dominated, and H(z) can be computed directly from the present-day matter density ρ m,0 , independently of the nature or behavior of dark energy.Note that there is some flexibility to the above definition of error; for example, the approach in Ref. [13] matches D M (z) over all z > 0 rather than H(z) over z < 4. In Sec. 3, we will discuss how this flexibility translates to an uncertainty in our best-fit results.
For some models of quintessence, H Q (z) can be easily calculated from analytically parameterized expressions for w Q (z) (see, e.g., Refs.[14,15]).However, these expressions are approximate and only valid in the regime where 1 + w Q (z) ≪ 1, which limits the range of models that can be studied.In this work, in order to find H Q (z) for a quintessence model driven by a canonical scalar field φ with potential V (φ), we switch variables from z to N ≡ − ln(1+z) and numerically solve the full equations of motion Here, V ,φ ≡ dV (φ)/dφ, w Q (N ) is the equation of state of the scalar field, Ω Q (N ) = 1 − Ω m (N ) is the fraction of total energy density attributable to the scalar field, and primes denote derivatives with respect to N .Note that these equations assume a spatially flat universe.We set our initial conditions deep in the matter-dominated past, with Ω Q = 10 −6 .The initial field value will depend on the particular model, but the field velocity φ ′ (N ) can be initially set to zero without loss of generality, as it will evolve toward an attractor trajectory while Ω Q ≪ 1.We end the simulation upon reaching a pre-selected fiducial value of Ω Q (0), which we take to be an input to the mapping procedure.Note that throughout this work, we will measure H Q (N ) in units of H Q (0), and in the examples we provide, we will assume that Ω Q (0) = 0.7.
Next, we can generate a candidate fit H fit (N ) by solving the system Here, we do not assume that the fitted fractional dark energy density Ω fit (0) equals Ω Q (0), nor do we assume that H fit (0) equals H Q (0).Instead, we treat these two initial values as free parameters, in addition to the combination (w 0 , w a ) specifying the equation of state w fit (N ).We do check, however, that the best-fit values for Ω m (0) ≡ 1−Ω fit (0), H fit (0), and the combination Ω m (0)H fit (0) 2 are within ∼ 1% of their values in the quintessence model.This ensures that the fit is compatible with CMB constraints in addition to measurements of the low-redshift universe.The final step in our procedure is to scan over the four free parameters and identify the best-fit combination of {w 0 , w a , Ω fit (0), H fit (0)} with the smallest error as defined in Eq. ( 5).We note that fitting all four parameters, rather than just w a as in the approach of de Putter and Linder, is necessary in order to properly identify the best fits for a wider range of quintessence models than those considered in Ref. [13].After projecting out the best-fit values of H fit (0) and Ω fit (0), we are left with a single point on the w 0 -w a plane that best represents the quintessence model we started with.More generally, this protocol can be used to map a one-parameter family of quintessence models onto a best-fit curve on the w 0 -w a plane, or a multi-parameter family of models onto a best-fit region.

Degeneracy and uncertainty
In our fitting process, it is instructive to identify a region of "acceptable" combinations of (w 0 , w a ) in addition to the central or "best-fit" combination.There is an approximate degeneracy between w 0 and w a when matching cosmological observables, causing this acceptable region to be highly eccentric.In particular, for any best-fit combination (w 0 , w a ) that matches a quintessence model's H(z) with reasonable accuracy, there will be a large set of combinations satisfying ∆w a /∆w 0 ≈ −5 that match the model with similar accuracy.The slope of this degeneracy can vary by ∼ 10% depending on the model being fitted, but its qualitative appearance is clear and distinct.
To demonstrate that this effect is fundamental to the (w 0 , w a ) parameterization and not dependent on any particular choice of quintessence model, Fig. 1 illustrates the degeneracy when fitting (w 0 , w a ) combinations to a fiducial model strictly obeying the parameterized equation of state in Eq. ( 1) with (w 0 , w a ) = (−0.827,−0.75); these are the central values of the DESI+CMB+PantheonPlus constraints [8].The darkest regions in the figure correspond to the (w 0 , w a ) combinations that most accurately match the fiducial model's H(z), and the dashed line with slope ∆w a /∆w 0 = −4.8indicates the axis of degeneracy.One can see that many combinations aligned along this axis are accurate fits to the central model within E ≲ 0.3%.This degeneracy is not only relevant for our theoretical fits to quintessence models, but it also makes a prediction about the orientation of constraint contours in the w 0 -w a plane produced by observational analyses.In particular, if the observed values of H(z) over a broad range of redshifts are well fit by one combination of (w 0 , w a ), then we expect that it will also be reasonably fit by other combinations along the axis of degeneracy.The resulting eccentricity and orientation of the constraint contours are indeed apparent in, e.g., Refs.[8][9][10][11], with the slope most closely matching our prediction in the combined BAO+CMB+SNe Ia results from DESI [8].Uncertainties associated with the fitting protocol.A line of best-fit (w 0 , wa) combinations for a hypothetical one-parameter quintessence model is shown in solid black.Each point on the line has an associated "acceptable-fit" region of variable size, depicted here qualitatively as blue elliptical contours satisfying E ≤ Emax for some Emax.The contours are generically largest near the ΛCDM limit, and they blend together to form a tapered acceptable-fit ribbon (shaded pink) around the best-fit line.A second type of uncertainty, representing the flexibility in the definition of error when determining the best-fit line itself, is illustrated qualitatively with two alternative best-fit lines (dashed) converging in the ΛCDM limit.
Unlike the orientation of the degeneracy, the size of the degenerate region in theoretical (w 0 , w a ) fits is model-dependent.In the example shown in Fig. 1, the best-fit error is E = 0 by construction, and there is a correspondingly large region of acceptable fits, depending on the precision E max with which H(z) can be measured.In general, a larger error E for the best-fit (w 0 , w a ) combination corresponds to a smaller region of acceptable fit.
When plotting a best-fit curve on the (w 0 , w a ) plane corresponding to a one-parameter family of quintessence models, the acceptable-fit regions of each point merge together to form an acceptable-fit ribbon, shown in Fig. 2 as the pink shaded region.The ribbon is widest at the ΛCDM limit (w 0 = −1, w a = 0) of any quintessence model, provided such a limit exists, but the rate at which the ribbon tapers off is model-dependent.
Finally, a second and entirely independent type of uncertainty reflects our confidence (or lack thereof) in the location of the best-fit curve itself on the w 0 -w a plane.This second uncertainty stems from a fundamental ambiguity in the notion of "best-fit," i.e., in the definition of error (Eq.5).We have observed that choosing a different definition of error-say, a mean-square error instead of a maximum, or matching D M 's instead of H's-can affect the slope of best-fit lines by ≲ 10%, depending on the model.The lines pivot about the ΛCDM limit (assuming such a limit exists in the models being fitted), where all definitions of error agree that the cosmological observables are best fit by (w 0 , w a ) = (−1, 0).This uncertainty in the slope is depicted by the dashed lines in Fig. 2.

Models
In this work, we illustrate the methods outlined in the previous section using three classes of thawing quintessence models as examples: exponential potentials, hilltops, and plateaus.Each of these models is driven by a scalar field φ with canonical kinetic energy density rolling down a potential V (φ).The scalar field is initially held constant during the radiation-and matter-dominated epochs by Hubble friction, such that w → −1 at large redshift.Then, as the dark energy density becomes comparable to the matter density, the Hubble friction decreases, the field accelerates down the potential, and in turn w increases as z decreases (or "thaws" away from −1).This behavior corresponds to a negative value of w a in the (w 0 , w a ) parameterization, as appears to be preferred by recent observational constraints [8][9][10][11].
We parameterize the exponential potential as where φ is in units of the reduced Planck mass M pl = ℏc/(8πG) and the initial field value is set to zero, such that V 0 ∼ O(H 2 0 M 2 pl ).A theoretical motivation for studying this class of models is that scalar fields with exponential potentials are ubiquitous in supergravity, modified gravity, and superstring theories; see for example [16][17][18].Additionally, the exponential potential probes two interesting limits of quintessence models: one where V ,φ /V is constant and much smaller than M −1 pl ("slow-roll" thawing quintessence [14,19]), and one where V ,φ /V is constant but not small (i.e., when λ ≳ M −1 pl ).Hilltop potentials with large, negative V ,φφ provide a complementary probe into the regime where V ,φ /V is small compared to M −1 pl but not constant.Potentials of this type are good approximations to axion models (or pseudo-Nambu Goldstone bosons generally) provided the scalar field is initially frozen by Hubble friction near the top of its potential [15].We take a quadratic approximation and write the hilltop potential as where again φ is in units of M pl and V 0 ∼ O(H 2 0 M 2 pl ).In this work, we analyze a relatively flat hilltop with k 2 M 2 pl = 1 and a more concave hilltop with k 2 M 2 pl = 100.For each case, we will use a variety of initial field values to map out the set of possible best-fit (w 0 , w a ) combinations to the given potential.
The plateau potential we consider is something of a hybrid between the previous two examples, pairing a region of smooth exponential decay with a sharp, cliff-like drop: This is an example of a case where V ,φ /V cannot be assumed to be either small or constant.We choose an exponential cliff (rather than a power-law) to match the model of Ref. [20], where it was shown that the potential ( 15) can lead to a transition from accelerated expansion to a regime of slow contraction when V plateau becomes negative, as can occur in a cyclic cosmology.For this model, we are free to set the initial field value to φ = 0 without loss of generality.As before, φ is in units of M pl , V 0 ∼ O(H 2 0 M 2 pl ), and κ ≲ 1 is a constant.This model has three free parameters, κ, M , and m, which leads to a best-fit region on the w 0 -w a plane rather than a curve.Note that in any limit where the second term in the potential becomes negligible, we recover the exponential model with λ = 1/M .

Results
The best-fit curves and regions for the three classes of models introduced in Sec. 4 are depicted in Fig. 3.As mentioned in Sec. 2, this analysis assumes a fiducial value of Ω Q (0) = 0.7.
The orange curve corresponds to the exponential model.The limit λ → 0 maps onto the ΛCDM parameters w 0 = −1 and w a = 0.Each fit along this curve has error E ≲ 0.1%, which is the greatest level of accuracy among the models we tested.
The region between the exponential curve (orange) and the plateau boundary (red) corresponds to the three-parameter family of plateau models.The upper boundary (the exponential curve) corresponds to the limit κ ≪ 1 in which the cliff is negligibly small.The lower boundary corresponds to large κ and a substantial cliff that causes a sudden increase in the equation of state as z decreases.This sudden increases in w(z) causes the models to be mapped onto more negative values of w a .The errors near this boundary can be somewhat greater but still satisfy E ≲ 0.5% for the cases we tested.We note that the lower boundary of this region is approximate, and other types of plateau models may not mimic the behavior of the specific potential examined in this work.
The two remaining curves correspond to hilltop models with k 2 M 2 pl = 1 (green) and k 2 M 2 pl = 100 (blue).The respective errors satisfy E ≲ 0.2% and E ≲ 0.7% within the region shown in the plot.The closer the initial field value is to the hilltop, the closer the best-fit (w 0 , w a ) model is to ΛCDM, and the smaller is the best-fit error.
These results are overlaid with the observational constraints obtained by the DESI collaboration (blue and brown 2σ contours) when combining measurements of BAO, CMB, and SNe Ia [8].This allows the plot to be used for comparing quintessence models to each other and to the observations at the same time.For example, the likelihood contours in Fig. 3 appear to favor quintessence potentials with sharp drops, like plateau models with the steepest cliffs or hilltop models with k 2 M 2 pl ∼ 100.However, we make this point mainly for the purpose of illustrating the utility of w 0 -w a plots in general, and we would not suggest drawing any strong conclusions based on currently available data.
Notably, all the quintessence models we tested obey the NEC and have w(z) > −1 for all z, yet they are best fit by curves or regions on the w 0 -w a plane that satisfy w 0 + w a < −1.This boundary is marked by the dashed line in Fig. 3.We therefore conclude that a preference for this sector of the plane, which would naively imply NEC violation at large z according to equation (1), is actually compatible with NEC-satisfying models of quintessence.This result is a manifestation of the fact that even large differences in w(z) over a wide range of z ≫ 1 (where Ω Q (z) is negligibly small) have a negligible impact on cosmological observables such as H(z).[Bottom panel] The evolution of w Q (z) for the same hilltop model (solid curve), compared with the evolution of w fit (z) for the best-fit (w 0 , wa) combination (dashed curve) based on matching H(z)'s.
At smaller redshifts, H(z) is more sensitive to differences in w(z), but only through an integral.As a result, two models can have very similar H(z)'s even if their w(z)'s differ noticeably over a narrow range of redshifts near z = 0.In turn, this means that the best-fit values of w 0 obtained when matching H(z)'s need not be equal to-or even close to-the actual value w Q (0) for the quintessence model being considered.This is specifically the case for the hilltop and plateau models, in which w Q (z) is increasing rapidly as z → 0. We illustrate an example of this behavior for a hilltop model with k 2 M 2 pl = 100 in Fig. 4, juxtaposing the fits to H(z) and D M (z), which have sub-percent level errors, with a comparison of w(z)'s between the model and the fit, which differ by O(100%).In fact, the best-fit curves drawn in Fig. 3 include cases where w Q (0) > −1/3, in which case the universe today is no longer accelerating, while the best-fit (w 0 , w a ) combination still has w 0 ≤ −0.65.
As discussed in Sec. 3, there is some uncertainty in our best-fit results due to the flexibility in how we define the error E. This uncertainty is negligible for the exponential curve, but the plateau boundary and hilltop curves should be interpreted as having a ∼ 10% error bar on the value of w a at any given w 0 .Additionally, extending the plateau and hilltop curves beyond the region shown would produce increasingly large best-fit errors E, ultimately reaching a point where there exists no good fit to the quintessence model within the (w 0 , w a ) parameter space.In cases like this-including more general models of dark energy or modified gravity that are not well fit by any (w 0 , w a ) combinationthe safer approach would be to perform a separate analysis specific to the model in question, rather than trying to map it onto, and then constrain, the w 0 -w a parameter space.
Finally, we note that plots like Fig. 3 carry no information about finetuning of either parameters or initial conditions for a given model.For example, in a hilltop model, only a small range of initial field values is simultaneously sufficiently close to the peak of the hilltop to produce a noticeable period of accelerated expansion and sufficiently far from the hilltop to be distinguishable from ΛCDM.This fine-tuning must be judged and weighed separately.

Discussion
In this work, we modified and illustrated a mapping protocol that assigns a best-fit value of (w 0 , w a ) to any given model of quintessence based on matching the evolution of the Hubble parameter H(z) in a spatially flat universe.For the cases we tested, whose maps are shown in Fig. 3, we confirmed that the (w 0 , w a ) parameterization was able to match the quintessence model's H(z) with less than 0.7% error at all z.We note that for each model considered in this work, the Swampland conjectures [18,21,22] (thought to be required for a consistent theory of quantum gravity) are satisfied for a substantial range of parameters.The plateau and hilltop potentials, which are sharply decreasing and can become negative, also have a natural place in cyclic bouncing cosmology [23,24].
Though less precise than Bayesian model comparisons using Markov-Chain Monte Carlo simulations, this protocol provides a quick and useful way to visually compare different models of quintessence to each other and to observational data in a common two-dimensional parameter space.For example, if we take the constraints reported by DESI at face value, we see that plateau and hilltop models fare better than exponential models, with highly concave hilltop models being most compatible with the data.We reiterate, however, that this comparison test does not take into account any fine-tuning of the parameters or initial conditions of models that land within the observational constraint contours, and it is subject to the uncertainties discussed in Sec. 3.This mapping protocol not only allows us to assess models of quintessence against observational data, but it also sheds new light on how to interpret the observational likelihood contours themselves.First, we have shown that the eccentricity and orientation of contours generated from measurements across a broad range of redshifts is a generic feature of the (w 0 , w a ) parameterization.Second, we have found that for some models of thawing quintessence, the best-fit value of w 0 based on matching H(z)'s can differ significantly from the true value of w Q (0) predicted by the model.Finally, we have pointed out that the thawing quintessence models analyzed in this work, all of which obey the NEC (i.e., have w Q (z) ≥ −1), are mapped onto (w 0 , w a ) combinations that satisfy w 0 > −1 and w 0 + w a < −1.An observational preference for this sector, therefore, does not require the kinds of exotic field theories needed to enable a transition from NEC violation at large z to NEC compliance at small z (see, e.g., Ref. [25]).
This last finding has an important corollary: contrary to the suggestion in Ref. [26], we have shown that it is not just reasonable but crucially important for observational analyses to include combinations of (w 0 , w a ) satisfying w 0 + w a < −1 in their priors with high credence.Otherwise, these analyses would be inadvertently excluding families of simple, well-motivated models of thawing quintessence from consideration.

Figure 1 :
Figure 1: A mild degeneracy among (w 0 , wa) models is illustrated by computing the error E of each (w 0 , wa) fit when matching the H(z) of a fiducial (w 0 , wa) model with w 0 = −0.827,wa = −0.75.The orientation of the degeneracy is indicated by the dashed line with slope ∆wa/∆w 0 = −4.8.

Figure 3 :
Figure 3: A w 0 -wa plot showing the predictions for three types of canonical quintessence potentials (exponential, plateau, hilltop) that are NEC-satisfying.For plateau models, the predictions map onto the region between the exponential curve and the plateau boundary curve indicated.The plot also shows the best-fit contours from DESI BAO+CMB+(PantheonPlus or Union3) that prefer the sector w 0 + wa < −1 (below the dark dashed line), naively suggesting NEC violation at large z.

Figure 4 :
Figure 4: [Top panel] Fractional difference in H(z) and D M (z) between a hilltop quintessence model with k 2 = 100 and the best-fit (w 0 , wa) model.Note that |∆D L |/D L,Q = |∆D M |/D M,Q .For this model, as for all models considered in this work, the Hubble parameter H(z) and the integrated quantities D M (z) and D L (z) are matched with sub-percent errors at all z.[Bottom panel] The evolution of w Q (z) for the same hilltop model (solid curve), compared with the evolution of w fit (z) for the best-fit (w 0 , wa) combination (dashed curve) based on matching H(z)'s.