Underdetermination of dark energy

There is compelling evidence that the Universe is undergoing a late phase of accelerated expansion. One of the simplest explanations for this behaviour is the presence of dark energy. A plethora of microphysical models for dark energy have been proposed. The hope is that, with the ever increasing precision of cosmological surveys, it will be possible to precisely pin down the model. We show that this is unlikely and that, at best, we will have a phenomenological description for the microphysics of dark energy. Furthermore, we argue that the current phenomenological prescriptions are ill-equipped for shedding light on the fundamental theory of dark energy.


I. INTRODUCTION
Evidence for the accelerating expansion of the universe has been steadily accumulating since its discovery in 1998 [1][2][3][4][5][6][7].While these observations are still consistent with a cosmological constant, there are other proposals that seek to account for this accelerating expansion by postulating the existence of dark energy [8].Dark energy typically refers to an exotic form of matter, represented by additional fundamental field(s) that dynamically evolve over time and may (or may not) modify gravity on cosmological scales [9], depending on the exact nature of the field.
Characterizing this accelerating expansion usually begins by defining the equation of state.
where P is the pressure and ρ is the energy density of whatever is sourcing the accelerating expansion.If expansion is driven by a cosmological constant, the equation of state is locked in at the value w = −1.If, on the other hand, expansion is driven by a dynamical field, it will evolve over time and the equation of state can then be described as a function of a, the scale factor of the universe, so that w = w(a).
Ultimately, one would like to be able to extract the detailed time dependence of the equation of state in the form of w(a) or w(z) where z is the redshift, 1 + z ≡ a −1 and, indeed, there have been attempts at doing so [10,11].In practice there is limited information that one can extract about the evolution of w(a) and, as a result, this has led to the wide-spread adoption of a very natural and popular parameterization known as the Chevallier-Polarski-Linder (CPL) parameterization [12,13], given by: w(a) = w 0 + w a (1 − a). ( This parameterization approximates the dark energy equation of state close to today (a ∼ 1); where w 0 gives the value of the equation of state now and w a characterizes the temporal evolution of dark energy.This allows for a very useful description of dark energy models in terms of the parameters (w 0 , w a ).Fig. (1) shows the current constraints on these parameters.A cosmological constant, Λ, is of course given by w a = 0 and w 0 = −1, while various dynamically driven dark energy possibilities occupy the rest of the parameter space.
One interesting thing to note is that, while a cosmological constant is still favoured, there is a tremendous amount of open parameter space and future work that must be done in order to measure (w 0 , w a ).This invites the following question: what can pinning down the values of (w 0 , w a ) teach us about the microphysical nature of dark energy?As others have noted [14], constraining down to w a ∼ 0 and w 0 ∼ −1 strongly points towards a cosmological constant, but can never fully eliminate dynamically driven dark energy as a possibility because it may simply not yet have entered a stage of temporal evolution that is observable to us.There is of course another possibility.What if future surveys indicate that w 0 ̸ = −1 and that there is significant temporal evolution captured in w a ?
Clearly, observing w a ̸ = 0 with statistical significance would rule out the basic cosmological constant scenario and point towards the existence of dark energy driven expansion.However, would such observations allow us to definitively say anything more specific about dark energy, other than inferring its dynamical nature?Most obviously, we would like to gain an understanding of the fundamental microphysics responsible of dark energy [14][15][16][17][18][19][20][21][22][23][24][25].This presumably (in keeping with the field theoretic paradigm of modern physics) would be captured in a Lagrangian expression; one that includes structural information concerning the relevant types of dynamical field(s), the couplings to gravity and other material fields, and the forms of the respective kinetic and potential terms.To this end, the shear volume of dark energy models that have been explored is extensive: a by no means exhaustive list includes canonical quintessence, k-essence, α-attractors, f(R) gravity, Horndeski scalar-tensor gravity, DHOST theories, Einstein-Aether theories, bi-metric theories, and many more [8,9,[26][27][28][29][30][31][32][33][34].Furthermore, many of these theories have a large number of specific realizations.
It is hard to overstate the value of gaining such access to the fundamental microphysics of dark energy.For example, if dark energy is caused by some exotic field, this could allow us to situate dark energy within the standard model of particle physics.Such information might give us crucial information regarding persistent conundrums such the Hierarchy problem or what further symmetry principles are at play in particle physics [35,36].If modifications to gravity are at play, this could give us information regarding gravity's renormalizability [29].Cosmologically speaking, it could indicate the ultimate fate of the universe itself: will the universe continue in an accelerating expansion forever (as in the cosmological constant scenario or particular "freezing" realizations of dark energy) or will the universe stop accelerating (as in "thawing" realizations of dark energy) [14]; or could it even potentially begin contracting [37,38]?
In this article, we will attempt to shed some further light concerning the degree to which constraining and determining the values of the parameters w 0 and w a will inform us about the fundamental microphysics driving dark energy.Our verdict is largely pessimistic: it is unlikely that constraining the (w 0 , w a ) plane will ever allow us to single out a specific theory of dark energy, even if we were to detect clear evidence of dynamical behavior (w 0 ̸ = −1 and w a ̸ = 0).
We support our argument in two primary ways: (i) we demonstrate that there is significant underdetermination in the (w 0 , w a ) plane in the sense that there is a complete degeneracy between multiple realizations of dark energy and their mapping over vast swaths of this parameter space; (ii) we highlight some shortcomings of using (w 0 , w a ) by demonstrating that this parameterization is remarkably sensitive to the properties of the data used to constrain it, introducing additional confounding factors.
On (i), it has been previously suggested in many places in the literature that typical realizations of "freezing" and "thawing" models of dark energy occupy small, welldefined regions in the (w 0 , w a ) plane; and that measurements of these parameters would clearly indicate whether or not dark energy was described by a simple realization (e.g.single field, minimally coupled, canonical, etc.) of one of these classes [14,[39][40][41][42]. Furthermore, the typical understanding holds that finding values for w 0 and w a outside these narrow regions would indicate more complicated dynamics, exotic physics, or highly unnatural fine-tuning.
We challenge this orthodoxy specifically regarding thawing dark energy.We do so by exploring arguably the simplest model of dark energy, a minimally coupled scalar field with a quadratic potential.We analytically demonstrate, using exact solutions in a matterdominated background, that this model can arbitrarily sweep across huge sections of the (w 0 , w a ) plane depending on simple choices of model parameters.Furthermore, this model is of particular interest because it can be understood as an effective field theory approximation of a significant number of other distinct models [43]; this means that it provides a map between any space on the (w 0 , w a ) plane that it sweeps and many distinct models.We then show that these conclusions hold more generally where we numerically integrate the equations for a mixed dark matter/dark energy universe thought to describe the actual universe we live in.Indeed, the difficulty of constructing a unique potential from observational data has been noted in other places (see e.g.[22]).Here we analytically and numerically illustrate that determining a unique dark energy is almost certainly impossible by using this simple model (and the many models that it approximates) to cover significant portions of the remaining viable parameter space.
Regarding (ii), we demonstrate that the w 0 , w a parameterization is very sensitive to the range of redshifts one fits over.That is, we show that this parameterization captures some w(a) evolutions for the quadratic model reasonably well, but fairs far less successfully with others; making the mapping between theoretical dark energy models and the (w 0 , w a ) phase space sensitively dependent on arbitrary choices regarding the range of redshifts one fits over.We illustrate this effect with a few different choices of survey parameters.
The paper proceeds as follows.Sect.(II) provides an overview of quintessence: we introduce two flavours of the quadratic potential ("slow-roll" and "hilltop" variations), derive their analytic solutions, and discuss how this model can effectively approximate a huge variety of models with distinct potentials.Sect.(III) derives analytic expressions relating the quadratic potential to the (w 0 , w a ) parameter space and demonstrates that this model can sweep huge portions of the (w 0 , w a ) plane.We then numerically solve the equations of motion and show that this conclusion holds in a realistic universe with both dark matter and dark energy.Sect.(IV) discusses these results in light of current surveys aimed at constraining w 0 and w a .Sect. (V) concludes.

II. QUINTESSENCE
Quintessence is perhaps the simplest theoretical proposal for dark energy (see e.g.[44,45] for early papers on the subject), and is characterized by a dynamical scalar field φ with a canonical kinetic term and a potential V (φ).The theory is given by an action of the following form: where M pl is the reduced Planck mass, g is the determinant of the metric g µν , R is the Ricci curvature scalar, and S m is the action for matter.The scalar field minimally couples to gravity through the metric determinant and is assumed to have no direct coupling to matter fields.
FIG. 1.Current constraints on the (w0, wa) parameter space.Broadly speaking, this parameter space can be divided up into a few regions."Freezing" quintessence corresponds to the region of the parameter space where wa > 0 and refers to models where the dark energy equation of state is evolving asymptotically towards wDE ≃ −1, while "thawing" quintessence corresponds to the region where wa < 0 and refers to models where the dark energy equation of state is evolving away from wDE ≃ −1 to less negative values.The region given by w < −1 is known as the phantom region.It requires more exotic physics to describe and will not be a focus of this study.As we can see, most of the viable parameter space is locating in the thawing region.
The equation of state for dark energy w DE driven by such a field is given by: where P φ is the pressure and ρ φ is the energy density of the scalar field φ.Here, we can see that w DE will dynamically evolve in time with the evolution of the scalar field.
The evolution of the scalar field is given by the scalar field equation of motion: where V ′ (φ) = dV /dφ and H gives the expansion rate of the universe through the first Friedmann equation: where a is the scale factor of the universe and ρ is the density of matter and radiation.Together, these equations completely determine the dynamics of quintessence.Solving for the dynamics of φ allows one to then determine the evolution of w DE through Eq. ( 4).
The most extensively studied quintessence models fall into two broad categories: freezing quintessence or thawing quintessence.As the names suggest, freezing quintessence describes dark energy evolution w DE which was different from w DE ≃ −1 in the past, but is now evolving asymptotically towards this value as the universe expands (i.e."freezing in") and broadly speaking falls in the w a > 0 part of the parameter space; thawing quintessence describes dark energy evolution where w DE has been close to w DE ≃ −1 in the past, but is now beginning to evolve towards larger values as the universe expands (i.e."thawing out") and broadly speaking falls in the w a < 0 region of the parameter space.Here, we will be particularly concerned with thawing models as this is the less constrained part of the parameter space.
So far everything has been completely general.One must then specify a form of the potential function V (φ), which will enable one to solve the equations of motion, as well as give a particular realization of the relevant microphysics driving dark energy and allow for a direct investigation of any observables associated with the specific model(s).Mirroring a similar situation with inflation [46,47], the number of distinct models that are distinguished by their potential function is large.The possible forms of the potential include all imaginable varieties and combinations of power laws, inverse power laws, exponentials, axions, trigonometric functions, and many more.See e.g.[8,27,36] for reviews on quintessence which mention many of these potentials and [46] for an exhaustive review on inflationary potentials, many of which can be similarly adapted to quintessence.Furthermore, the exact form of the potential holds important information regarding the microphysics responsible for dark energy, the overarching structure and integration of such a field into the standard model of particle physics, and the ultimate fate of the universe's evolutionary trajectory.
However, in this article, we will only consider arguably the simplest model of quintessence: dark energy driven by a quadratic m 2 φ 2 potential.The reason for this is that this model can be understood from an effective field theory perspective as the leading order expansion of many other distinct models.In other words, an arbitrary (analytic) single-field model can be represented by an expansion: ) However, this can often be cast in a form that resembles the quadratic potential.For example, consider another well-studied model described by an exponential potential [45,[48][49][50].Expanding, we have that, However, a simple field re-definition φ → φ − φ 0 where φ 0 is constant allows us to cast this in a quadratic form resembling the m 2 potential.That is, where φ 0 can be chosen such that φ 0 = 1/λ and the linear term vanishes, and any constant terms can be absorbed into the definition of V 0 .Thus, we see that an exponential potential can, under certain circumstances, be well-approximated as a quadratic potential of the form , where in this case the "effective mass" is given by λ 2 V 0 = V ′′ .Similarly, one could consider the potential of pseudo Nambu-Goldstone boson given by [54], and find that under this expansion these models are also approximated by a quadratic potential with a negative sign in front of the mass term Given that this is an effective parameterization in a local region of the potential, we do not need to concern ourselves with the fact that V (φ) may be unbounded from below.
If the m 2 φ 2 model can cover large swaths of the parameter space in dark energy observables, this also implies a significant underdetermination in the form of the potential and the underlying microphysics driving dark energy.That is, any dark energy trajectories that can be obtained by an m 2 φ 2 model can also trivially be recast as any number of other dark energy models (admitting such an expansion of course) with distinct functional forms of their potentials and different fundamental microphysics.As indicated above, in the following we will consider two variations of the quadratic model: (i) the positive m 2 φ 2 model, which is representative of the standard slow-roll thawing quintessence models and (ii) the negative m 2 φ 2 model, which is representative of so-called "hilltop" quintessence models.Let us now flesh out our understanding of these two branches of the quadratic potential.
We begin with the positive quadratic model as this is representative of the most standard kind of thawing quintessence.Standard thawing models are characterized by a set of "slow roll" conditions (not to be confused with the inflationary slow roll conditions) [55].They are: These conditions are important because they ensure that the potential is flat enough to yield w DE ≃ −1 at early times, before the field starts slowly rolling, driving w DE to larger values at late time.
As mentioned earlier, our representative of standard, slowly rolling thawing models will be given by the quadratic potential: Another benefit of this particular model is that, depending on the dominant background energy component, Eqs.
(5-6) can be solved exactly (see e.g.[56][57][58]).Here and throughout this article, we work in dimensionless where H 0 is the Hubble constant today.We also take the ansatz a(t) ∝ t p which is true in both radiation and matter domination, leading to exact analytic expressions for φ: where J n and Y n are Bessel functions of the first and second kinds respectively and n = (1/2) 9p 2 − 6p + 1.
To simplify, we work in a matter dominated background (p = 2/3, n = 1/2).As we will show later though, our results are generalisable to the accelerated expansion era.The general solution for φ(t) is of the form: The particular solution φ(t) for the initial value problem with an initial value φ i , an initial velocity φi , at some early initial time t i (where we have written this using the small parameter ξ = mt i ≪ 1) is: Expanding in ξ, to leading order we find φ(t): where in the last line we have expanded again (this time in mt) as we have used that mt < 1 so that the scalar field has not entered the oscillatory regime (i.e. this region of the potential will not produce w DE ∼ −1).Eq. ( 16) will allow us to write down an analytic expression for Eq. ( 4) for this quintessence model.
While the slow roll conditions given in Eqs.(10)(11)) are sufficient to generate a typical thawing quintessence model, they are not both strictly necessary.As described first in [51], one can relax the condition given in Eq. (11) and still maintain the qualitative features of thawing quintessence, but with some notable differences.In this scenario, the scalar field rolls down a local maximum of the potential.Here, the model is such that the field remains close enough to the local maximum that Eq. ( 10) is still valid, but not necessarily Eq. ( 11).This has been dubbed hilltop quintessence [51,59] in analogy with similar hilltop models of inflation (see e.g.[60,61]).
We will examine a hilltop quintessence model given by a quadratic potential both because (as before) there are simple analytic expressions for the scalar field dynamics and this model approximates a tremendous variety of distinct dark energy models as their leading order expansion.The potential is given by: Eqs.
(5-6) can again be solved exactly for matter and radiation dominated backgrounds [62].The only difference is the presence of the minus sign in front of m 2 ; resulting in similar analytic solutions this time given by where I n and K n are modified Bessel functions of the first and second kinds respectively.In the matter dominated era, the general solution for φ is: The particular solution φ(t) for the initial value problem with an initial value φ i , an initial velocity φi , at some early initial time t i (where we have written this using the small parameter ξ = mt i ≪ 1) is: Expanding in ξ, to leading order we find: Notice that these are hyperbolic functions so we do not need to demand that mt < 1 to avoid the oscillatory regime.Similarly, Eqs. ( 21) will allow us to write down an analytic expression for Eq. ( 4) for this hilltop model.Having obtained solutions for the scalar field, we will now proceed to analyzing the evolution of the dark energy equation of state in these respective models.

III. DARK ENERGY IN THE (w0, wa) PLANE
A. The w0-wa parameterization revisited As mentioned earlier, the favoured parameterization for dark energy is given by: Before we proceed, we need to be clear about how w 0 and w a are actually determined.In practice, a given choice of w 0 and w a will determine the time evolution of the dark energy density, ρ DE or, in the case of quintessence, ρ φ .This density, via the Friedmann equations, will determine the expansion rate of the Universe as a function of time.Given a set of observations -typically measurements of standards candles over a redshift range -one can pin down values of the expansion rate at different times, or redshift.One then finds the w 0 and w a (and their associated uncertainties) that best fit the observations.In other words, in practice, w 0 and w a arise from fitting the data over a range of redshifts.
We can consider another parametrization of the equations of state, One can then calculate the values of ( w0 , wa ) for different choices of potentials, background expansions and initial conditions.This can give us a useful, often analytic, understanding of their features and how they relate to the underlying theory.But it should be clear, from the outset that the values obtained in this way, while indicative, will not be the values one obtains through the fitting procedure described above.We will bear this in mind as we proceed in what follows and we will be particularly careful, by using this notation, to distinguish between the two different "types" of (w 0 , w a ).
Let us now explore how the dynamics for different dark energy models (as well as their mapping into the ( w0 , wa ) parameters) can be straightforwardly understood from the scalar field equation of motion, φ + 3H φ + V ,φ = 0, where we have an "acceleration" term, a "friction" term, and a "potential" term.These terms will be important at different epochs depending on the classification of dark energy model.For example, in freezing models at early times the evolution of the potential term is significant as w DE ̸ = −1.Yet, as a freezing model approaches w DE ≃ −1, the potential becomes very flat and the dynamics are dominated by the friction term.The situation is different with thawing models as at early times the equation of state is locked in at w DE ≃ −1 and the friction term is dominant.As a thawing model evolves away from w DE ≃ −1, the potential term becomes more important as it is the field's evolution through its potential that causes dark energy to thaw.This leads to well-defined bounds that typical (i.e.slow-roll) freezing and thawing models respectively are thought to live in throughout the course of their evolution, where these bounds only span a tiny subset of the broader freezing (w a > 0) or thawing (w a < 0) regions.In particular, the thawing bounds are given by −1 ≲ wa /(1 + w0 ) ≲ −3 and it is commonly accepted in the literature that thawing quintessence models are constrained to live within this range (see e.g.[14,[39][40][41][42]); which clearly represents only a tiny fraction of the broader thawing region depicted in Fig. (1).These bounds are largely determined by the background in which the scalar field evolves in.Essentially, one can examine the ratios between the friction and acceleration terms and the potential and acceleration terms, and conclude that a slow-roll thawing model will be one for which this lower bound obtains when in a radiation dominated background and this upper bound obtains when in a matter dominated background (see e.g.[14,39,42] for further discussion).Thus, a thawing quintessence model in a matter-dominated universe will consequently evolve along a narrow line given by wa /(1 + w0 ) ∼ −3.
However, most of the available parameter space in the ( w0 , wa ) plane clearly lies outside of this given range.The question is, can we find a way to use the V 0 ± m 2 φ 2 models to fill in any of the remaining parameter space outside of the typically quoted thawing bounds?The answer, as we shall see below, is yes.This is because hilltop models of the type considered above have notably different evolutionary trajectories than the more standard slow-roll thawing models [51].
B. Analytic expressions for w(a) In order to determine the evolution of w(a) for the positive quadratic model, we utilize Eqs.(12,16) and replace them in Eq. ( 4) to get, In order to determine the expression as a function of the scale factor, we recall that during matter domination t = t 0 a 3/2 with respect to some reference time t 0 .
As expected, the slow roll m 2 φ 2 model evolves along this previously quoted line in the ( w0 , wa ) plane with a slope of −3.
Also utilizing t = t 0 a 3/2 and defining ϵ = (mφi) 2 2V0 , we now have that 1 , w(a) = −1 + 2ϵ a 6 (mt 0 ) 4 sinh mt 0 a 3/2 −mt 0 a 3/2 cosh mt 0 a 3/2 2 (27) and that, At a = 1 and using Eq. ( 2), the expression simplifies considerably, giving: The behavior of the hilltop model will be notably different than that of the slow roll model.For example, the behavior for small m (or equivalently V ′′ ≪ 1) approaches wa ≈ −3(1 + w0 ).This is not surprising as this is the regime in which the hilltop model approximates the slow roll conditions in the previous model.However, as one violates these conditions with larger m, the trajectory of w(a) can change substantially.A consequence of this is that these models can now look significantly different in the ( w0 , wa ) plane and one can quite easily tune the slope to be steeper depending on the choice of m.For example, Eq. ( 29) indicates that the slope would span between −3 ≲ wa /( w0 + 1) ≲ −15.5 for .01≤ mt 0 ≤ 6.This (i) indicates that we can use the quadratic hilltop model to essentially find any value we want in the ( w0 , wa ) plane simply by selecting appropriate mass/V ′′ and V 0 parameters.Furthermore, 1 See [51,59] for derivations of related, but more complicated expressions for w(a) in hilltop models.The expressions derived there represent approximate analytic solutions for w(a) assuming that the universe evolves in a dark energy dominated background.In this paper, we will focus on the exact analytic solution in a matter dominated background and the full numerical solution for a mixed dark matter/dark energy universe.
due to this model's ability to approximate any other model that admits a Taylor expansion of the type considered in Eq. ( 7), this (ii) similarly indicates that we can map any distinct dark energy model within this family (quadratic, exponentials, axions, etc) to any point in the ( w0 , wa ) plane that the quadratic model can reach through any one of these models' "effective" mass and V 0 terms.Taken together, these point to a potentially serious underdetermination in the microphysics of dark energy with respect to our observable parameterization of w(a).It seems like no matter where we end up in the ( w0 , wa ) plane within the broader thawing region, there are always a multitude of models that can be easily mapped to any point with a simple choice of parameters.

C. Covering the ( w0, wa) plane
Let us now see more concretely how these models map into the ( w0 , wa ) plane.To begin with, it is instructive to plot w(a) for different choices of parameters.If we first consider the models with a positive quadratic term (depicted in Fig. (2a)), we see that, for a fixed w0 , the evolution of the different solutions all converge to a single trajectory and seem to have the same slope as a decreases.This is qualitatively different for what happens for models with a negative quadratic term; in Fig. (2b), we plot a few models with the same w0 but with a range of different slopes, as a decreases (see also [51] for further excellent discussion on the w(a) evolutionary trajectories of hilltop models).This is a strong indication of what one should expect when mapping each of these theories onto the ( w0 , wa ) plane.The positive quadratic model will be locked in a narrow line along this plane as its parameters are varied, while the negative quadratic model will sweep out across the plane as its parameters are varied because these variations will create substantially different w(a) trajectories.
In order explore the ( w0 , wa ) plane, we first use the analytic expressions we have derived, and sample over a range different parameter values to fill in the ( w0 , wa ) plane.As expected, the slow-roll model lies on the line with a slope −3, depicted in Fig. (3).This model (and any others that it approximates) are locked in this narrow region of the w0 , wa parameter space.
The hilltop version of the quadratic model, on the other hand, picks up where the slow-roll model leaves off and sweeps over the steeper ranges of the ( w0 , wa ) plane as Fig. (3) shows.This model (and similarly any of the other many models that it approximates) can occupy large swathes of parameter space and in principle can sweep up to the "phantom" line ( w0 < −1), at which point more exotic physics would clearly be required.
Our analytic results for a matter dominated universe are already a clear indication that we can obtain almost any value of ( w0 , wa ) with this simple model.Given that, under general conditions for many models of thawing quintessence, one can express the Taylor expansions of  their potentials as the potential described by this simple model, this is confirmation that these models are underdetermined.In other words, many models will lead to the same observational results.
Until now we have based our reasoning on analytic expressions for the Taylor expansion of w(a) in the matter era; they tell us exactly what drives the changes in slope in terms of the parameters of the microphysics theory.But, of course, the universe is no longer matter dominated so it will be important to make sure that these conclusions hold up under a more realistic scenario.To check that is the case, we numerically integrate Eq. ( 5) for a mixed dark matter/dark energy universe.Instead of assuming a pre-determined evolution for a, we solve the Friedman equation in the presence of the scalar field so that We choose initial conditions (at some suitably early time) such that φ ≃ 0, and a value of φ i that leads to Ω DE ≃ .7 and H = 1 at the end of the integration (at a = 1).
We then generate w(a) numerically for a wide variety of parameters and find ( w0 , wa ) for both models.The numerical integration for the positive quadratic, slow-roll model yields a similar result: this model lives on a very thin strip of parameter space.However, the line it traces out is wa /( w0 + 1) ∼ −1.5 (again at a = 1), rather than wa /( w0 + 1) ∼ −3 as in the matter dominated case.This is because in a universe with a substantial dark energy component, the Hubble friction term will be enhanced; meaning that the evolution of w(a) will be suppressed when compared to the matter dominated case for the same choice of parameters (see Fig. (4a)).Notice also, that this numerical integration shows that for a mixed dark matter/dark energy universe, the evolution of w(a) for the positive quadratic model is far more linear than the matter dominated universe.Thus, the resulting line behavior in the ( w0 , wa ) plane will be less than the matter dominated value.This wa /( w0 + 1) ∼ −1.5 result for thawing quintessence in a mixed dark matter/dark energy universe is consistent with what was found in [16,25,55,63] for similar kinds of models that can be said to represent more typical realizations of thawing quintessence.
The numerical integration for the negative quadratic, hilltop model in a mixed dark matter/dark energy universe also yields a similar result when compared with the analytic solution for a matter dominated universe: this model sweeps out a wide swath of the ( w0 , wa ) plane.As with the positive quadratic model, the friction term is enhanced meaning that the evolution of w(a) is suppressed FIG. 4. w(a) numerical integration for both the slowroll/positive quadratic model and the hilltop/negative quadratic model in a mixed dark matter/dark energy universe compared to their respective matter dominated analytic solution for the same choice of parameters.In both cases, we see that the numerical solution does not evolve quite as much.This is because the Hubble friction term is enhanced when dark energy becomes a significant part of the scalar field equation of motion when compared with the matter dominated solution.
when compared with the matter dominated case for the same choice of parameters (see Fig. (4b)).Here, the evolution of w(a) for the negative quadratic model is still non-linear, but the evolution is not quite as steep.Thus, we expect the resulting area that the model sweeps in the ( w0 , wa ) to be somewhat smaller to the matter dominated case.Here, we find, for this numerical evolution and mt 0 ∈ [.01, 6], −1.5 ≲ wa /( w0 + 1) ≲ −10 (Fig. (5)).
Thus we have shown, both analytically and with numerical solutions, that we can generate an incredibly broad family of possible behaviours for the equation of state, w(a).Specifically, we have shown that we can get nearly any arbitrary value of w0 and wa with a quadratic model for the dark energy potential, both analytically in a matter dominated universe and numerically in a mixed dark matter/dark energy universe thought to be a good description of our current universe.This means that our w0 wa FIG. 5.This is the result for w0 ≡ w(a = 1) and wa ≡ −dw/da(a = 1) determined after numerically integrating Eq. ( 30) for a universe with ΩDE ≃ .7 for both the slow-roll/positive quadratic model and the hilltop/negative quadratic model.Both models are qualitatively similar in their behavior, with the slow-roll model now living on wa ∼ −1.5( w0 + 1) while the hilltop model takes values roughly between −1.5( w0 + 1) ≲ wa ≲ −10( w0 + 1) and sweeps across the parameter space.
conclusions implying a significant underdetermination of the microphysics underlying dark energy with respect to the observables w0 and wa also hold in a realistic description of the universe, as these numerical solutions for the quadratic model will similarly map many distinct dark energy models onto the exact same regions of the ( w0 , wa ) plane just as well as the analytic solutions in a matter dominated universe.Furthermore, it is also worth emphasizing that these results can be understood as being broadly consistent with some other studies in the literature that have considered the ( w0 , wa ) parameterization from different perspectives than the one we have adopted here.For example, [25] takes the CPL parameterization as a starting point in order to determine which scalar field models can reasonably correspond to it and finds that a number of potentials V (φ) are consistent with wa ∼ −1.5( w0 + 1) because many potentials will exhibit sufficiently linear behavior when φ only rolls over a small enough region of the potential.Another interesting study [24] adapts the flow equation formalism from inflation to a very generic description of quintessence and also concludes that quintessence models can be found in many regions in the ( w0 , wa ) plane outside the typical freezing and thawing bounds; where they calculate ( w0 , wa ) from the best-measured principal components (PCs) (see [64]) of the w(a) trajectory.

IV. FITTING w0 AND wa
We do not directly measure w0 and wa .In practice, as discussed above, we fit w 0 and w a over a range of redshifts.The parameter space depicted in Fig. ( 1) is determined in this way.Typically one has a range of distance measurements -such as the angular diameter distance or luminosity distance of structures or objects -which are integrals of the expansion rate and which, in turn, are a function of the energy density of dark energy.One then infers the properties of the dark energy component from how well a given model fits the distance measurements.
There are a number of Stage IV surveys that have been proposed to further our understanding of dark energy [65][66][67].To simplify, we emulate a typical stage IV survey and assume that the results are measurements of the Hubble parameter as a function of redshift.Our conclusions would be no-different if we had used distance measurements themselves.As an example, we take the redshift range and forecast uncertainties in [68].
The model we need to fit is the the Friedman equation in the form 3(1+w0+wa) .
(31) And we approximate the negative logarithm of the likelihood as where H obs (z) is the observed H, H(z) is the computed H from Eq. ( 31) as a function of the w 0 and w a parameters, and σ i refers to the uncertainty at the redshift bin i. Assume now that the observed H(a) corresponds to a particular model of thawing quintessence for which we have full numerical solutions.Let us then generate numerical evolutions for w(a) and determine the best w 0 , w a fit using Eqs.(2,31,32) for redshifts z ∈ [.15, 1.85] (corresponding to the redshift bins from [68]).We begin with the positive quadratic model first and then proceed to the negative quadratic model.
For the positive quadratic model, this fitting procedure produces results that are very similar to what we found earlier: the w 0 , w a values for the quadratic model lie on a very narrow strip of the (w 0 , w a ) plane right around w a /(w 0 + 1) ∼ −1.5 (again consistent with other studies such as [16,25,55,63]).Upon reflection, this should not be terribly surprising.After all, as we have seen in Fig. (4a), the numerical solution w(a) for the positive quadratic model in a mixed dark matter/dark energy dominated universe like our own is highly linear.Thus, its mapping into the (w 0 , w a ) plane will not change substantially between only taking these values at a = 1 or fitting them over a quite significant range of redshifts.
The negative quadratic model, on the other hand, maps quite differently into the (w 0 , w a ) plane when this fitting procedure is employed.Qualitatively, we still see that various parameter choices for this model will produce a 'cloud' in the parameter space (rather than a narrow line as with the positive quadratic model).However, this cloud is significantly more narrow than we saw in the a = 1 case as it very roughly lies between −1.5 ≲ w a /(w 0 +1) ≲ −2.5.5).This indicates that there is an ambiguity in how these hilltop models are represented in the (w0, wa) plane as the exact size of the swept region (for the same choice of parameters) will sensitively depend on the range of redshifts that one fits over due to the highly nonlinear evolution of w(a) in these hilltop models.Fitting over a more restricted range of more recent redshifts would cause the fitted results to more closely resemble the a = 1 results.By contrast, the slow-roll models (dark blue line) still lie in the same narrow strip as their location is insensitive to the choice of fitting procedure as their evolution is highly linear.plane for the same numerical evolution as Fig. (5), except that w 0 and w a have been determined by fitting over z ∈ [.15, 1.85] rather than by determining the values w0 and wa at a = 1; where as we saw this choice of parameters previously swept across −1.5 ≲ wa /( w0 + 1) ≲ −10.Fig. (7) depicts how the w0 , wa values map onto the fitted parameters w 0 and w a ; and there we can clearly see how this fitting procedure 'squeezes' the (w 0 , w a ) phase space that these hilltop models live in.
While this does not change the broader conclusionthat there is a significant amount of parameter space for which the microphysics of dark energy is severely underdetermined because many distinct microphysical models can sweep significant parts of this parameter spacethis does reveal that exactly how a model of dark energy maps into the (w 0 , w a ) plane can potentially be sensitively dependent on somewhat arbitrary choices for fitting procedure.This surprising sensitivity can be traced to the fact that the w(a) evolution in the hilltop model has significant non-linearities.Consequently, as a linear parameterization, w 0 and w a may not necessarily provide a good description of a dark energy model that does not feature a fairly linear evolutionary trajectory w(a) (see also [25,51]).This immediately allows us to understand that the w 0 , w a parameterization will not provide a good description over that range of redshifts for the hilltop model, because the vast majority of w(a)'s evolution is Clearly the redshift range of the data plays a significant role in the range of w 0 , w a one obtains from thawing quintessence.Paradoxically, but not surprisingly, Fig. (8) shows that, if one narrows the range of redshifts (in this case from z ∈ [.15, 1.85] to z ∈ [0.15, 0.25]) we can see that we get a closer approximation with the w0 , wa parameterization to the full w(a) evolution.Consequently, one can repeat the procedure we have done in this section, but fit over a smaller range of redshifts.One then finds that doing so will result in a wider sweeping of the (w 0 , w a ) plane that much more closely resembles ( w0 , wa ) the more one restricts the range of redshifts to those where the most significant w(a) evolution occurs as the fit now "sees" the steeper part of the trajectory.The reason one might want to do this is that, as Fig. (8) clearly shows, this is a far better parameterization of the actual w(a) trajectory.However, in practice, the attendant costs of doing so would require utilizing less survey data and significantly increasing the uncertainties of the observations.With all these considerations though, it is clear that there is some ambiguity in how hilltop models of the  6) is squeezed when fit over a large range of redshifts: the fit for w0 and wa in that case does not fully capture the steeper part of the w(a) trajectory that occurs at the more recent redshifts and consequently maps these model into a less steep part of the (w0, wa) plane.The smaller and more recent the range of redshifts included in the fit, the closer that the fitted (w0, wa) values will be to the ( w0, wa) values determined at a = 1.
type considered here will map into the (w 0 , w a ) plane as different fitting procedures will result in the models sweeping different areas of the plane.Consequently, we should seek alternative parameterizations for the microphysics of dark energy that are not so dependent on such a choice of fitting procedure.

V. CONCLUSION
One of the goals of modern cosmology is to determine the microphysical model that underpins the accelerated expansion of the Universe.The most popular proposal is that it is driven by some form of dark energy which can be characterized by an equation of state.The guiding principle has been that, with current and future cosmological observations, we will be able make accurate measurements of the equation of state and, as a result, pin down the right microphysical model of dark energy.In this paper, we have focussed on a very broad class of models of dark energy -thawing quintessence -and showed that future observations will inevitably underdetermine the microphsyics of dark energy.
By focusing on a widely used parametrization of the cosmological dark energy -w 0 and w a -we have shown that we can get almost any value of these parameters with a simple quadratic model for the potential.Thus, contrary to what has often been claimed in the literature, simple realizations of thawing quintessence are not con-fined to a small region of the (w 0 , w a ) phase space.Furthermore, this does not require any highly exotic physics, unusual fine-tuning, or inordinately complicated dynamics.Essentially, one just needs a simple single field model of canonical quintessence with a quadratic potential, and in particular the hilltop flavour of this model which was first investigated in detail by [51].As we have seen, this model can find nearly any spot in the thawing region of the (w 0 , w a ) plane with a judicious choice of basic model parameters.We showed this using exact solutions for φ in the matter dominated universe as well as numerically integrating the scalar field equations of motion for a mixed dark matter/dark energy universe.
This indicates that there is a significant underdetermination with respect to this model and any other model that can be placed anywhere within the (w 0 , w a ) plane that this quadratic model can sweep.As we demonstrated earlier, we already know of several distinct microphysical models of dark energy that can be mapped into this exact same region of the parameter space because the quadratic model approximates all of the many models that admit of a Taylor expansion of the form Eq. ( 7) where the potential is expanded to quadratic order.Thus, these models will be indistinguishable from the quadratic model considered here from the point of view of observables in the (w 0 , w a ) plane.
Given that we can fill out the parameter space with the quadratic model and map between that model's predictions for the parameter space and many other models, this deflates some of the motivation for investigating models with different potentials.There could certainly be interesting or even compelling theoretical reasons or non-empirical motivations for pursuing specific models (coherence, explanatory power, aesthetics, problemsolving capabilities etc; see e.g.[69][70][71][72][73][74][75][76][77]).For example, there are specific potentials that are pursued for their ability to resolve outstanding fine-tuning problems [78] or that possess particularly attractive theoretical qualities such as having radiative stability due to their symmetries [53].However, it is unlikely that strictly empirical methods will single out a unique potential.
Our work does not rule out the possibility of structure formation [79][80][81][82][83], gravitational waves [84][85][86][87][88][89][90], fifth force tests [91][92][93] etc. pointing us towards more specific microphysical realizations of dark energy.For example, non-minimal couplings could conceivably show up in any of these different types of measurements and could be used to narrow down the possible microphysical models one might consider [18].Furthermore, if we are interested in learning something about the the microphysics of dark energy, we should seek alternative parameterizations of w(a).As we have seen with the hilltop models considered here, the w 0 , w a parameterization of dark energy is ambiguous in terms of how it represents certain classes of dark energy models in its parameter space due to these parameter's sensitive dependence on distance measurements.Parameterizing in terms of (V 0 , V ′′ ) could potentially be a more powerful approach and shed further light on the fundamental mechanism at play which is driving accelerated expansion [11,22].

FIG. 2 .
FIG.2.w(a) evolution for different choices of model parameters for the slow-roll/positive quadratic model and the hilltop/negative quadratic model.The parameter mt0 was varied and X was chosen so that w0 ≃ .90.For the slow-roll models, all converge on the same evolutionary track and will have the same slope in the ( w0, wa) plane.All the hilltop models have very different evolutions and the trajectories become steeper as mt0 is increased.Consequently, they will be represented differently in the ( w0, wa) plane.

FIG. 6 .
FIG. 6.The light blue shaded region depicts the result for w0 and wa determined by numerically integrating w(a) for the hilltop model in a universe with ΩDE ≃ .7 and finding the best fit for Eqs.(31-32) over z ∈ [.15, 1.85].This is overplotted against the w0 ≡ w(a = 1), wa ≡ −dw/da(a = 1) result (dark gray shaded region) for the same hilltop model and choice of parameters depicted in Fig. (5).This indicates that there is an ambiguity in how these hilltop models are represented in the (w0, wa) plane as the exact size of the swept region (for the same choice of parameters) will sensitively depend on the range of redshifts that one fits over due to the highly nonlinear evolution of w(a) in these hilltop models.Fitting over a more restricted range of more recent redshifts would cause the fitted results to more closely resemble the a = 1 results.By contrast, the slow-roll models (dark blue line) still lie in the same narrow strip as their location is insensitive to the choice of fitting procedure as their evolution is highly linear.

FIG. 7 .
FIG. 7. Comparison between the wa and w0 values depicted in the (w0, wa) plane of Fig. (6).This shows how the values determined by fitting H(z) through Eqs.(31-32) and those determined by w0 and wa map onto each other.In other words, this fitting procedure squeezes the representation of the hilltop model in the (w0, wa) plane.The exact degree to which the swept region is squeezed depends on the range of redshifts included in the fit.

FIG. 8 .
FIG. 8. Evolution of w(a) compared with two different fits for the w(a) = w0 + wa(1 − a) parameterization.The yellow line was determined by fitting H(z) through Eqs.(31-32) for z ∈ [.15, 1.85], while the green line was determined by fitting from z ∈ [0.15, 0.25].We can see that fitting over the more recent range of redshifts better captures the w(a) evolution for the steeper hilltop models considered here.Consequently, we can understand why the swept area of Fig. (6) is squeezed when fit over a large range of redshifts: the fit for w0 and wa in that case does not fully capture the steeper part of the w(a) trajectory that occurs at the more recent redshifts and consequently maps these model into a less steep part of the (w0, wa) plane.The smaller and more recent the range of redshifts included in the fit, the closer that the fitted (w0, wa) values will be to the ( w0, wa) values determined at a = 1.