{'fraction_non_alphanumeric': 0.033001409940077546, 'fraction_numerical': 0.01588385618611209, 'mean_word_length': 4.247745664739885, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New general parametrization of quintessence fields and its observational constraints

We present a new parameterization of quintessence potentials for dark energy based directly upon the dynamical properties of the equations of motion. Such parameterization arises naturally once the equations of motion are written as a dynamical system in terms of properly defined polar variables. We have identified two different classes of parameters, and we dubbed them as dynamical and passive parameters. The dynamical parameters appear explicitly in the equations of motion, but the passive parameters play just a secondary role in their solutions. The new approach is applied to the so-called thawing potentials and it is argued that only three dynamical parameters are sufficient to capture the evolution of the quintessence fields at late times. This work reconfirms the arbitrariness of the quintessence potentials as the recent observational data fail to constrain the dynamical parameters.


I. INTRODUCTION
One of the most famous unsolved mysteries in modern Cosmology is the accelerated expansion of the universe, an observation that has been widely confirmed ever since its discovery in 1998 [1][2][3][4][5][6]. The accelerated expansion is commonly attributed to a mysterious matter component generically dubbed as dark energy (DE). The most accepted DE model is the cosmological constant [7][8][9], which is in fact part of the so-called standard model of Cosmology [6]. From this point of view, a cosmological constant represents a constant vacuum energy which can explain the accelerated expansion very well, but its existence is problematic from the theoretical point of view [7,[10][11][12].
It seems then more natural to consider dynamical models where the DE component could be explained by extra fields in the matter budget or by modifications and/or extensions to our current understanding of the gravitational field [13]. The latter possibility has been just recently weakened by the detection of gravitational waves produced during the collision of binary system of neutron stars, mainly because of the exquisite measurement that confirms that gravitational waves propagate at the speed of light [14][15][16][17][18]. Among the still surviving dynamical models of DE we find in particular those of quintessence scalar fields, which have been present in the literature for almost three decades [19][20][21]. In a quintessence model, a scalar field is minimally coupled to gravity and a potential supply the required negative pressure to drive the accelerated expansion of the universe.
A wide range of quintessence potentials has been proposed in the literature [22][23][24][25][26][27][28] but none of these have a confirmation from the observational point of view. Depending on the evolution of the equation of state parameter of the scalar field quintessence scalar field models are crudely classified into two classes [29][30][31][32] (i) thawing models and (ii) freezing models. For thawing models, the potential becomes shallow at late times and the field gradually slows down. For freezing models, during the early cosmological time, the field is almost frozen due to the presence of Hubble friction and the scalar field starts to slowly roll-down the potential as the field mass becomes lower than the Hubble expansion rate. For a more detailed discussion of the quintessence dynamics we refer to [28,[33][34][35].
In this work we propose a general method to study the evolution of quintessence scalar field models with a general form of the scalar field potential. Using a suitable variable transformation, the equations of motion are written as a set of autonomous equations, which directly suggests a general parametrization of the quintessence potentials without having to know their precise form. Such a dynamical systems analysis of DE models is already popular in the study of cosmology, for examples and references see [35][36][37][38], but so far they have mostly used the original change of variables firstly introduced in Ref. [22]. The particular, polar form of the transformation into a dynamical system which we use in this work was first proposed for dark matter models and the inflationary scenario in [39,40], but see also [41,42] for other related works. As mentioned above, we shall show that there is a general parametrization from which almost all the popular quintessence potentials can be derived. The new parameters are the responsible for the dynamical behavior of the quintessence models, and we use this property to put observational constraints on their values and from this infer the functional form of the potentials that seem to be preferred by the data.
A summary of the paper is as follows. In Sec. II, we setup the mathematical background of the system. This section includes the formation of the autonomous system, its polar transformation, and a description of the general parametrization of the quintessence potentials. In addition, we provide a generic method to infer the potentials from the reverse integration of the given parametrization. In Sec. III we find approximate solutions to the equations of motion in their polar form to follow the evolution of the quintessence variables from the radiation dominated era up to the present time. As a result, we obtain analytical expressions that link initial values of the variables with present quantities of physical interest that can be used reliably in numerical solutions. Section IV is devoted to the numerical analysis of the quintessence solutions and their implementation in an amended version of the Boltzmann code CLASS. We also propose a new parametrization of the DE equation of state that suits well the behavior of the quintessence models, we study this by making a comparison with full numerical simulations of the equations of motion. The comparison with diverse cosmological observations is presented in Sec. IV by means of a full Bayesian analysis, in order to put constraints on the dynamical parameters of the quintessence models. Finally, we present a summary of our results and an outlook for future research in Sec. VI.

II. MATHEMATICAL BACKGROUND
We consider a flat Friedman-Lemaître-Robertson-Walker universe which is dominated by the standard matter fluids plus a quintessence scalar field. We also consider that all the component of the Universe are barotropic in nature, i.e. the pressure p j and the density ρ j are related one to each other by the expression p j = w j ρ j , where w j are the corresponding (constant) equation of state (EOS) parameter of each component. The Einstein field equations, the continuity equation for each matter fluid, and the (wave) Klein-Gordon equation of scalar field can be written, respectively, as where κ 2 = 8πG, H ≡ȧ/a is the Hubble parameter and a the scale factor of the Universe, V (φ) is the scalar potential and a dot means derivative with respect to cosmic time. The scalar field energy density ρ φ and pressure p φ are expressed in terms of the field variables, respectively, as In contrast to the corresponding quantities of the barotropic perfect fluids, the quintessence density and pressure cannot be handled independently and one necessarily requires to find separate solutions for φ andφ from Eq. (1d). In the sections below we will present new variables that help to solve Eqs. (1) more easily.

A. Dynamical system approach
To write Eq.(1d) as a set of autonomous equations, we introduce a new set of dimensionless variables [22,40,43] x As a result, Eq.(1d) is transformed into where now a prime is the derivative with respect to the number of e-foldings, N ≡ ln(a/a i ) and a i is the initial scale factor of the universe. In writing Eqs. (4) we have used the Friedmann constraint (1a) in the form Ω r +Ω m + Ω φ = 1, where the density parameters of the different matter fields are defined in the standard way as Ω j = κ 2 ρ j /(3H 2 ). In addition, the total EoS is given by Equations (4) have been thoroughly used in the literature to study the properties of quintessence potentials, see for instance Refs. [22,35] and references therein. Their main advantage is the possibility to consider a compact phase space for the variables x and y, so that all trajectories and critical points of interest can be studied without the intrinsic difficulties in the standard variables (φ,φ). One main limitation of this approach is that the form of the new potential variables y 1 and y 2 have to be calculated individually for each quintessence potential, and in this respect they do not offer a clear advantage over the direct solution of the KG equation (1d), and the solution of the latter still is the popular approach in the numerical studies of quintessence models.

B. Polar form of the equations of motion
We now introduce the following polar transformations of the variables: x = Ω 1/2 φ sin(θ/2) and y = Ω 1/2 φ cos(θ/2), where Ω φ = κ 2 ρ φ /3H 2 is the density parameter associated to the quintessence field and θ represents an angular degree of freedom. The system of equations (4), after simple manipulations, reduces to θ = −3 sin θ + y 1 , (6a) The EoS of the scalar field in terms of the polar variable is w φ = p φ /ρ φ = (x 2 − y 2 )/(x 2 + y 2 ) = − cos θ, which tells us of the direct relation between the two variables. Likewise, the ratio of kinetic and potential energies is given by tan and (6c) are the same for any kind of potential, and it is only Eq. (6b) that changes for different cases because of the presence of y 2 .

C. General form of the quintessence potentials
To find a solution of Eqs. (6) one needs to close the system of equations, and this can be done whenever the second potential variable y 2 can be written in terms of the variables θ, y 1 and Ω φ . This is equivalent, in the standard approach, to the fixing of the scalar field potential. Then, one possibility is to choose the functional form of the potential V (φ) and to derive from it the form of y 2 following the prescriptions in Eqs. (3). In Table I we give a list of thawing and freezing quintessence potentials that are very familiar in the current literature and their corresponding closed form in terms of y 2 . For those potentials, the dynamical system (6) becomes an autonomous one upon which we can use the known mathematical tools of such systems [44]. It must be noticed that our classification in thawing and freezing is based upon on the behavior of the solutions as described in the corresponding references, but such classifications cannot be read directly from the final form of y 2 , as also a proper choice of initial conditions must be taken into account. More details can be found in Sec. III below.
One can see that for the examples in Table I the functional forms of y 2 can be expressed in terms of the variables y and y 1 , or more precisely, as a polynomial in terms of the ratio y 1 /y. It is then natural to consider that there exists a more general function of y 2 in the form where α i are constant coefficients. As also shown in the examples in Table I, the constant coefficients α i will then drive the dynamics of the quintessence field, although they will not be directly related to other free parameters in the potential, which are denoted with Latin capital letters in the examples of Table I 1 .
For completeness, we show in Table II the inverse process that can be used upon Eq. (7) to obtain from it different quintessence potentials. The simplest possibility is α j = 0, for which Eq. (7) can be written as ∂ 2 φ V 1/2 = 0 and then upon integration we obtain V (φ) = (A + Bφ) 2 , where A, B are integration constants. This is precisely one particular example (with α 2 = 0) of Class Ia in Table II.
In the most general case, Eq. (7) can be written as a differential equation by means of the definitions of the variables y, y 1 and y 2 in Eqs. (3). Hence, where the derivatives are calculated with respect to the dimensionless variable κφ. Using the auxiliary function Thus, the inverse process to find the quintessence potential if the dynamical parameters α j are given consists in the integration of the fundamental equation (9). In general terms, Eq. (9) can be integrated by the method of partial fraction decomposition, for which we require first to find the roots of the polynomial on the right hand side. Once a solution is found for the auxiliary function λ = λ(κφ), the corresponding quintessence potential is obtained from V (φ) = exp −λ(κφ)/ √ 6 . The cases in Table II are those that correspond to the quadratic expansion (α j = 0 for j ≥ 3) in Eq. (9). It can be verified that there is a direct correspondence between the general cases in Table II with the particular examples  shown in Table I, as long as the constants A, B and C are adjusted accordingly. From the numerical point of view, the most general form of the potential is obtained from α j = 0, and all other forms should be a subclass of this. But analytically this is not achievable as the integration scheme is different for different choices of α j , and this is why we find it more natural to classify the quintessence potentials in the four classes shown in Table II.

D. Dynamical and passive parameters
We said before that the α-parameters are the only dynamical ones, and then their allowed values suggest natural classifications of the potentials in general classes. We TABLE I. List of quintessence potentials and their corresponding closed form of y2 in terms of the potential variables y and y1. In general, we see that y2/y is represented by a polynomial form in terms of the variable y1/y. Also, it should be noticed that some of the potential free parameters, indicated here by the capital Latin letters, as in the case of the scale energy A 4 , do not appear in the final form of y2. The only dynamical parameters in the potentials, that end up in the final form of y2, are indicated by the Greek letter λ (although it should not be confused with the variable defined in Eq. (9)).

Ref.
Potential V (φ) Closed form of y2 Thawing potentials [43,45] A 4 (1 + Bφ) 2λ A 4 e λ 1 κφ + e λ 2 κφ 6λ1λ2y + √ 6(λ1 + λ2)y1 + 1 2 y 2 1 /y hereafter dub them dynamical parameters. Apart from this, there will be constants of integration (A, B and C in the examples of Table II), which are then redundant from the dynamical point of view and do not have any influence of the behavior of the field solutions, except in the set up of the initial conditions. We will refer to them as passive parameters.
To briefly illustrate the difference with respect to the dynamical ones, let us consider the well known example in Class Ia which is the quadratic potential V (φ) = (m 2 φ /2)φ 2 (A = 0, B = m φ / √ 2 and α 2 = 0), where m φ is the mass of the scalar field [40]. It can be shown that for this case y 1i = 2 √ 2B/H i = 2m φ /H i , where H i is the initial value of the Hubble parameter (see also Sec. IV below). The initial values of the other two dynamical variables, θ i and Ω φ,i must be fixed by taking additional considerations, like the expected contribution of the quintessence field at the present time. Thus, the value of parameter B plays a role in the set up of the initial conditions of the field variables, even though it does not affect at all their evolution and dynamics. More about the free parameters in the potentials, both dynamical and passive, is discussed in the sections below.

III. APPROXIMATE SOLUTION OF THE EQUATIONS OF MOTION
Here we will show how to obtain a solution of the equations of motion (6) that is of general applicability to any kind of quintessential potential of the (monotonic) thawing type. This will in turn be useful also to obtain appropriate initial conditions for the general numerical solutions that will be used in Sec. IV.
It must be stressed out that the thawing condition for the quintessence models requires that, initially, the EoS is w φ −1 and also that y 1 > 0. This means that the quintessence EoS will deviate from the cosmological constant value at late times. Notice that the thawing condition is here chosen by hand, but our formalism also allows other possibilities (freezing, tracker, skater, etc.), which we leave for future studies.
We start by noting that from observations we expect the present value of the quintessence EoS to be w φ −1, which is equivalent in terms of the polar variables to θ < 1. Moreover, at the epoch when the universe entered into the matter dominated phase from a radiation dominated phase, the scalar field energy density was still very subdominant Ω φ 1. By taking into account the approximations Ω φ 1 and θ 1, we neglect the second term in Eq. (6b) and find separate solutions during the radiation and matter domination eras. We shall then match the separate solutions at the time of the radiationmatter equality, and with this we shall try to make a reasonably good guess of the initial conditions of the universe which can lead to the present accelerated universe.
One final note is that in the radiation and the matter dominated cases the e-foldings N are different. For radiation domination N r = ln(a/a i ) where a i is the initial value of the scale factor, whereas for matter domination N m = ln(a/a eq ), where a eq is the scale factor of the universe at the epoch of radiation-matter equality.

A. Radiation dominated era
As the universe is dominated by radiation the total EoS simply is ω tot = 1/3, and due to the smallness of θ we can use the following approximations: sin θ θ and cos θ 1, so that also w φ −1. Hence, Eqs. (6) reduce to  (7) and (9). Notice that we only considered the expansion in Eq. (7) up to the second order, as the inverse process is not analytical if higher order terms are included.

No
Structure of y2/y Form of the potentials V (φ) The growing solution of Eqs. (10), within the radiation domination era, are given by where a subindex r denote the solution during radiation domination and a subindex i denote the initial value of the corresponding variable. In addition, we also find that y 1 = 5θ, which is just the attractor solution for these variables during radiation domination.

B. Matter dominated era
As the universe is dominated by matter we now consider that ω tot = 0, and after using the same approximations as in Eqs. (10), Eqs. (6) now become After solving Eqs. (12) we obtain the matter dominated solutions Here, a subindex m denote the solution during matter domination and a subindex eq denote the initial value of the corresponding variable at the time of radiationmatter equality. In contrast to the previous radiation dominated case, we are not neglecting the decaying solution in Eq. (13) as it will be required below to handle the transition between the two cosmological eras. We matched the approximate solutions (11) and (13) at the time of radiation-matter equality a eq = Ω r0 /Ω m0 so that we can find a solution at matter domination that carries information about the initial conditions set up in radiation domination. From Eqs. (11) we find the values of the variables at radiation-matter equality: θ eq = θ i (a eq /a i ) 2 , y 1eq = 5θ eq = y 1i (a eq /a i ) 2 and Ω φeq = Ω φi (a eq /a i ) 4 , which we substitute in Eqs. (13) to obtain (14c)

C. Estimation of initial conditions
Equations (14) can be used to estimate the initial conditions on the dynamical variables from the present values of θ and Ω φ . We will assume that the matter domination solutions (14) are valid up to the present time. This is not correct from a formal point of view, but we have verified the appropriateness of Eqs. (14) by a direct comparison with numerical solutions, see Sec. IV below. Hence, by taking a = 1 in Eqs. (14), the initial conditions for the quintessence dynamical variables can be estimated from where we have only taken the leading term in Eq. (14a) for the expression of θ i . The same procedure applied to Eq. (14b) gives y 1i = 5θ i , which is exactly the attractor solution expected during radiation domination, see Eqs. (11).

D. Final considerations
We learn from the first of Eqs. (15), because of the direct relation of the polar variable θ to the scalar field EoS, that the present value of the latter w φ0 is an output value that is directly determined by the initial value θ i . This would be in agreement with the standard field approach to quintessence, in which one has to try different initial conditions, for both φ i andφ i , to explore a given range of values for w φ0 . The main difficulty in the latter is that there is not a straightforward relation between the pair (φ i ,φ i ) and the present value w φ0 , and then the search of initial conditions for a proper sampling of w φ0 must be done differently for each potential V (φ). One big advantage of our approach in this respect is that we can avoid such a hassle and use generic initial conditions for all cases, irrespectively of the particular form of the potential.
On the other hand, the relation between the present and initial values of the scalar field density parameter Ω φ is just the one that is obtained for a cosmological constant; this is hardly surprising as one assumption was that the scalar field EoS was close to −1 for most of the evolution of the Universe.
One final quantity of interest is the ratio y 1 /y at the present time; from the solutions presented above we find that y 10 y 0 Its present value is basically set up by the initial conditions, or equivalently, as suggested by the last equality in Eq. (16), by the present values of the quintessence parameters. Hence, the ratio y 1 /y should remain small for most of the evolution of the Universe if the quintessence field is to be the dark energy, i.e. if 1 + w φ0 ∼ 0. This reinforces our assumption that it is just enough to consider an expansion up to the second order in Eq. (7), and then the dynamics of the quintessence fields will be represented, in general, by the values of the first three coefficients α 0 , α 1 and α 2 .

IV. GENERAL PROPERTIES OF QUINTESSENCE MODELS
In this section, we shall study the evolution of the universe considering the general form of y 2 = α 0 y + α 1 y 1 + α 2 y 2 1 /y, and we shall explain the general procedure to constrain the dynamical parameters α.

A. Class Ia and the cosmological constant
We explain here the correspondence, within our approach, between quintessence and the cosmological constant. Let us start with the simplest possibility which is y 2 = 0, in terms of the expansion of y 2 in Eq. (7) this also means that: α 0 = α 1 = α 2 = 0. As shown in Table II this case corresponds to the parabolic potential V (φ) = (A + Bφ) 2 , where A and B are integration constants. The potential has its minimum at φ c = −A/B, and then the mass of the quintessence field is simply given by m 2 φ ≡ ∂ 2 φ V (φ c ) = 2B 2 . Thus, parameter B gives us information about the mass scale of the quintessence field, whereas parameter A tells us about the displacement of the minimum away from the origin at φ = 0. Moreover, there is now a straightforward interpretation of one of the potential parameters: y 1 = 2 √ 2B/H, see Eq. (3b), and then we find that at all times y 1 = 2m φ /H. The parameter A is left undetermined as it plays no role in the dynamics of the quintessence field.
Let us in addition impose B = 0, which also means that y 1 = 0 for the whole evolution. Equation (4b) is identically satisfied, whereas Eq. (4a) provides the solution tan(θ/2) = tan(θ i /2) (a/a i ) −3 ; we find that θ → 0 as (a/a i ) → ∞, and then also that w φ → −1 at late times. This case corresponds to the case V (φ) = const., that is, to the so-called skater models discussed in [56,57]. Skater models then belong to our Class Ia of quintessence potentials under the condition y 1 = 0. 2 We now revise the case of a constant EoS w φ = w φ0 . Although this can be seen as the simplest generalization of the cosmological constant, from Eqs. (6) we find that a constant EoS in the quintessence case could be obtained if, apart from y 1 = 0, we also impose θ i = − cos −1 (−w φ0 ) and θ = 0. However, the latter condition cannot be sustained if θ = 0, see Eq. (6a), and then we get the same situation for skater models described in the above paragraph in which θ → 0. That is, the quintessence EoS cannot remain constant and must evolve towards the value w φ → −1. These same calculations show that the only consistent conditions for a constant EoS are θ = 0 and θ = 0, which correspond exactly to the cosmological constant case.
The morale from this discussion is then twofold. First, it is not possible to find a quintessence solution that emulates a constant EoS for DE apart from the cosmological constant case. And second, in terms of the parameters in our approach, the cosmological constant is just the null hypothesis (θ = 0 and y 1 = 0), and then any deviation from the null value of the dynamical variables and parameters will be a measurement of the preference of the data catalogs on the quintessence models.

B. The quintessence EoS
One of the primary cosmological parameters in the studies of DE is the EoS of the DE field. In our approach, the EoS is, through the relation w φ = − cos θ, 2 The evolution equation for the variable θ in this case simply is θ = −3 sin θ, which in terms of the EoS can be rewritten as w φ = −3(1 − w 2 φ ). The foregoing equation is exactly the one for skater models [56]. also one of the dynamical variables to describe the evolution of the quintessence field. Here we will discuss the influence of the dynamical parameters α and in doing so we will determine the behavior of the EoS under general quintessence potentials.
We show in Fig. 1 the plots of 1 + w φ as a function of the redshift z for different values of the dynamical parameters α 0 , α 1 and α 2 , respectively. The plots are for the quantity 1 + ω φ instead of ω φ as we are interested about the deviations of the quintessence EoS from the cosmological constant case. But also because at present we can make the approximation 1 + w φ0 θ 2 0 /2, and we see that it is variable θ that provides such deviations. The numerical solutions are grouped according to the Class in Table II they belong to and for the indicated values of the dynamical parameters.
From these plots, one can clearly see that the variation of the EoS is more sensitive to α 2 , and less sensitive to the variation of α 0 . This is just the expected result as α 2 is the partner coefficient of y 2 1 /y 2 in the series expansion (11). It is interesting to note from Fig. 1 that a desired value of EoS of the dark energy can be obtained for a wide range of α parameters. Recent cosmological observations can only constrain the present value of the EOS, but unless there is any constraint on the evolution of the EOS it will not be possible to choose from the solutions for different α parameters. Hence, our expectation is that the statistical analysis using cosmological observations in the next section will not be able to constrain the α parameters. Additionally, we also see that the curves show a monotonic growth if the dynamical parameters are positive, but the curves develop a bump (ie a maximum appears) if the parameters are negative enough. We can only speculate that this latter effect seems to be an indication for the possible appearance of oscillations in the evolution of the EoS, but we will leave this topic for a future study.
As for the monotonic growing behavior at late times that is found for positive values of the dynamical parameters, it can be fit by the following expression, where w 1 and γ are free parameters. Notice that w φ0 is explicitly present in Eq. (17) to ensure that the actual value of the EoS is obtained at a = 1. The second term on the rhs of Eq. (17) corresponds to the expected behavior during matter domination: from the leading solution in Eq. (14a), and together with Eq. (15), we obtain that (1 + w φ ) m (1/2)θ 2 m (1/2)θ 2 0 a 3 . Hence, for those cases in which this approximation is good enough until the present time we expect that w 1 ∼ 1 and γ ∼ 0. Any difference with respect to these values will signal the transition from matter to quintessence domination and of the presence of the dynamical parameters α.
The results from a least-squares fitting of the parametrization (17) to the numerical solutions obtained from CLASS in some selected cases are shown in Ta-  III. The values of the parameters γ and w1 obtained from a least-squares fit of the parametrization (17) to some of the numerical solutions in Fig. 1. Notice that in general γ 1, which means that the leading power in the parametrization (17) is larger than a 3 . The standard errors around the obtained values of the parameters are 1%. dard errors around the obtained values of the parameters are 1%. Not surprisingly, it is consistently found that γ 0, which indicates that the EoS accelerates its growth from −1 as the quintessence field starts to dominate the matter budget. Equation (17) can be compared with other parameterizations of the dark energy EoS, like the famous Chevalier-Polarski-Linder one: w = w 0 +w 1 (1−a), which is clearly inappropriate to describe the evolution of the quintessence models in this work. There exist other parameterizations, see for instance [58][59][60][61][62][63] and references therein, but they usually have a more complicated form than Eq. (17). Although they may serve to test more complicated DE models, they are certainly not the best options to test one DE model as simple as quintessence.
Like in the case of the CPL parametrization, notice that Eq. (17) uses the present value of the EoS as an explicit parameter, but one clear advantage of our approach is that we are parameterizing the underlying dynamical variable θ, and then we are recovering the right behavior of the EoS at early times. Notwithstanding this, we will not pursue a study of the dynamics represented by the parametrized EoS (17) because of the obvious degeneracies with the dynamical parameters: one can see from Table III that different combinations of the α's will result in similar values of the free parameters γ and w 1 . Also, our parameterization (17) is only valid for redshifts z 10, as for larger redshifts we need to take into account the full solutions for radiation domination and the radiation-matter transition, see Sec. III above. All of this makes any reconstruction of the quintessence potential from the EoS parametrization fruitless, and then it is more convenient to work directly with the dynamical variables α extracted from the potentials.

V. OBSERVATIONAL CONSTRAINTS AND RESULTS
Here we discuss our general strategy to put observational constraints on the dynamical parameters that characterize the quintessence field.

A. General setup and datasets
We use an amended version of the Boltzmann code CLASS [64] and the Monte Carlo code Monte Python [65,66]. Amendments to CLASS includes those necessary for MontePhyton to be able to sample the parameters that we describe next. There are 6 parameters that we want to constrain: θ 0 , y 10 , Ω φ0 , α 0 , α 1 , α 2 , but only 5 of them are required as input parameters, namely: Ω φ0 , θ 0 , α 0 , α 1 , α 2 , because the value of y 10 is to be inferred from the full numerical evolution. It must be stressed out that as we sample the values of θ 0 we will also be sampling the present values of the quintessence EoS w φ0 through the relation w φ0 = − cos θ 0 . In practice, the present EoS is then an input value and we will have full control of its sampling, which is another advantage of our method and variables over the standard approach to quintessence fields.
As in many other instances, we still need to finely tune the initial values of the dynamical variables at the beginning of every numerical run. For that, we write y 1i = 5θ i , θ i = P ×Eq. (15a) and Ω φi = Q×Eq. (15b), where the values of P and Q are adjusted with the shooting method already implemented within CLASS for scalar field models. A few iterations of the shooting routine are enough to find the correct values of θ i , y 1i and Ω φi that lead to the desired Ω φ0 and w φ0 with a very high precision; in all instances it has been found that P, Q = O(1), which indicates that Eqs. (15) are good approximations to the required initial conditions. Here we only consider the background dynamics of the quintessence fields and leave the study of their linear perturbations for a future work.
In doing a full sampling of the dynamical parameters α 0 , α 1 , α 2 , we will also be sampling the general form of the potentials shown in Table II. This way we expect to be able to impose constrains on the dynamical parameters but not on the passive ones of the potential V (φ). As explained before, these other parameters are related and can obtained from the dynamical variables θ 0 , y 10 and Ω φ0 , although this would have to be done case by case for each one of the potentials in Table II. For purposes of generality, we will focus on the constraints to the dynamical parameters and consider only two examples of constraints on passive parameters.

B. General results
The general constraints on the parameters of the quintessence models are shown in Fig. 3, where the 1D and 2D posterior distributions are represented in a triangle plot;it is also shown the Mean Likelihood Estimate (MELE) (dashed lines), which is another output from the Montepython code. In what follows we report our results using the median values of the 1D Posterior (not the mean likelihood) plus/minus a confidence interval, which is defined as the range containing 90% of the samples. This particular choice is because some of the posterior parameters are not Gaussian. All of our analyses achieved a convergence ratio (Gelman-Rubin criteria) of R − 1 ≈ 0.005 for the standard cosmological parameters, although for some of the scalar field parameters the best convergence ratio we got did not go below R − 1 ≈ 0.01. Notice that all parameters are well constrained to some region of the parameter space, in both the posteriors and the MELE, except for the dynamical parameters α whose posterior distributions are plainly flat along the full prior range. This means that the data sets considered does not show any preference for a particular quintessence potential.
The present contribution of the quintessence field to the matter budget results in Ω φ0 = 0.719 +0.015 −0.015 , which is in agreement with previous studies [6]. Similarly, we find that the EoS 1 + w φ0 < 0.107 and that 0 < y 1 < 2.24 (95% C.L.), values that are in agreement with the cosmological constant value w φ = −1 and y 1 = 0. These results together show that the quintessence models revolve around the cosmological constant values.
We now go back to the flat posteriors of the dynamical parameters α. It means that a solution of the quintessence field compatible with the observational dataset can always be found for any value of the dynamical parameters, and the reason behind such result is that the initial conditions of the quintessence variables can be finely tuned accordingly to compensate for any α = 0. For instance, for larger values of any of the dynamical parameters we can start the field evolution closer to the cosmological constant case, so that initially w φ → −1 as much as necessary.
Moreover, the flat posteriors in Fig. 3 also imply that there is not clear preference for any of the classes of potentials in Table II. Given this situation, it may be reasonable to just consider the most economic possibility which is Class Ia in Table II: V (φ) = (A + Bφ) 2 . As discussed in Sec. II D, one actually recover the quadratic potential if A = 0 and B = m φ / √ 2, and then we can say that for practical purposes no quintessence potential can fit the data any better than the quadratic potential.
We now turn our attention to the passive parameters in the quintessence potentials. In contrast to dynamical parameters, we shall argue that the passive ones can be subjected to observational constraints. It must be noticed that passive parameters, in their role as integration constants in Table II, can only be determined if we fix either the initial or the final conditions in the solutions of the equations of motion. In the cosmological context we are interested in the final conditions as it is necessary to adjust the parameters to recover the present values of different observables.
Taking as a reference the quadratic potential again, for which the mass of the scalar field m φ is a passive parameter, we show in the top panel of Fig. 4 the posteriors of different cosmological quantities. The fit indicates that Ω m0 = 0.304 +0.018 −0.016 , Ω φ0 = 0.695 +0.016 −0.018 , 1 + w φ0 < 0.129, and m φ < 1.36 × 10 −33 eV (95% C.L.). It can be seen that the preferred value of the EoS is close to −1, and that the scalar field mass m φ has an upper bound. The latter constraint can be easily understood if we recall that the field mass can be calculated from the expression m φ = (1/2)y 10 H 0 , and then any bounds on the scalar field mass are directly obtained from those on the present values y 10 and H 0 .
Other cases, though, are not as clear as the quadratic one. Let us consider the axion potential, which we write as where m φ is the mass of the axion field and f φ is its so-called decay constant. The axion potential belongs to Class IIa with α 0 = 3/(κ 2 f 2 φ ) (together with α 1 = 0 = α 2 ), which indicates that, according to our classification, f φ is a dynamical parameter, whereas the combination m 2 φ f 2 φ forms a passive one. The immediate result is that f φ cannot be constrained from cosmological observations, which is at odds with previous results in the literature (see [74] and references therein). Our interpretation here is that previous studies were not able to sample all possible values of f φ , mostly because of the intrinsic difficulties in solving the quintessence field equations in their normal form, see Eq. (1d). Smaller values of f φ require the field to start closer to the top of the potential, and this is a tough numerical task even in our approach.
As for the passive parameter in the axion potential, it can be shown that it can be determined from [40] The passive parameter of the quintessence potential cannot be written solely in terms of the cosmological observables (H 0 ) and the dynamical variables (w φ0 , y 10 , ρ φ0 ), and then it cannot be clearly constrained because of the presence of the decay constant f φ . Notice however that if f φ 1 (in appropriate units), which corresponds to where ρ φ0 is the present quintessence density, and in this limit there appears an upper bound for the passive parameter basically inherited from the one on ρ φ0 . Likewise, if f φ 1, corresponding to α 0 → 0, we find that m φ f φ (y 10 H 0 /2)f φ , which shows that the passive parameter in this limit will be unbounded from above and that m φ (y 10 H 0 /2), which is exactly the result obtained for the quadratic potential.
The posterior distributions for the axion case are shown in the bottom panel of Fig. 4. The fit indicates that Ω m0 = 0.288 +0.037 −0.035 , Ω φ0 = 0.712 +0.035 −0.037 , 1 + w φ0 < 0.154, and m 2 φ f 2 φ < 3.66 × 10 −10 eV 4 (95% C.L.). Apart from the upper bound for the passive parameter m φ f φ , the preferred values of the other parameters are similar to those of the quadratic potential (top panel in Fig. 4) and also to those of the general case shown in Fig. 3. Hence, the study of particular cases does not provide stronger bounds for the cosmological parameters.

VI. CONCLUSIONS
In this work, we have presented a general study of quintessence dark energy models that allows a general comparison with observational data without the need to specify their functional form. This is possible because the equations of motion of the quintessence field are written as an autonomous system and later transformed to a polar form that automatically satisfies the Friedmann constraint. Moreover, one of the new dynamical variables   in the polar form is directly related to the quintessence EoS, which then means that the latter is no longer a parameter derived from the field equations but rather one that controls the evolution of the quintessence field.
One interesting finding of this work is the general form of the quintessence potentials. To close the polar system of equations one needs the information about a second potential variable that we called y 2 . The functional form of y 2 depends on the particular choice of quintessence potentials, but by observing the results obtained from different potentials we proposed a series form of y 2 that covers a wide range of models. We have correspondingly identified four different classes of quintessence potentials in terms of the series coefficients of y 2 , which is integrated back to get the functional form of the quintessence potentials V (φ) that belong to the four classes.
We have found a general solution of the equation of motion in their polar form by taking into account the fact that the quintessence EOS is very close to −1 and subdominant in both the matter and radiation dominated eras. This solution is particularly interesting as it estimates the information about the initial conditions of the quintessence variables deep inside the radiation era by using the present values of the cosmological parameters. We have incorporated the expressions of the initial conditions in an amended version of the Boltzmann code CLASS, with which we have worked out the numerical solutions of the polar equations of motion. This has allowed us to find a parameterization for the evolution of the EoS that seems to suit better the case of thawing quintessence than others proposed in the literature. The parametrization works well because is based on the analytical solutions found for the polar variables.
However, we did not consider the new parametrization of the EoS for a comparison with observations, but we rather worked directly with the polar equations of motion. According to our study of the quintessence potentials, we distinguished two separate set of parameters in them: the dynamical ones and the passive ones. The dynamical ones appear explicitly in the equations of motion and then have a direct influence on the evolution of the field variables. In contrast, the passive ones are integration constants that can be expressed as combinations of the polar variables and other cosmological variables like the Hubble parameter. The comparison with observations showed that the passive variables can in principle be constrained, but that is not the case of the dynamical parameters in the quintessence potentials, whose posteriors are fully flat. We have verified that this is in agreement with other results already published in the literature. This is one of our main results: that observations cannot establish a preference for a given functional form of the quintessence potential.
The dynamical variables were constrained but their allowed values are close to those of a cosmological constant and in this sense our analysis does not show any preference for quintessence models over a constant dark energy density. As a side result, we have also argued that the results on the dynamical variables can be used to put constraints on the passive parameters of the field potentials. This was done for a couple of particular examples, but our methods can be used for other types of potentials as well.
In all, our results indicate that there will be always a set of dynamical parameters which will satisfy the observational constraints for any given potential. According to our method, this is because our current observations can only put an upper bound on the present value of the DE EoS, 0 ≤ 1+w φ0 < 0.107 (in the general case, see Sec. V B above). The degeneracy in our results could be broken if there were any indication of a non-zero lower value in the EoS (which would, in turn, rule out a cosmological constant), as this will narrow the possible evolutionary paths of the quintessence variables and in consequence the allowed values of the dynamical variables. But given the current state of affairs, we cannot but to conclude that the problem with the arbitrariness of the functional form of the quintessence potential still remains unsolved. To show more about the convenience of Eq. (7) as a general representation of quintessence potentials, we first write Eq. (4c) as and in combination with Eqs. (4b) and (7) we find If we use again the function defined in Sec. II C, y 1 /y = −( √ 6/κ)∂ φ ln(V ) = λ, Eq. (A2) can be written in the form 3 which is known as the tracking parameter [53,[75][76][77][78]. A direct comparison between Eqs. (A2) and (A3) gives which shows the direct relation between our new potential variable y 2 and the tracking parameter. Some previous works have considered that for selected scalar field potentials there is a closed form of the tracking parameter Γ(λ) in terms of λ [53], and for those same potentials our dynamical system (4) becomes an autonomous one because y 2 = y 2 (θ, y 1 , Ω φ ). Our method in this paper suggests that we may as well consider not the complete form but just a series expansion of Γ(λ) to find general solutions of quintessence potentials.
Finally, Eqs. (4a) and (4b) are also rewritten as which resemble the dynamical system of an exponential potential firstly studied in Ref. [22].