New constraints on the expansion rate of the Universe seem to favor evolving dark energy in the form of thawing quintessence models i.e. models for which a canonical minimally coupled scalar field has at late times begun to evolve away from potential energy domination. We scrutinize the evidence for thawing quintessence by exploring what it predicts for the equation of state. We show that in terms of the usual ChevalierPolarskiLinder parameters inlineformulammlmath displayinlinemmlmrowmmlmsubmmlmrowmmlmiwmmlmimmlmrowmmlmrowmmlmn0mmlmnmmlmrowmmlmsubmmlmrowmmlmathinlineformula inlineformulammlmath displayinlinemmlmsubmmlmiwmmlmimmlmiammlmimmlmsubmmlmathinlineformula thawing quintessence is in fact only marginally consistent with a compilation of the current data. Despite this we embrace the possibility that thawing quintessence is dark energy and find constraints on the microphysics of this scenario. We do so in terms of the effective mass inlineformulammlmath displayinlinemmlmsupmmlmimmmlmimmlmn2mmlmnmmlmsupmmlmathinlineformula and energy scale inlineformulammlmath displayinlinemmlmsubmmlmiVmmlmimmlmn0mmlmnmmlmsubmmlmathinlineformula of the scalar field potential. We are particularly careful to enforce uninformative flat priors on these parameters so as to minimize their effect on the final posteriors. While the current data favors a large and negative value of inlineformulammlmath displayinlinemmlmsupmmlmimmmlmimmlmn2mmlmnmmlmsupmmlmathinlineformula when we compare these models to the standard inlineformulammlmath displayinlinemmlmi mathvariantnormalmmlmimmlmiCDMmmlmimmlmathinlineformula model we find that there is scant evidence for thawing quintessence.
['Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Scant evidence for thawing quintessence
{'Scant evidence for thawing quintessence': "William J. Wolf, 1, ∗ Carlos Garc'ıa-Garc'ıa, 1, † Deaglan J. Bartlett, 2, ‡ and Pedro G. Ferreira 1, §  \n1  \nAstrophysics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH, United Kingdom 2 CNRS & Sorbonne Universit'e, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France  \nNew constraints on the expansion rate of the Universe seem to favor evolving dark energy in the form of thawing quintessence models, i.e., models for which a canonical, minimally coupled scalar field has, at late times, begun to evolve away from potential energy domination. We scrutinize the evidence for thawing quintessence by exploring what it predicts for the equation of state. We show that, in terms of the usual Chevalier-Polarski-Linder parameters, ( w 0 , w a ), thawing quintessence is, in fact, only marginally consistent with a compilation of the current data. Despite this, we embrace the possibility that thawing quintessence is dark energy and find constraints on the microphysics of this scenario. We do so in terms of the effective mass m 2 and energy scale V 0 of the scalar field potential. We are particularly careful to enforce un-informative, flat priors on these parameters so as to minimize their effect on the final posteriors. While the current data favors a large and negative value of m 2 , when we compare these models to the standard ΛCDM model we find that there is scant evidence for thawing quintessence.", 'I. INTRODUCTION': "A combination of the latest measurements of baryon acoustic oscillations (BAO) [1], luminosity distances of distant supernovae (SNe) [2-4], and measurements of the Cosmic Microwave Background (CMB) [5-8] have provided hints of deviations from the Λ-Cold Dark Matter (ΛCDM) paradigm. One possible interpretation is that the accelerated expansion of the Universe at late times is due to some form of evolving dark energy, such as a quintessence scalar field, with a time evolving equation of state.  \nMuch of the discussion in the literature has focused on the Chevallier-Polarski-Linder (CPL) parameterization [9, 10] of the equation of state of dark energy, w ≡ P φ /ρ φ (where P φ and ρ φ are the pressure and energy density of the dark energy), in which  \nw ( a ) = w 0 + w a (1 -a ) , (1)  \nwhere a is the scale factor of the Universe. w 0 gives the value of w today and w a characterizes its temporal evolution. ΛCDM is given by w a = 0 and w 0 = -1, while a variety of dynamically driven dark energy possibilities occupy the rest of the parameter space defined by the ( w 0 , w a ) plane.  \n̸  \nThe CPL parametrization is somewhat of a blunt instrument, a form of data compression which is amenable to over-interpretation [11]. Nevertheless, the constraints on w 0 and w a seem to tell an intriguing story about the microphysics of dark energy as is apparent on close inspection of Fig. 1. For a start, they seem to show that w a = 0, which would represent evidence for a time evolving equation of state. Furthermore, the locale in the ( w 0 ,  \nFIG. 1. Constraints on CPL parameters (68% and 95% C.L.) using the combined DESI BAO data [1], Pantheon+ SNe data [3], and CMB (Planck and ACT lensing) data [5-8].  \n<!-- image -->  \nw  \nw a ) plane which is singled out by the current data seems to favour 'thawing' models of dark energy, i.e., models which will have started off with w ( a ) ≃ -1 before evolving to the larger values that they have now. And finally, the inclination of the degeneracy of the ( w 0 , w a ) constraints penetrates any area corresponding to 'phantom' dark energy, i.e., dark energy in which the equation of state is w < -1.  \nWhile the new constraints on ( w 0 , w a ) are not yet definitive, with particular choices of data sets being questioned in the analysis [12, 13], their qualitative trend has led to an avalanche of results concerning dark energy [1447]. We highlight a few important comments and results  \nthat have arisen. First of all, and following the point made in Wolf and Ferreira [11], a few authors argued that the ( w 0 , w a ) determined from a given survey are tied to the specific properties of that survey. In other words, while ( w 0 , w a ) are, naively, often interpreted as the first terms of a Taylor expansion of the equation of state, they are in fact fitting parameters which depend on the survey properties. Cortˆes and Liddle [19] made the further point that the CPL parameterization was in fact showing that w ≃ -1 in the redshift range where there is most data. Shlivko and Steinhardt [18] echoed the point made in Wolf and Ferreira [11], that one cannot infer the microphysics of the dark energy model from the inferred ( w 0 , w a ) parameters and, for example, whether it is phantom or not because this involves extrapolating the best fit w ( a ) far outside the redshifts at which the data are used to determine the parameters. On the other hand, more general, model agnostic studies such as Ye et al. [34], Lodha et al. [20], and Calderon et al. [45] point towards a (marginally) statistically significant detection of phantom-like behavior for z > 0 . 6, where z is the cosmological redshift.  \nIn [11], it was argued that a particularly simple quintessence model, with a scalar field, φ , and a potential  \nV ( φ ) = V 0 + 1 2 m 2 φ 2 , (2)  \ncan, in principle, cover an incredibly wide region of the 'thawing' quadrant of the ( w 0 , w a ) plane. In essence, this potential represents the leading order terms in an effective field theory expansion for a huge variety of scalar field models in the region of field space in which it can have an impact on the dynamics of the Universe. Thus, it was argued, the vast majority of dark energy models under consideration (i.e. thawing models) are undetermined. In other words, with the types of cosmological measurement available, it will be impossible to identify the specific microphysics of the dark energy model, beyond that given by Eq. (2).  \nIn this paper, we study what the new constraints on the expansion rate of the Universe tell us about quintessence. More specifically, we directly map the dark energy model given by Eq. (2) into the ( w 0 , w a ) plane considering the data provided by BAO, SNe, and CMB measurements. We find that, from the point of view of the CPL parameterization, and once all of the data is accounted for, thawing quintessence (as represented by Eq. (2)) is not favored by the BAO, SNe and CMB data. Thus, thawing quintessence more generally lies mostly outside the favoured regions once it is directly fit to the data. Nevertheless, we can still examine the implications for thawing quintessence, and to this end, we directly constrain the microphysics of thawing quintessence in terms of V 0 and m 2 .  \nThe paper proceeds as follows. In Section II we review dynamical dark energy and thawing quintessence, and describe the different regimes it has, depending on whether m 2 is negative or positive. This allows us to  \nunderstand how the potential from Eq. (2) is able to cover such a broad range of thawing behaviour. In Section III we describe, in detail, how thawing quintessence maps onto the ( w 0 , w a ) plane, and how this mapping is dependent on the survey characteristics one is using to measure the CPL parameters. We show how, in the case that we consider a compilation of all data sets (including, in particular, the CMB) thawing quintessence is only marginally consistent with the current constraints on ( w 0 , w a ). In Section IV, we find constraints on the microphysics of thawing quintessence, taking particular care to understand the priors on V 0 and m 2 , and how they affect the posteriors we obtain. An important, novel aspect, of the analysis is our use of an analytic ansatz to find the initial conditions for the scalar field evolution, using Symbolic Regression techniques described in Appendix B. In Section V we discuss the implications of our findings. In a series of Appendices we lay out some of the details and, more significantly, new techniques developed for this analysis.", 'II. THAWING DARK ENERGY': "In this paper we will focus on dark energy driven by a quintessence scalar field and described by the action:  \nS = ∫ d 4 x √ -g [ 1 2 M 2 Pl R -1 2 g µν ∂ µ φ∂ ν φ -V ( φ ) ] + S m . (3)  \nM Pl is the reduced Planck mass, g is the determinant of the metric g µν , R is the Ricci curvature scalar, V ( φ ) is the scalar field potential, and S m is the action for matter. The scalar field has a canonical kinetic term and is minimally coupled to gravity through the metric determinant.  \nSuch a theory leads to a dark energy equation of state w ( a ),  \nw ( a ) = ˙ φ 2 2 -V ( φ ) ˙ φ 2 2 + V ( φ ) , (4)  \nwhich evolves in time with the evolution of the scalar field, where overdots are derivatives with respect to cosmic time. The evolution of the scalar field is given by the scalar field equation of motion,  \n¨ φ +3 H ˙ φ + V ' ( φ ) = 0 , (5)  \nwhere V ' ( φ ) ≡ dV/dφ and H is the expansion rate of the Universe given by the first Friedmann equation,  \nH 2 ≡ ( ˙ a a ) 2 = 1 3 M 2 pl ( ρ + ρ φ ) , (6)  \nwhere ρ φ and ρ are, respectively, the scalar field and ordinary matter energy densities and the Hubble parameter today is H 0 = 100 h km s -1 Mpc -1 with h ≃ 0 . 67.  \nThe behavior of w ( a ) is the most straightforward way to characterize dark energy, and one can identify different broad classes of dark energy models, such as thawing (when dw/da > 0) or freezing (when dw/da < 0) models [48]. In the case of a quintessence scalar field with a canonical kinetic energy, as in Eq. (3), the dynamical behavior of w ( a ) is largely determined by the potential V ( φ ). While a great number of potentials have been considered in the literature, the physics of thawing quintessence can be effectively captured by Eq. (2) where, throughout this paper, m 2 will be in units of ( H 0 /h ) 2 and V 0 in units of M 2 Pl ( H 0 /h ) 2 . This has do with the fact that any arbitrary analytic potential can be represented by an effective Taylor expansion in φ , combined with a field redefinition to get rid of the linear term [11]. For example, exponential tracker potentials [49-51] are well-described at leading order by the standard ( m 2 > 0) quadratic version of Eq. (2), while pNGB and axion potentials [5254] are well-described by the hilltop ( m 2 < 0) version of Eq. (2) (see e.g. [11, 54-56] for more on the quadratic hilltop potential). Consequently, the potential given in Eq. (2) can be understood as the leading order description for many potentials of great theoretical interest, and characterizes all of these distinct theories in terms of an energy scale V 0 and an effective mass m 2 = d 2 V/dφ 2 in the leading order quadratic term.  \nFurthermore, the potential of Eq. (2) captures the entire phenomenology of thawing quintessence if one thinks of it, purely, in terms of w ( a ). The m 2 > 0 branch of the model phenomenologically captures standard 'slowroll' thawing quintessence, where all of these models converge to a highly linear, universal dynamical behavior in the dark energy equation of state w ( a ) given by dw/da | a =1 ≈ 1 . 5[1 + w ( a = 1)] (see e.g. [11, 57-61]). On the other hand, the m 2 < 0 'hilltop' branch behaves quite differently as there are an exceedingly wide range of dynamically possible behaviors for w ( a ). For example, at small ∣ ∣ m 2 ∣ ∣ , the hilltop models approach the linear evolution of the m 2 > 0 models, while at larger ∣ ∣ m 2 ∣ ∣ the hilltop models can evolve in a highly non-linear manner, where w ( a ) is flat across most of the Universe's history but shoots sharply upwards at more recent redshifts (see [11, 18, 54, 62] for further discussion). Consequently, this model can in principle be mapped across the ( w 0 , w a ) plane depending on the choices of mass, energy scale, and initial conditions [11, 18]. Thus, the potential of Eq. (2) is an ideal model to consider constraints on thawing quintessence microphysics in terms of an energy scale and effective mass.  \n̸  \nIt is useful to understand the dynamics of thawing quintessence driven by the potential in Eq. (2) in slightly more detail. Without loss of generality, we will assume that ˙ φ = 0 initially. If m 2 > 0, we have that φ will start off with φ ini = 0 and a potential energy given by V = V 0 + 1 2 m 2 φ 2 ini . Initially, the scalar field will remain frozen due to Hubble friction, slowly developing a kinetic energy over time. Once the scalar field has sped up enough, w ( a ) will start to deviate from -1, and, at  \nsome point, the scalar field will end up oscillating around the minimum of the potential. For time scales such that mt ≫ 1, the equation of state is that of dust, w ( a ) ≃ 0 [63], but for the choice of parameters we will consider here, that is well beyond the time scales we will be looking at. Note that, for m 2 > 0, we will still recover a positive energy density and w ( a ) = -1 even if we set V 0 = 0. With both V 0 and m 2 > 0, we expect to have a degeneracy at early times between the two parameters. Indeed, we can define the effective cosmological constant at early times to be Λ = ( V 0 + 1 2 m 2 φ ini ) / (3 M 2 Pl ), making manifest the degeneracy between the two parameters. This degeneracy will play an important role in the understanding the posterior constraints on thawing quintessence in Section IV.  \nIf m 2 < 0, the situation is markedly different [11, 54, 62]. The potential is now a 'hilltop' and the solutions are unstable, exponential in | m | t . The larger | m | , the faster the instability. Now, one can now tune the initial conditions, placing φ ini arbitrarily close to 0 leading to an extended period in which w ( a ) = -1. But once the scalar field has grown enough, the onset of kinetic energy domination can be rapid. This translates in an equation of state which remains at -1 over most of the Universe's expansion history and then rapidly transitions to w > -1. The higher the mass, the later, but also more rapid, the transition away from w ∼ -1.  \n̸  \nThe case of m 2 = 0 merits a brief note. We are considering the case where ˙ φ ini = 0. If that is the case, and m 2 = 0, we have that V = V 0 throughout its evolution and we recover the case of a cosmological constant. We can assume that ˙ φ ini = 0 but we then have ˙ φ ∝ a -3 and is thus very rapidly redshifting away. So we have that m 2 = 0 effectively corresponds to w ( a ) = -1 throughout.", 'III. THAWING QUINTESSENCE FROM THE POINT OF VIEW OF ( w 0 , w a ).': "In this section we explore how thawing quintessence maps on to the ( w 0 , w a ) plane of the CPL parametrization. To do so, we need to recap how ( w 0 , w a ) are actually inferred from cosmological data. In practice, we do not directly measure w ( a ), but infer ( w 0 , w a ) from data that probe H ( z ) at different redshifts. In the CPL parameterization, H is given by  \nH 2 ( a ) = H 2 0 [ Ω m a -3 +(1 -Ω m ) e 3 w a ( a -1) a -3(1+ w 0 + w a ) ] , (7)  \nwhere Ω m is the fractional energy density in nonrelativistic matter today and H 0 is the Hubble constant. A dark energy model's representation in terms of ( w 0 , w a ) depends on how Eq. (7) compares with actual measurements of cosmological quantities sensitive to H ( z ).  \nWe now map all of thawing quintessence as given by Eq. (2) into the ( w 0 , w a ) plane given by the CPL parameterization. We do so by mimicking what is done with  \nFIG. 2. Equation of state, w ( a ), for a scalar field model with m 2 < 0 (( V 0 , m 2 ) = (0 . 9725 , -40 . 0), black solid line) alongside the corresponding best fit CPL models when using the BAO, SNe, and CMB survey characteristics (as described in Section III) both individually and combined, showing that the best fit ( w 0 , w a ) parameters are clearly sensitive to the survey parameters.  \n<!-- image -->  \nreal data. For a given choice of ( V 0 , m 2 ) we generate the associated H ( z ) and then the equivalent observables which are measured when one uses BAO, SNe and CMB data. Given these observables we then find the best-fit CPL parameters. In this paper we will focus on the BAO data from DESI [1], the SNe data from Pantheon+ [3], and the CMB data from Planck and ACT [5-8].  \nTo be more specific we determine the corresponding CPL parameters, ( w 0 , w a ), given the characteristics of a particular survey, 'S'. We characterize a given survey in terms of a set of observables O i at different redshifts, z i , and their associated covariance, Cov. We construct the distance between two models in terms of  \nχ 2 S = ( O CPL -O φ ) T ( O data O φ ) T Cov -1 O data O φ ( O CPL -O φ ) (8)  \nwhere O CPL and O φ are the observables computed with the CPL parameterization ( w 0 , w a ), or the scalar field ( V 0 , m 2 ), respectively, and the factors O data / O φ ensure that we fit the observables using the surveys measurements' relative errors, and not the absolute ones. Thus, for a set of thawing quintessence parameters, ( V 0 , m 2 ), we minimize Eq. (8) to find the corresponding CPL parameters, ( w 0 , w a ). So, for example, mapping ( V 0 , m 2 ) into the CPL parameters ( w 0 , w a ) for the DESI survey would involve minimizing χ 2 DESI from Eq. (8), while doing so for a combination of data sets can be done by minimizing χ 2 tot = χ 2 DESI + χ 2 SNe + χ 2 CMB .  \nFor this process to be at all feasible, we consider a reduced or compressed version of the data. In Appendix A, we explain how the data can be compressed into a set of 'measured' physical parameters and their covariance. In  \n,  \nbrief, the CMB data can be compressed into an acoustic scale ℓ A and shift parameter R ( z ∗ ) at the photon decoupling redshift z ∗ [64-66], the SNe data are sensitive to the luminosity distance and can be compressed into measurements of E ( z ) ≡ H ( z ) /H 0 at different redshifts [67], and the BAO data are directly sensitive to the angular diameter distance or Hubble parameter, which can be directly used in our fit (see Appendix A for details).  \nIn [11] we showed that, due to wide range of dynamical behavior in the evolution of w ( a ) in some of these models, the exact representation of these models in terms of ( w 0 , w a ) can potentially depend sensitively on which redshifts are probed in distance measurements sensitive to H ( z ). In other words, due to the potential non-linearity in the evolution of the equation of state for this model, the best fit linear CPL parameterization can look different depending on which redshifts the surveys are probing and thus which data points and observations are considered in the fitting procedure.  \nTo illustrate this point, consider the example in Fig. 2 where we depict the cosmological evolution of the equation of state for a particular (and somewhat generic) thawing model from the m 2 < 0 branch of Eq. (2) with parameters ( V 0 , m 2 ) = (0 . 9725 , -40 . 0), as well as its corresponding best fit CPL models for various data sets, echoing what we have said in [11]. In this article, we are considering data from BAO, SNe, and CMB observations, all of which can be understood to be sensitive to H ( z ) at different redshift epochs. The best fit ( w 0 , w a ) can differ quite substantially depending on the redshifts the observables are sensitive to, with the Pantheon+ SNe [3] data probing redshifts up to z ≃ 0 . 07 giving ( w 0 , w a ) ≃ ( -0 . 91, -0 . 40), DESI BAO data probing redshifts up to z ≃ 0 . 295 [1] giving ( w 0 , w a ) ≃ ( -0 . 93, -0 . 17), CMB data locked in at the photon decoupling epoch z ≃ 1090 [5] giving ( w 0 , w a ) ≃ ( -0 . 96, -0 . 04), and the combined use of all the data sets giving ( w 0 , w a ) ≃ ( -0 . 93, -0 . 15).  \nThe more recent the epoch probed, the steeper the region inferred in the ( w 0 , w a ) plane, which is due to the sharp time variation in w ( a ) at more recent epochs that this particular model exhibits. Models given by the other ( m 2 > 0) branch of Eq. (2) do not inherit this ambiguity to the same degree because they evolve far more linearly, and are thus not quite as sensitive to exactly which redshifts are probed by a survey. Furthermore, the inferred values of the CPL parameters seemingly imply that this dark energy model becomes phantom in the past. However, one can clearly see that this is just an artifact of extrapolating the linear parameterization of w ( a ) outside the range where the relevant distance observations are fit as the actual equation of state remains non-phantom across all of cosmic history (a point noted by [18] as well). Notice also that the fitted w 0 can differ quite substantially from the true value of the equation of state today, w ( a = 1).  \nTo identify the region of the ( w 0 , w a ) plane which is occupied by thawing quintessence models, we use the Boltzmann solver hi class [68-70] to generate cosmological  \nFIG. 3. 68% and 95% C.L. posterior distribution of the CPL parameters (dashed line) from the BAO (upper left), BAO+SNe (lower left), BAO+CMB (upper right), and BAO+CMB+SNe (lower right) data, where we use the data sets coming from DESI data [20], Pantheon+ [3], and Planck+ACT lensing [5-8]. Overlayed are the CPL parameters for thawing quintessence, obtained fitting the surveys characteristics as explained in Section III. One can see that there is reasonable overlap when only the DESI or DESI+SNe data is considered. Once the CMB is included and all of the data is taken together, thawing quintessence is barely viable.  \n<!-- image -->  \n0  \nmodels with dark energy described by the potential given in Eq. (2) by sampling over m 2 , V 0 , Ω m , and h (note that φ ini is fixed by imposing the Friedmann equation). We then fit the predictions of these models to the CPL parameterization in Eq. (1) using the combined data from all of these probes, and following the procedure described in Appendix A.  \nWe sample 100,000 points in latin hypercube in the ranges m 2 ∈ [ -150 . 0 , 5], V 0 ∈ [0 . 0 , 2 . 0], Ω m ∈ [0 . 29 , 0 . 33] and h ∈ [0 . 65 , 0 . 70]. Given both Ω m and h are highly constrained by the combined cosmological data we are using, these ranges encompass the allowed regions for these parameters. While we will fully broaden and open these priors when we proceed to the MCMC sampling, these tighter constraints will serve us for the present section.  \nWhen it comes to the range of mass and energy scales considered, the reasoning is as follows. For m 2 > 0, as previously noted, these models converge on a universal evolution for w ( a ) (and cosmological behavior more gen-  \nerally); therefore, we only need a very limited range of mass values in the positive branch to completely capture their phenomenology. Furthermore, at large m 2 > 0 these models enter the oscillatory regime, which would render them unviable candidates for dark energy. Thus, the cutoff at m 2 = 5.  \nThe m 2 < 0 regime is, in principle, unbounded in terms of the range of masses we could look at. This is because, as mentioned above, when one increases the mass, one can always fine-tune the initial condition φ ini such that the field starts closer to the top of the hill. However, as discussed in [71], there is a plausible effective field theory argument for cutting off the mass range at a certain point. We will be overly generous and consider a lower bound of -150 < m 2 but, in practice, we will see that, given the characteristics of the surveys, we are insensitive to differences in potentials with m 2 < -few × 10.  \nWe are now ready to map thawing quintessence directly into the latest constraints on w 0 and w a . We sam-  \nw  \nFIG. 4. Dependence of the CPL representation of thawing quintessence on survey properties. All surveys have the same upper bound (solid black line) but differing lower bounds for the corresponding ( w 0 , w a ) regions.  \n<!-- image -->  \nple the dark energy model described by the potential in Eq. (2) as detailed in this section to generate the observables and then infer, for several collections of survey data, the corresponding ( w 0 , w a ) as described in Appendix A. In Fig. 3, we plot our results for a few different data choices. Here, we see how the thawing quintessence models lie in the ( w 0 , w a ) plane relative to the constraints derived from both the BAO measurements and the combination of BAO and SNe measurements. As we can see, there is no evidence for deviations from ΛCDM, and there is reasonable overlap between what one might expect from thawing quintessence and the data. The situation changes substantially if we include the CMB. The CMB substantially tightens the constraints on ancillary cosmological parameters (such as Ω m ), leading to tighter constraints on ( w 0 , w a ). Furthermore, it probes much deeper redshifts and, as we saw in Fig. 2, this pulls the region covered by thawing quintessence towards the top right hand corner of the ( w 0 , w a ) plane, away from the region constrained by the data.  \nThe situation is most dramatic if we combine all the data sets together. There we find that the constraints on the ( w 0 , w a ) plane pull away somewhat significantly from the ΛCDM point. But we also find that the region covered by thawing quintessence models is now almost completely disjoint from the region favoured by data (albeit with a small overlap). Thus, with a combination of data sets which give the tightest constraints on ( w 0 , w a ) we seem to find that thawing quintessence is only marginally consistent with the data. Our results echo those found in Ye et al. [34]. In their analysis they simplified the mapping to ( w 0 , w a ) by fitting the equation of state directly; as we have shown in this paper, and argued in [11], the mapping depends on the characteristics of the survey one is considering.  \nIn Fig. 4 we can see that representing thawing quintessence models in terms of ( w 0 , w a ) is clearly sensitive to the properties of the survey considered. They all share a common upper boundary given by the universal behavior of the m 2 > 0 models, but how far the boundary extends for the m 2 < 0 models depends on the survey properties. That is, the redshifts at which the distance measurements are taken will impact what the best fit ( w 0 , w a ) parameters are for a particular dark energy model, with the SNe data at recent times pulling the contour down into steeper regions of the plane, while CMB data from early times pulls the contours into the flatter regions of the plane.  \nIt is useful to understand the role of m 2 in the mapping between thawing quintessence and the ( w 0 , w a ) plane. We do so by focusing on the slope of quintessence models in the ( w 0 , w a ) plane, which we define to be α = w a / (1+ w 0 ). For a given survey, this slope is primarily set by m 2 , with a small variance due to the other parameters. In Fig. 5, we plot the slope as a function of m 2 for the combination of data which includes BAOs, SNe and the CMB, following the procedure of Appendix A. We can see that as m 2 becomes increasingly negative, the slope tends towards a constant, α ≃ -2 . 25 and the lower m 2 models all eventually become indistinguishable from each other. We can compress the contours from Fig. 1 into the red horizontal band in Fig. 5. Again, as we saw in the bottom right hand panel of Fig. 3, through the lens of the CPL parametrization, thawing quintessence is only marginally consistent with current data.  \nIn summary, the CPL parametrization is a somewhat flawed form of data compression, but it is useful to gain an understanding of how a particular model compares to constraints from data. By finding the ( w 0 , w a ) for thawing quintessence (for particular survey specifications) we are, to some extent, determining the prior distribution of this class of theories in the CPL parameter plane. We can then compare it to the likelihood for ( w 0 , w a ) derived from the data and see if there is significant or very little overlap. We see that, for the case in which we use a combination of BAO, SNe and CMB data there is, indeed, very little overlap between the thawing quintessence prior and the likelihood.", 'IV. CONSTRAINTS ON THAWING QUINTESSENCE': "So far we have seen that thawing quintessence is only marginally consistent with a compilation of current data through the lens of the CPL parametrization. But one can look at the problem in a different way and take thawing quintessence as the unique cause for late time acceleration and constrain its parameters. This is tantamount to considering a prior on dark energy which is that of thawing quintessence. Thus, instead of constraining ( w 0 , w a ) from the data and then comparing these constraints to what is predicted by thaw-  \nFIG. 5. Blue: The predicted slope w a / (1 + w 0 ) of thawing quintessence models in the ( w 0 , w a ) plane for survey characteristics of BAO+CMB+SNe data. Red: CPL parameterization data constraints (as Fig. 1) . The shaded areas represent 2 σ regions for both thawing quintessence and the data.  \n<!-- image -->  \ning quintessence for those parameters, one constrains the thawing quintessence parameters directly.  \nThe cosmological model one is now constraining is effectively described in terms of a set of parameters. ( V 0 , m 2 , Ω m , . . . ) and, in principle, one would like to assume uninformative, uniform priors on all these parameters. This is particularly true of m 2 as this is the parameter that is responsible for time evolution in the equation of state. One needs to be sure that the prior is such that it does not artificially enhance or suppress evidence for non-zero m 2 .  \n<!-- image -->  \nTABLE I. Prior distributions used in the cosmological parameter inference. U ( a, b ) stands for an uniform distribution in the range [ a, b ]. For Planck's nuisance parameters, we use the Cobaya default values, specified in [6].  \nU  \n-  \nThere are a few aspects of the analysis which need to be highlighted. First of all, we saw that, as m 2 becomes ever more negative, the more indistinguishable w a / ( w 0 +1) becomes between models of different mass. This can be understood from the dynamics of the scalar field described  \nω  \nFIG. 6. 68% and 95% C.L. posterior distributions for thawing quintessence from the combination of DESI BAO data [1], Pantheon+ SNe data [3], and CMB data [5-8] obtained from sampling uniformly in the parameters given by Table I. See Fig. 14, in Appendix D, the full corner plot with all cosmological parameters.  \n<!-- image -->  \nabove - for large, negative values of m 2 , the scalar field evolves very rapidly but only at very late times. So, from the point of view of the redshift ranges being covered by the surveys which may not probe such late times, the cosmological evolution of the thawing quintessence models will become indistinguishable. Indeed, we shall see that, for a fixed compilation of data, the likelihood of the thawing quintessence model will tend to a constant as m 2 →-∞ .  \nWe have already highlighted the fact there should be a degeneracy for m 2 > 0 which, in turn, may lead to an enhanced prior for positive masses. This could, potentially, be a source of concern in interpreting the posterior constraints we obtain on m 2 as 'volume effects' could possibly bias the constraints towards m 2 ≥ 0. As we will see, that is not the case - the likelihood more than compensates for the effect of the prior there and the posterior favors m 2 < 0.  \nA crucial part of our analysis is ensuring that, for m 2 < 0, the prior is flat. This is non-trivial given how one generates the observables with hi class . Given a value of m 2 , the Boltzmann-solver has to iterate to find the corresponding initial condition, φ ini , which will lead to the cosmology with the correct values of the remaining cosmological parameters. This is a time consuming procedure for each value of m 2 and highly dependent on a 'good' choice for the iterative procedure.  \nTo speed up the iterative search, we have implemented  \nFIG. 7. Minimum χ 2 = -2 ln L (where L is the likelihood) as a function of m 2 . Blue dots are the result of the numerical minimization while the red line is the best fit to the dots for m 2 < 0, shown for illustrative purposes.  \n<!-- image -->  \na novel procedure, described in Appendix B: we have used Symbolic Regression to determine an analytic expression for φ ini as a function of the remaining parameters, most notably m 2 , which is accurate to within 1%. This is then used as the initial guess for the solver and allows us to efficiently generate a flat prior for m 2 . We confirm that our procedure leads to flat priors by running the inference without any data (in practice, by setting the likelihood to a constant) as can be see in Fig. 13 in Appendix C.  \nWith the priors correctly under control, we proceed to infer the posterior on m 2 for a range of data sets. We use the same data as in the DESI analysis [1]; i.e., DESI 2024 BAO measurements [1, 72, 73], the SNe Ia Pantheon+ samples [3, 74], the TTTEEE CMB measurements from Planck 2018 [5, 6], and the CMB lensing from ACT DR6 [7, 8].  \nDESI 2024 BAO measurements consist of a set of six measurements of the angular diameter distance, D M , the Hubble parameter, D H = c/H , and the angle-averaged quantity D V = ( zD 2 M D H ) 3 , relative to the sound horizon, obtained from measurements of the BAO scale, using 6 million galaxies in the range of redshifts 0 . 1 < z < 4 . 2.  \nThe Pantheon+ sample is made of 1701 light curves from 1550 different SNe Ia ranging z = 0 . 001 to 2 . 26, obtained with observations of DES [75, 76], Foundation [77], Pan-STARRS [78], Supernova Legacy Survey [79], SDSS [80], HST [67, 81-85] and multiple lowz samples (see [3] for the full list).  \nFrom Planck, we use the auto- and cross-correlation of the CMB temperature and E -mode polarization fields. We use the lowℓ power spectra, obtained with the 'Commander' component separation algorithm, in the range 2 ≤ ℓ ≤ 29 and the highℓ plik likelihood covering 30 ≤ ℓ ≤ 2508 for the temperature auto-correlation ( TT ) and 30 ≤ ℓ ≤ 1996 for the TE and EE components [86]. We marginalize over the nuisance parameters related to the foregrounds calibration.  \nFinally, we use CMB lensing from ACT DR6 as provided by the official likelihood 1 with the variant act baseline [7, 8], consisting of the ACT CMB reconstructed lensing power spectra between the scales 40 < L < 763.  \nWe use the likelihoods as implemented in Cobaya [87, 88], which we then use to explore the parameter space using a Metropolis-Hasting algorithm [89-92] by sampling the parameters in Table I and imposing uniform priors. We address convergence using the GelmanRubin criterion [93], imposing that R -1 < 0 . 02 in the diagonalized parameter space. Finally, we compute the theory predictions with hi class and model the nonlinear matter power spectrum with the halofit algorithm [94, 95] 2 .  \nIn Fig. 6, we show the posteriors for the cosmological parameters, including ( V 0 , m 2 ) for the combination of BAOs, SN and CMB data sets. The posterior for m 2 is qualitatively, as expected, favoring m 2 < 0 and flat out to the minimum (most negative) m 2 considered. Indeed, the data favors a thawing quintessence model over the m 2 = 0, ΛCDM model. Unfortunately, given the improper nature of the prior on m 2 (i.e., it is in principle, constant for m 2 < 0 and unbounded from below), it is dependent on our choice of the lower bound, and thus impossible to make a clear (prior independent) statement on how much we can rule out ΛCDM vis-'a-vis thawing quintessence. This situation is not too dissimilar to what one find when constraining Jordan-Brans-Dicke (JBD) models [97], where one can show how constraints on the JBD parameter which characterizes the non-minimal coupling of the scalar field to gravity is completely dependent on the choice of priors.  \nIt is useful, nevertheless, to try and assess if the data has a preference for thawing quintessence relative to ΛCDM. The simplest approach is to look at the χ 2 = -2 ln ˆ L ( m 2 ), where ˆ L is the likelihood of the model with mass m 2 and maximized over the other parameters 3 . From Fig. 7, where we plot the χ 2 of the model as a function of m 2 , we can see that the χ 2 assymptotes to χ 2 f ≃ 2434 for sufficiently negative values of m 2 . If we now compare the model with m 2 = 0 (where there is no thawing) with one with a suitably negative value of m 2 , we find a ∆ χ 2 ≃ 2 . 8 showing a very marginal preference for thawing quintessence.  \nWe can further refine the comparison by considering the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) as these incorporate more information about the model [98-100]. The AIC is given by AIC = 2 k + χ 2 , where k is the number of parameters  \nin the model. We can now compare the ΛCDM model with a thawing model with negative m 2 . Apart from the cosmological parameters, ΛCDM has { V 0 } , whereas thawing quintessence has { V 0 , m 2 , φ ini } . In the case of thawing quintessence, φ ini (or any of the others) is fixed by the Friedmann equation. In the case of ΛCDM, this is also true: V 0 is fixed by the Friedmann equation. As a consequence, ∆ k ≡ k ΛCDM -k Thawing = -2 and ∆AIC ≃ -1 . 2, rendering them effectively (statistically) indistinguishable, albeit a mild preference towards ΛCDM. If we consider the BIC, which is given by BIC = k ln( n ) + χ 2 , where n ≃ 2400 is the number of data points, we have that ∆BIC ≃ -12 . 8 which strongly favours the ΛCDM model. Thus, we can conclude that there is little evidence for thawing quintessence with current cosmological data.  \nWe can compare these results with what we find if, instead, we look at the CPL parametrization. There is a more marked difference between the best fit ( w 0 , w a ) model as compared to ΛCDM. There we find that ∆ χ 2 ≃ 6 . 8, ∆AIC ≃ 2 . 8 and ∆BIC ≃ -8 . 8 indicating that, from the point of view of two of our measures, the best fit CPL model is slightly preferred while, again, from the point of view of the BIC, it is not.  \nFinally, we take a look again at the equation of state for thawing quintessence, given the current data. In Fig. 8, we show constraints on the value of w ( z ) across the recent expansion history going out to z ≃ 0 . 5 assuming that dark energy is described by the present model of thawing quintessence. This shows that, until quite recently, if dark energy is described by thawing quintessence, w ( a ) would be indistinguishable from a cosmological constant, before becoming quite unconstrained at recent times. Even here though, from the perspective of constraining directly in terms of ( V 0 , m 2 ), a value for the equation of state indistinguishable from ΛCDM is still allowed. It is also interesting to see, in Fig. 9, constraints on the value of the equation of state w ( a = 1) today in terms of m 2 . Unsurprisingly, the larger the negative hilltop mass becomes, the more unconstrained this value becomes as the equation of state can evolve very fast at these masses.", 'V. CONCLUSION': "In this paper we have looked at the evidence for evolving dark energy in the form of quintessence, i.e., a scalar which is minimally coupled to gravity and has a canonical kinetic energy. Current constraints point, if anything, towards thawing quintessence, where the scalar field starts to evolve at late times. Following [11], we consider the potential which encompasses all possible thawing quintessence-like behaviour.  \nFirst, we transform its equation of state w ( a ) into the CPL parametrization, ( w 0 , w a ), and compare its constraints with current data. We find that the evidence for thawing quintessence, from the point of view of CPL, is marginal at best, in agreement with what was found in Ye  \nFIG. 8. Reconstruction of the equation of state (1- and 2σ levels) from the combination of DESI BAO data [20], Pantheon+ SNe data [3], and CMB data [5-8]. The equation of state is indistinguishable from a cosmological constant for most of the Universe's history, becoming unconstrained at recent times.  \n<!-- image -->  \nFIG. 9. 68% and 95% C.L. posterior distribution of w ( a = 1) (not to be confused with the CPL parameter w 0 ) and m 2 in thawing quintessence from the combination of DESI BAO data [20], Pantheon+ SNe data [3], and CMB data [5-8].  \n<!-- image -->  \net al. [34]. An essential aspect of our analysis is that we determine the ( w 0 , w a ) parameters of the quintessence model using exactly the same procedure as was used to obtain them from the data. This means we are comparing like with like and highlights the fact that a given quintessence model does not map on to a unique value of ( w 0 , w a ).  \nWe then directly constrain the parameters of thawing quintessence, ( V 0 , m 2 ), and find there is some evidence for evolution, distinguishing it from a cosmological constant. The posterior distribution of m 2 is such that it is difficult to quantify how strong the evidence is. Nevertheless, by looking at the relative likelihoods of the best fit thawing quintessence model relative to ΛCDM, and using other information criteria, we find that thawing quintessence is between mildly (∆AIC ≃ -1 . 2) and strongly (∆BIC ≃ -12 . 8) disfavored with respect to ΛCDM. In other words, there is no preference for thawing quintessence.  \nAs we have argued, the quadratic potential we are considering here covers, in our view, a vast range of thawing quintessence models. As shown in Wolf and Ferreira [11], this means that we can span a large wedge in the ( w ( a = 1) , -w ' ( a = 1)) parameters. But this is not the same as the ( w 0 , w a ) plane, as explained above. One might, however, consider thawing like equations of state which deviate from the ones considered here and that could have a more significant overlap with the ( w 0 , w a ) constraints. Qualitatively they would have to behave as follows: the equation of state would start off at -1, then, at some intermediate redshift (coincident with the redshift ranges of the data we are considering here) it would thaw rapidly to then stabilize again at some constant value until it reaches today. Even so, it would be difficult to obtain the slopes detected in the ( w 0 , w a ) given the role of the CMB data (see Fig 2). Furthermore, it would require some creative engineering of a potential that would lead to that behaviour.  \nThe fact that thawing quintessence does so poorly relative to the view of CPL parameterization for dark energy is intriguing. Does it mean we need to consider a broader family of scalar field models, with non-minimal couplings or non-canonical kinetic energy? This might be the case, although as argued above and by others, the evidence for phantom evolution from ( w 0 , w a ) is misleading. On the other hand, Ye et al. [34] have come up with a proof of concept that it is possible to construct a (simple) well behaved non-minimally coupled model that leads to phantom behaviour and is consistent with the data. Dubbed 'Thawing Gravity' it will be interesting to see if it is consistent with other constraints on fifth forces which arise in non-minimally coupled theories. At the very least, it raises the possibility that we are detecting evidence of non-minimal coupling on cosmological scales.  \nAs stated in the introduction, the CPL model is a blunt instrument and clearly too blunt to make any strong inferences about the underlying microphysics (or even just the qualitative behavior) of dark energy. Its strength comes from the fact that it is very simple as it only depends on two parameters. Other proposals have been put forward which are as minimal but run into the same problems as the CPL parameters. On the other hand, there are a number of more complicated parameterizations which track the evolution of dark energy more accurately, but these run the risk of overfitting and be-  \ning plagued by degeneracies between the parameters. Clearly, an all purpose, simple yet accurate parameterization is still lacking.  \nThe analysis in this paper hinges on the data being correct; this might well not be the case as pointed in [12, 13]. It would be wise to reserve judgement until there are more, and better, constraints on the expansion history of the Universe. We expect this to be the case over the coming years.", 'ACKNOWLEDGEMENTS': "We acknowledge the hi class developers, Emilio Bellini, Miguel Zumalac'arregui and Ignacy Sawicki, for sharing a private version of the code. We thank David Alonso, Harry Desmond, Constantinos Skordis and Maria Vincenzi for useful discussions. WW is supported by St. Cross College, Oxford. CGG is supported by the Beecroft Trust. DJB is supported by the Simons Collaboration on 'Learning the Universe.' PGF is supported by STFC and the Beecroft Trust.  \nSoftware : We made extensive use of hi class [68-70], Cobaya [87, 88] and the numpy [101, 102], scipy [103], matplotlib [104], GetDist [105] and pyoperon [106] python packages.  \nFor the purposes of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.", 'Appendix A: Compressing the data.': 'To account for the surveys constraining power when estimating ( w 0 , w a ), for simplicity, we compress the CMB and SNe samples into cosmological background quantities that we can then fit following Section III.', '1. CMB': "We compress the CMB data (Planck's TTTEEE [5, 6] and ACTDR6 lensing [7, 8]) into the background quantities [64-66]: the angular scale of the sound horizon  \nℓ A = π D A ( z ∗ ) r s ( z ∗ ) = π θ ∗ , (A1)  \nand the CMB shift parameter  \nR ( z ∗ ) = √ Ω m H 2 0 D A ( z ∗ ) , (A2)  \nat last scattering (i.e. the photon decoupling redshift z ∗ ). In addition, one also needs to add constraints on ω b = Ω b h 2 and n s to recover the full likelihoods constraints. We show in Fig. 10 that a Gaussian likelihood using these four parameters, with mean and covariance shown in Table II, is good enough to reproduce the CMB posterior distributions.  \nFIG. 10. Posterior distribution (68% and 95% C.L.) obtained with the full CMB likelihood (blue) in comparison with the compressed version (red) for the ΛCDM model. The compressed CMB data recovers the posterior.  \n<!-- image -->", '2. SNe IA': 'We compress the SNe Ia data following [67]. Briefly, the SNe data measures  \nD L = 1 + z H 0 ∫ d z E ( z ) , (A3)  \nwhere we have defined the expansion rate E ( z ) ≡ H ( z ) /H 0 . One can compress the SNe data in model independent measurements of 1 /E ( z ): instead of obtaining 1 /E ( z i ) from a cosmological model we can directly compute it as a spline: 1 /E ( z ) = spline(1 /E ( z i )), with nodes 1 /E ( z i ). We use a cubic spline with nodes defined at z = { 0 . 07 , 0 . 2 , 0 . 35 , 0 . 55 , 0 . 9 , 1 . 5 , 2 , 2 . 5 } , which are the same as in [67], extended to higher redshift due to the extra SNe in the Pantheon+ sample. Pantheon+ is not able to constrain the highest redshift, so we drop that node from the compressed data set. Furthermore, we approximate the compressed likelihood as a Gaussian distribution with the mean and covariance shown in Table II. The theory code to compress the data and the Cobaya likelihood can be found in https: //gitlab.com/carlosggarcia/mcmc-likelihoods . Finally, we show in Fig. 11 that we can accurately recover the posterior distribution with the compressed likelihood.  \nCombining the results of our data compression for both the CMB and SNe IA, we show in Fig. 12 that the compressed data can successfully recover the headline results from [1] and depicted in Fig. 1. That is, when combined  \nFIG. 11. Posterior distribution (68% and 95% C.L.) obtained with the Pantheon+ SNe Ia likelihoods and its compressed version. One can see the later agrees well with the full likelihood.  \n<!-- image -->  \nwith DESI BAO data, the compressed CMB and SNe IA variables reproduce the results of the full data to a very good approximation.', 'Appendix B: Analytic approximation for initial conditions': 'The history of a thawing quintessence model is uniquely determined by the parameters of the potential, ( V 0 , m 2 ) and the initial value of the scalar field φ ini (recall that we are assuming that ˙ φ ini = 0) at some early time. This will determine the final fractional energy density in the scalar field Ω φ = 1 -Ω m . We are interested, however, in choosing Ω φ , or more specifically, Ω m . This means that a set of ( V 0 , m 2 ) and Ω m should uniquely determine φ ini . In practice, this is done by using a rootfinding algorithm in hi class , which finds the correct φ ini . This is a time consuming process and sufficiently onerous that it makes an MCMC inference prohibitively slow. It is preferable to have an accurate predictor of φ ini which can greatly reduce evaluations of hi class .  \nFor both computational efficiency and for interpretability, we wish to find symbolic approximations for φ ini as a function of V , m 2 , and the background cosmological parameters Ω m and h . To obtain such an expression, we utilize the supervized machine learning technique of Symbolic Regression (SR) [107]: an automated  \n±  \nTABLE II. Left. Pantheon+ SNe data compressed in data driven measurements of 1 /E ( z ). Right. CMB compressed variables. We show in Fig. 10 and Fig. 11 that a Gaussian likelihood with the mean and covariance provided here is good enough to recover the posterior distribution, in both cases.  \n<!-- image -->  \nw  \nFIG. 12. Posterior distribution (68% and 95% C.L.) obtained with DESI BAO and the full CMB and Pantheon+ SNe likelihoods (blue), compared to that obtained with DESI BAO and the compressed likelihoods of the CMB and SNe. The compressed data, combined with DESI, can clearly recover the constraints on ( w 0 , w a ) that come from the full complement of data to a high degree of accuracy.  \nsearch through function space to find simple mathematical expressions which well approximate a given dataset.  \nWe use the genetic programming based SR code operon [106], which uses natural selection inspired processes such as crossover (breeding) and mutation to evolve an initial population of functions over time, where well-fitting and simple expressions are preferentially kept and worse-fitting or overly-complex functions are discarded. As such, after several iterations, a population of expressions which well approximate the dataset emerge. We choose operon over other SR codes due to its strong performance in benchmark tests [108, 109] and as it has  \nalready been demonstrated to be effective in cosmological and astrophysical studies [110-113].  \nTo obtain a symbolic approximation, we first need a training set. It is easier to obtain m 2 as a function of φ ini and the other parameters than vice versa, hence we generated ∼ 10,000 points on a Latin hypercube in the range φ ini ∈ [0 . 1 , 10], V 0 ∈ [ -10 , 10], Ω m ∈ [0 . 01 , . 99] and h ∈ [0 . 2 , 1 . 0]. We then obtain the corresponding value of m 2 for each set of parameters by using hi class . Any points which do not provide converged solutions are removed. We use 90% of these data for training and the remainder as a validation set.  \nDue to the large dynamic range of φ ini , we choose to fit for log φ ini instead of φ ini directly. We apply a multiobjective optimization scheme, where we attempt to minimize both the mean absolute error of log φ ini (and thus the fractional error on φ ini ) and the model length, which is approximately equal to the number of symbols appearing in the expression. operon relies on the concept of ϵ -dominance [114], such that these objective values are considered equal in the search if they fall within some hyper-parameter ϵ of each other. After some experimentation, we find that ϵ = 10 -2 is an appropriate choice for balancing accuracy and simplicity.  \nWe fit log φ ini as a function of V , m 2 , Ω m and h , where we allow symbolic expressions containing the operators + , -, × , /, √ · , pow , log , exp and free parameters, which are optimized [115] using the Levenberg-Marquardt algorithm [116, 117]. We run operon for up to 24 hours, allowing a maximum model length of 140. We visually inspect the candidate solutions to obtain one which we deem to balance accuracy and simplicity, enforcing that our solution has a mean absolute error of at least 1% on both the training and validation sets. These models are often over-parameterized, so we manually remove superfluous parameters and attempt to round those which are close to integers, zero or simple rational numbers to try to simplify the expression.  \nGiven the very different behaviors of the m 2 < 0 and m 2 > 0 branches, we choose to use a different expression for each. After some manipulation, for m 2 < 0, we obtain  \n±  \nthe following expression of length 23 for the initial φ :  \nlog φ ini ≈ m 2 e b 0 h + 1 2 log ( b 1 h 2 (1 -Ω m ) -V 0 m 2 ) , (B1)  \nwhere b 0 ≈ -1 . 5 and b 1 ≈ 3 . 0. This approximation gives  \nan error on the training and test sets of 8 × 10 -3 and 9 × 10 -3 , respectively, so is sufficiently accurate for our purposes.  \nWe find that the functional form for m 2 > 0 has to be more complex, with an expression of length 48 fulfilling our selection criteria. This is given by  \nlog φ ini ≈ c 0 + c 1 ( Ω m c 2 V 0 + c 3 h + e c 4 Ω m + c 5 V 0 ) -1 + c 6 ( Ω m c 7 V 0 + c 8 h + e c 9 Ω m -c 10 m 2 + c 11 ( c 12 Ω m -m 2 ) c 13 ( Ω m + c 14 V 0 -c 15 m 2 ) ) ( ( c 16 Ω m ) c 17 h c 18 -e c 19 V 0 ) -1 , (B2)  \nwhere c = { 0 . 006 , 0 . 0558 , 0 . 8222 , -3 . 1121 , 1 . 4964 , -1 . 2683 , -1 . 775 , 0 . 9655 , -3 . 3049 , 1 . 3418 , -0 . 0816 , -1 . 1722 , 0 . 0291 , -0 . 3267 , 0 . 2304 , -0 . 9416 , 54 . 7118 , 0 . 3744 , 0 . 6663 , -0 . 5982 } . We find that this expression gives a mean absolute error of 0.01 on the training set, and 8 × 10 -3 on the test set, so is able to give a percentlevel approximation for φ ini , so is again sufficiently accurate for our purposes. The small errors in our approximations are corrected for by the root-finding algorithm in hi class , but this requires fewer iterations than if we used a random starting point for the optimizer.', 'Appendix C: m 2 prior': 'As discussed in Section II, an important aspect of our work is ensuring that we have effective flat priors in the parameters that we sample in the inference. This is particularly relevant for m 2 if we want to compare in a meaningful way thawing quintessence with ΛCDM. We show in Fig. 13 the effective prior and the 1D-posterior distribution from the combination of all data (DESI BAO, Pantheon+ SNe Ia and CMB). As mentioned in Section II, there is a degeneracy between m 2 and V 0 that enhances the prior on the m 2 > 0 branch slightly. However, the prior on m 2 is effectively flat across the range of masses considered. We find that this extra volume at positive m 2 does not drive the constraints and data actually disfavors this regime.', 'Appendix D: Full posterior distribution for thawing quintessence': "Fig. 14 gives the full corner plot featuring all parameters varied as per Table. I when constraining thawing quintessence in terms of the parameters ( V 0 , m 2 ).  \nFIG. 13. 1D-posterior distribution for m 2 from the combination of DESI BAO data [1], Pantheon+ SNe data [3], and CMB data [5-8] (blue) and the effective prior imposed (gray).  \n<!-- image -->  \nm  \nω  \nFIG. 14. 68% and 95% C.L. marginalized posterior distributions for the full cosmological parameter space of thawing quintessence from the combination of DESI BAO data [1], Pantheon+ SNe data [3], and CMB data [5-8].  \n<!-- image -->  \n- [1] A. G. Adame et al. (DESI), DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations, (2024), arXiv:2404.03002 [astro-ph.CO].\n- [2] D. Rubin et al. , Union Through UNITY: Cosmology with 2,000 SNe Using a Unified Bayesian Framework, (2023), arXiv:2311.12098 [astro-ph.CO].\n- [3] D. Scolnic et al. , The Pantheon+ Analysis: The Full Data Set and Light-curve Release, Astrophys. J. 938 , 113 (2022), arXiv:2112.03863 [astro-ph.CO].\n- [4] T. M. C. Abbott et al. (DES), The Dark Energy Survey: Cosmology Results With ˜1500 New High-redshift Type Ia Supernovae Using The Full 5-year Dataset, (2024), arXiv:2401.02929 [astro-ph.CO].\n- [5] N. Aghanim et al. (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].\n- [6] N. Aghanim et al. (Planck), Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641 , A5 (2020), arXiv:1907.12875 [astro-ph.CO].\n- [7] F. J. Qu et al. (ACT), The Atacama Cosmology Telescope: A Measurement of the DR6 CMB Lensing Power Spectrum and Its Implications for Structure Growth, Astrophys. J. 962 , 112 (2024), arXiv:2304.05202 [astroph.CO].\n- [8] M. S. Madhavacheril et al. (ACT), The Atacama Cosmology Telescope: DR6 Gravitational Lensing Map and Cosmological Parameters, Astrophys. J. 962 , 113 (2024), arXiv:2304.05203 [astro-ph.CO].\n- [9] E. V. Linder, Exploring the expansion history of the universe, Phys. Rev. Lett. 90 , 091301 (2003), arXiv:astroph/0208512.\n- [10] M. Chevallier and D. Polarski, Accelerating universes with scaling dark matter, Int. J. Mod. Phys. D 10 , 213 (2001), arXiv:gr-qc/0009008.\n- [11] W. J. Wolf and P. G. Ferreira, Underdetermination of dark energy, Phys. Rev. D 108 , 103519 (2023), arXiv:2310.07482 [astro-ph.CO].\n- [12] G. Efstathiou, Challenges to the Lambda CDM Cosmology (2024) arXiv:2406.12106 [astro-ph.CO].\n- [13] G. Efstathiou, Evolving dark energy or supernovae systematics? (2024), arXiv:2408.07175 [astro-ph.CO].\n- [14] Y. Tada and T. Terada, Quintessential interpretation of the evolving dark energy in light of DESI observations, Phys. Rev. D 109 , L121305 (2024), arXiv:2404.05722 [astro-ph.CO].\n- [15] K. V. Berghaus, J. A. Kable, and V. Miranda, Quantifying Scalar Field Dynamics with DESI 2024 Y1 BAO measurements, (2024), arXiv:2404.14341 [astroph.CO].\n- [16] C.-G. Park, J. de Cruz Perez, and B. Ratra, Using nonDESI data to confirm and strengthen the DESI 2024 spatially-flat w 0 w a CDM cosmological parameterization result, (2024), arXiv:2405.00502 [astro-ph.CO].\n- [17] W. Yin, Cosmic Clues: DESI, Dark Energy, and the Cosmological Constant Problem, (2024), arXiv:2404.06444 [hep-ph].\n- [18] D. Shlivko and P. J. Steinhardt, Assessing observational constraints on dark energy, Phys. Lett. B 855 , 138826 (2024), arXiv:2405.03933 [astro-ph.CO].  \n- [19] M. Cortˆes and A. R. Liddle, Interpreting DESI's evidence for evolving dark energy, (2024), arXiv:2404.08056 [astro-ph.CO].\n- [20] K. Lodha et al. (DESI), DESI 2024: Constraints on Physics-Focused Aspects of Dark Energy using DESI DR1 BAO Data, (2024), arXiv:2405.13588 [astroph.CO].\n- [21] B. R. Dinda, A new diagnostic for the null test of dynamical dark energy in light of DESI 2024 and other BAO data, (2024), arXiv:2405.06618 [astro-ph.CO].\n- [22] Y. Carloni, O. Luongo, and M. Muccino, Does dark energy really revive using DESI 2024 data?, (2024), arXiv:2404.12068 [astro-ph.CO].\n- [23] D. Wang, The Self-Consistency of DESI Analysis and Comment on 'Does DESI 2024 Confirm ΛCDM?', (2024), arXiv:2404.13833 [astro-ph.CO].\n- [24] K. S. Croker, G. Tarl'e, S. P. Ahlen, B. G. Cartwright, D. Farrah, N. Fernandez, and R. A. Windhorst, DESI Dark Energy Time Evolution is Recovered by Cosmologically Coupled Black Holes, (2024), arXiv:2405.12282 [astro-ph.CO].\n- [25] P. Mukherjee and A. A. Sen, Model-independent cosmological inference post DESI DR1 BAO measurements, (2024), arXiv:2405.19178 [astro-ph.CO].\n- [26] N. Roy, Dynamical dark energy in the light of DESI 2024 data, (2024), arXiv:2406.00634 [astro-ph.CO].\n- [27] H. Wang and Y.-S. Piao, Dark energy in light of recent DESI BAO and Hubble tension, (2024), arXiv:2404.18579 [astro-ph.CO].\n- [28] H. Wang, Z.-Y. Peng, and Y.-S. Piao, Can recent DESI BAO measurements accommodate a negative cosmological constant?, (2024), arXiv:2406.03395 [astro-ph.CO].\n- [29] I. D. Gialamas, G. Hutsi, K. Kannike, A. Racioppi, M. Raidal, M. Vasar, and H. Veermae, Interpreting DESI 2024 BAO: late-time dynamical dark energy or a local effect?, (2024), arXiv:2406.07533 [astro-ph.CO].\n- [30] A. Notari, M. Redi, and A. Tesi, Consistent Theories for the DESI dark energy fit, (2024), arXiv:2406.08459 [astro-ph.CO].\n- [31] L. Orchard and V. H. C'ardenas, Probing Dark Energy Evolution Post-DESI 2024, (2024), arXiv:2407.05579 [astro-ph.CO].\n- [32] H. Wang, G. Ye, and Y.-S. Piao, Impact of evolving dark energy on the search for primordial gravitational waves, (2024), arXiv:2407.11263 [astro-ph.CO].\n- [33] W. Giar'e, M. Najafi, S. Pan, E. Di Valentino, and J. T. Firouzjaee, Robust Preference for Dynamical Dark Energy in DESI BAO and SN Measurements, (2024), arXiv:2407.16689 [astro-ph.CO].\n- [34] G. Ye, M. Martinelli, B. Hu, and A. Silvestri, Nonminimally coupled gravity as a physically viable fit to DESI 2024 BAO, (2024), arXiv:2407.15832 [astroph.CO].\n- [35] B. R. Dinda and R. Maartens, Model-agnostic assessment of dark energy after DESI DR1 BAO, (2024), arXiv:2407.17252 [astro-ph.CO].\n- [36] J.-Q. Jiang, D. Pedrotti, S. S. da Costa, and S. Vagnozzi, Non-parametric late-time expansion history reconstruction and implications for the Hubble tension in light of DESI, (2024), arXiv:2408.02365 [astroph.CO].\n- [37] B. Ghosh and C. Bengaly, Consistency tests between SDSS and DESI BAO measurements, (2024), arXiv:2408.04432 [astro-ph.CO].\n- [38] O. Luongo and M. Muccino, Model independent cosmographic constraints from DESI 2024, (2024), arXiv:2404.07070 [astro-ph.CO].\n- [39] A. C. Alfano, O. Luongo, and M. Muccino, Cosmological constraints from calibrated E p -E iso gamma-ray burst correlation by using DESI 2024 data release, (2024), arXiv:2408.02536 [astro-ph.CO].\n- [40] J. a. Rebou¸cas, D. H. F. de Souza, K. Zhong, V. Miranda, and R. Rosenfeld, Investigating Late-Time Dark Energy and Massive Neutrinos in Light of DESI Y1 BAO, (2024), arXiv:2408.14628 [astro-ph.CO].\n- [41] Y.-H. Pang, X. Zhang, and Q.-G. Huang, Constraints on Redshift-Binned Dark Energy using DESI BAO Data, (2024), arXiv:2408.14787 [astro-ph.CO].\n- [42] O. F. Ramadan, J. Sakstein, and D. Rubin, DESI Constraints on Exponential Quintessence, (2024), arXiv:2405.18747 [astro-ph.CO].\n- [43] J. J. Heckman, O. F. Ramadan, and J. Sakstein, First Constraints on a Pixelated Universe in Light of DESI, (2024), arXiv:2406.04408 [astro-ph.CO].\n- [44] S. Bhattacharya, G. Borghetto, A. Malhotra, S. Parameswaran, G. Tasinato, and I. Zavala, Cosmological constraints on curved quintessence, (2024), arXiv:2405.17396 [astro-ph.CO].\n- [45] R. Calderon et al. (DESI), DESI 2024: Reconstructing Dark Energy using Crossing Statistics with DESI DR1 BAO data, (2024), arXiv:2405.04216 [astro-ph.CO].\n- [46] W. J. Wolf, P. G. Ferreira, and C. Garc'ıaGarc'ıa, Matching current observational constraints with non-minimally coupled dark energy, (2024), arXiv:2409.17019 [astro-ph.CO].\n- [47] R. Arjona and S. Nesseris, A swampland conjecture DESIder'atum?, (2024), arXiv:2409.14990 [astroph.CO].\n- [48] R. R. Caldwell and E. V. Linder, The Limits of quintessence, Phys. Rev. Lett. 95 , 141301 (2005), arXiv:astro-ph/0505494.\n- [49] P. G. Ferreira and M. Joyce, Cosmology with a primordial scaling field, Phys. Rev. D 58 , 023503 (1998), arXiv:astro-ph/9711102.\n- [50] E. J. Copeland, A. R. Liddle, and D. Wands, Exponential potentials and cosmological scaling solutions, Phys. Rev. D 57 , 4686 (1998), arXiv:gr-qc/9711068.\n- [51] R. R. Caldwell, R. Dave, and P. J. Steinhardt, Cosmological imprint of an energy component with general equation of state, Phys. Rev. Lett. 80 , 1582 (1998), arXiv:astro-ph/9708069.\n- [52] J. A. Frieman, C. T. Hill, A. Stebbins, and I. Waga, Cosmology with ultralight pseudo Nambu-Goldstone bosons, Phys. Rev. Lett. 75 , 2077 (1995), arXiv:astroph/9505060.\n- [53] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 , 1 (2016), arXiv:1510.07633 [astro-ph.CO].\n- [54] S. Dutta and R. J. Scherrer, Hilltop Quintessence, Phys. Rev. D 78 , 123525 (2008), arXiv:0809.4441 [astro-ph].\n- [55] L. Boubekeur and D. H. Lyth, Hilltop inflation, JCAP 07 , 010, arXiv:hep-ph/0502047.\n- [56] W. J. Wolf, Minimizing the tensor-to-scalar ratio in single-field inflation models, Phys. Rev. D 110 , 043521 (2024), arXiv:2407.00358 [astro-ph.CO].  \n- [57] C. Garc'ıa-Garc'ıa, E. Bellini, P. G. Ferreira, D. Traykova, and M. Zumalac'arregui, Theoretical priors in scalar-tensor cosmologies: Thawing quintessence, Phys. Rev. D 101 , 063508 (2020), arXiv:1911.02868 [astro-ph.CO].\n- [58] R. J. Scherrer and A. A. Sen, Thawing quintessence with a nearly flat potential, Phys. Rev. D 77 , 083515 (2008), arXiv:0712.3450 [astro-ph].\n- [59] E. V. Linder, Quintessence's last stand?, Phys. Rev. D 91 , 063006 (2015), arXiv:1501.01634 [astro-ph.CO].\n- [60] R. J. Scherrer, Mapping the Chevallier-Polarski-Linder parametrization onto Physical Dark Energy Models, Phys. Rev. D 92 , 043001 (2015), arXiv:1505.05781 [astro-ph.CO].\n- [61] D. J. E. Marsh, P. Bull, P. G. Ferreira, and A. Pontzen, Quintessence in a quandary: Prior dependence in dark energy models, Phys. Rev. D 90 , 105023 (2014), arXiv:1406.2301 [astro-ph.CO].\n- [62] T. Chiba, Slow-Roll Thawing Quintessence, Phys. Rev. D 79 , 083517 (2009), [Erratum: Phys.Rev.D 80, 109902 (2009)], arXiv:0902.4037 [astro-ph.CO].\n- [63] M. S. Turner, Coherent Scalar Field Oscillations in an Expanding Universe, Phys. Rev. D 28 , 1243 (1983).\n- [64] E. Komatsu, J. Dunkley, M. R. Nolta, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. Limon, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut, S. S. Meyer, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, Five-Year Wilkinson Microwave Anisotropy Probe Observations: Cosmological Interpretation, ApJS 180 , 330 (2009), arXiv:0803.0547 [astroph].\n- [65] P. A. R. Ade et al. (Planck), Planck 2015 results. XIV. Dark energy and modified gravity, Astron. Astrophys. 594 , A14 (2016), arXiv:1502.01590 [astro-ph.CO].\n- [66] L. Chen, Q.-G. Huang, and K. Wang, Distance Priors from Planck Final Release, JCAP 02 , 028, arXiv:1808.05724 [astro-ph.CO].\n- [67] A. G. Riess et al. , Type Ia Supernova Distances at Redshift > 1.5 from the Hubble Space Telescope Multi-cycle Treasury Programs: The Early Expansion Rate, Astrophys. J. 853 , 126 (2018), arXiv:1710.00844 [astroph.CO].\n- [68] M. Zumalac'arregui, E. Bellini, I. Sawicki, J. Lesgourgues, and P. G. Ferreira, hi class: Horndeski in the Cosmic Linear Anisotropy Solving System, J. Cosmology Astropart. Phys. 2017 , 019 (2017), arXiv:1605.06102 [astro-ph.CO].\n- [69] E. Bellini, I. Sawicki, and M. Zumalac'arregui, hi class background evolution, initial conditions and approximation schemes, J. Cosmology Astropart. Phys. 2020 , 008 (2020), arXiv:1909.01828 [astro-ph.CO].\n- [70] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes, J. Cosmology Astropart. Phys. 2011 , 034 (2011), arXiv:1104.2933 [astro-ph.CO].\n- [71] S. Tsujikawa, Quintessence: A Review, Class. Quant. Grav. 30 , 214003 (2013), arXiv:1304.1961 [gr-qc].\n- [72] A. G. Adame et al. (DESI), DESI 2024 III: Baryon Acoustic Oscillations from Galaxies and Quasars, (2024), arXiv:2404.03000 [astro-ph.CO].\n- [73] A. G. Adame et al. (DESI), DESI 2024 IV: Baryon Acoustic Oscillations from the Lyman Alpha Forest, (2024), arXiv:2404.03001 [astro-ph.CO].\n- [74] D. Brout et al. , The Pantheon+ Analysis: Cosmological Constraints, Astrophys. J. 938 , 110 (2022), arXiv:2202.04077 [astro-ph.CO].\n- [75] D. Brout et al. (DES), First Cosmology Results Using Type Ia Supernovae From the Dark Energy Survey: Photometric Pipeline and Light Curve Data Release, Astrophys. J. 874 , 106 (2019), arXiv:1811.02378 [astroph.IM].\n- [76] C. B. D'Andrea et al. (DES), First Cosmology Results using Supernovae Ia from the Dark Energy Survey: Survey Overview, Performance, and Supernova Spectroscopy, Astron. J. 160 , 267 (2020), arXiv:1811.09565 [astro-ph.CO].\n- [77] R. J. Foley et al. , The Foundation Supernova Survey: Motivation, Design, Implementation, and First Data Release, Mon. Not. Roy. Astron. Soc. 475 , 193 (2018), arXiv:1711.02474 [astro-ph.HE].\n- [78] D. M. Scolnic et al. (Pan-STARRS1), The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample, Astrophys. J. 859 , 101 (2018), arXiv:1710.00845 [astro-ph.CO].\n- [79] M. Betoule et al. (SDSS), Improved cosmological constraints from a joint analysis of the SDSS-II and SNLS supernova samples, Astron. Astrophys. 568 , A22 (2014), arXiv:1401.4064 [astro-ph.CO].\n- [80] M. Sako et al. , Photometric Type IA Supernova Candidates from the Three-Year SDSS-II SN Survey Data, Astrophys. J. 738 , 162 (2011), arXiv:1107.5106 [astroph.CO].\n- [81] R. L. Gilliland, P. E. Nugent, and M. M. Phillips, HighRedshift Supernovae in the Hubble Deep Field, ApJ 521 , 30 (1999), arXiv:astro-ph/9903229 [astro-ph].\n- [82] A. G. Riess et al. (Supernova Search Team), The farthest known supernova: support for an accelerating universe and a glimpse of the epoch of deceleration, Astrophys. J. 560 , 49 (2001), arXiv:astro-ph/0104455.\n- [83] A. G. Riess et al. (Supernova Search Team), Type Ia supernova discoveries at z > 1 from the Hubble Space Telescope: Evidence for past deceleration and constraints on dark energy evolution, Astrophys. J. 607 , 665 (2004), arXiv:astro-ph/0402512.\n- [84] A. G. Riess et al. , New Hubble Space Telescope Discoveries of Type Ia Supernovae at z > =1: Narrowing Constraints on the Early Behavior of Dark Energy, Astrophys. J. 659 , 98 (2007), arXiv:astro-ph/0611572.\n- [85] N. Suzuki et al. (Supernova Cosmology Project), The Hubble Space Telescope Cluster Supernova Survey: V. Improving the Dark Energy Constraints Above z > 1 and Building an Early-Type-Hosted Supernova Sample, Astrophys. J. 746 , 85 (2012), arXiv:1105.3470 [astroph.CO].\n- [86] N. Aghanim et al. (Planck), Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641 , A5 (2020), arXiv:1907.12875 [astro-ph.CO].\n- [87] J. Torrado and A. Lewis, Cobaya: code for Bayesian analysis of hierarchical physical models, J. Cosmology Astropart. Phys. 2021 , 057 (2021), arXiv:2005.05290 [astro-ph.IM].\n- [88] J. Torrado and A. Lewis, Cobaya: Bayesian analysis in cosmology, Astrophysics Source Code Library, record ascl:1910.019 (2019).\n- [89] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation  \n- of state calculations by fast computing machines, The Journal of Chemical Physics 21 , 1087 (1953), https://doi.org/10.1063/1.1699114.\n- [90] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57 , 97 (1970), https://academic.oup.com/biomet/articlepdf/57/1/97/23940249/57-1-97.pdf.\n- [91] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys. Rev. D66 , 103511 (2002), arXiv:astroph/0205436 [astro-ph].\n- [92] A. Lewis, Efficient sampling of fast and slow cosmological parameters, Phys. Rev. D87 , 103529 (2013), arXiv:1304.4473 [astro-ph.CO].\n- [93] A. Gelman and D. B. Rubin, Inference from Iterative Simulation Using Multiple Sequences, Statistical Science 7 , 457 (1992).\n- [94] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, Stable clustering, the halo model and non-linear cosmological power spectra, MNRAS 341 , 1311 (2003), arXiv:astro-ph/0207664 [astroph].\n- [95] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Revising the Halofit Model for the Nonlinear Matter Power Spectrum, ApJ 761 , 152 (2012), arXiv:1208.2701 [astro-ph.CO].\n- [96] J. C. Hill et al. , Atacama Cosmology Telescope: Constraints on prerecombination early dark energy, Phys. Rev. D 105 , 123536 (2022), arXiv:2109.04451 [astroph.CO].\n- [97] S. Joudaki, P. G. Ferreira, N. A. Lima, and H. A. Winther, Testing gravity on cosmic scales: A case study of Jordan-Brans-Dicke theory, Phys. Rev. D 105 , 043522 (2022), arXiv:2010.15278 [astro-ph.CO].\n- [98] A. R. Liddle, Information criteria for astrophysical model selection, Mon. Not. Roy. Astron. Soc. 377 , L74 (2007), arXiv:astro-ph/0701113.\n- [99] H. Akaike, A New Look at the Statistical Model Identification, IEEE Transactions on Automatic Control 19 , 716 (1974).\n- [100] G. Schwarz, Estimating the Dimension of a Model, Annals of Statistics 6 , 461 (1978).\n- [101] T. E. Oliphant, A guide to NumPy , Vol. 1 (Trelgol Publishing USA, 2006).\n- [102] S. Van Der Walt, S. C. Colbert, and G. Varoquaux, The numpy array: a structure for efficient numerical computation, Computing in Science & Engineering 13 , 22 (2011).\n- [103] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, ˙ I. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . . Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17 , 261 (2020).\n- [104] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9 , 90 (2007).\n- [105] A. Lewis, GetDist: a Python package for analysing Monte Carlo samples (2019), arXiv:1910.13970 [astro-", 'ph.IM].': "- [106] B. Burlacu, G. Kronberger, and M. Kommenda, Operon c++: An efficient genetic programming framework for symbolic regression, in Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion , GECCO '20 (Association for Computing Machinery, New York, NY, USA, 2020) p. 1562-1570.\n- [107] G. Kronberger, B. Burlacu, M. Kommenda, S. M. Winkler, and M. Affenzeller, Symbolic Regression (Chapman & Hall / CRC Press, 2024).\n- [108] W. G. L. Cava, P. Orzechowski, B. Burlacu, F. O. de Fran¸ca, M. Virgolin, Y. Jin, M. Kommenda, and J. H. Moore, Contemporary symbolic regression methods and their relative performance, CoRR abs/2107.14351 (2021), 2107.14351.\n- [109] B. Burlacu, Gecco'2022 symbolic regression competition: Post-analysis of the operon framework, in Proceedings of the Companion Conference on Genetic and Evolutionary Computation , GECCO '23 Companion (Association for Computing Machinery, New York, NY, USA, 2023) p. 2412-2419.\n- [110] D. J. Bartlett, L. Kammerer, G. Kronberger, H. Desmond, P. G. Ferreira, B. D. Wandelt, B. Burlacu, D. Alonso, and M. Zennaro, A precise symbolic emulator of the linear matter power spectrum, A&A 686 , A209 (2024), arXiv:2311.15865 [astro-ph.CO].\n- [111] D. J. Bartlett, B. D. Wandelt, M. Zennaro, P. G. Ferreira, and H. Desmond, SYREN-HALOFIT: A fast, in-  \n- terpretable, high-precision formula for the ΛCDM nonlinear matter power spectrum, A&A 686 , A150 (2024), arXiv:2402.17492 [astro-ph.CO].\n- [112] E. Russeil, F. Olivetti de Fran¸ca, K. Malanchev, B. Burlacu, E. E. O. Ishida, M. Leroux, C. Michelin, G. Moinard, and E. Gangler, Multi-View Symbolic Regression, arXiv e-prints , arXiv:2402.04298 (2024), arXiv:2402.04298 [cs.LG].\n- [113] S. AbdusSalam, S. Abel, and M. Crispim Romao, Symbolic Regression for Beyond the Standard Model Physics, arXiv e-prints , arXiv:2405.18471 (2024), arXiv:2405.18471 [hep-ph].\n- [114] M. Laumanns, L. Thiele, K. Deb, and E. Zitzler, Combining convergence and diversity in evolutionary multiobjective optimization, Evol. Comput. 10 , 263-282 (2002).\n- [115] M. Kommenda, B. Burlacu, G. Kronberger, and M. Affenzeller, Parameter identification for symbolic regression using nonlinear least squares, Genetic Programming and Evolvable Machines 21 , 471-501 (2020).\n- [116] K. Levenberg, A method for the solution of certain non - linear problems in least squares, Quarterly of Applied Mathematics 2 , 164 (1944).\n- [117] D. W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, Journal of the Society for Industrial and Applied Mathematics 11 , 431 (1963), https://doi.org/10.1137/0111030."}