Constraints for the thawing and freezing potentials

We study the accelerating present universe in terms of the time evolution of the equation of state $w(z)$ (redshift $z$) due to thawing and freezing scalar potentials in the quintessence model. The values of $dw/da$ and $d^2w/da^2$ at scale factor of $a = 1$ are associated with two parameters of each potential. For five type scalar potentials, the scalar fields $Q$ and $w$ as function of time $t$ and/or $z$ are numerically calculated under the fixed boundary condition of $w(z=0)=-1+\Delta$. The observational constraint $w_{obs}$ (Planck 2015) is imposed to test whether the numerical $w(z)$ is in $w_{obs}$. Some solutions show the thawing features in the freezing potentials. Mutually exclusive allowed regions in $dw/da$ vs. $d^2w/da^2$ diagram are obtained in order to identify the likely scalar potential and even the potential parameters by future observational tests.


I. INTRODUCTION
Although almost two decades have passed since the observation of the accelerating universe, dark energy is still left as one of the most mysterious problems in physics [1]. We investigate the scalar field in quintessence scenario how relevant it is to the dark energy [2][3][4]. In this scenario, if the potential of the scalar field has n independent parameters, the first, second, · · · and n-th time derivatives observation of the equation of state w are enough to specify the scalar potentials and to predict the higher derivatives [5,6]. It has been said that the potentials could be classified to freezing type and thawing type [7]. The equation of state w will approach to −1 in the freezing type and w will depart from −1 in the thawing model. We have adopted three types for freezing potential as V = M 4+α /Q α (inverse power law; α > 0) [8], V = M 4 exp(βM/Q) (exponential; β > 0), and V = M 4+γ /Q γ exp(4πQ 2 ) (mixed; γ > 0) [9], and two types for thawing potential as V = M 4 (cos(Q/f ) + 1)) (cosine) and V = M 4 exp(−Q 2 /σ 2 ) (Gaussian) [10]. They are rather simple potentials and have two parameters to specify the form of potential, respectively.
By adopting the values of the derivatives dw/da and d 2 w/da 2 at the scale factor of a = 1, the characteristic form parameters of each potential could be fixed [5,6] and then the theoretical time development of the scalar field Q in the potential could be numerically calculated from the present to the past (and then the reverse).
Through the investigation of the adopted potentials, we have found that there is a forbidden region in the first and second derivatives of w for each potential. Although there is a lot of uncertainty at the moment, it is estimated the allowed regions for potentials in the derivative space of w such as dw/da − d 2 w/da 2 under the comparison with the numerical results and the observations. The main essence of the allowed region for each potential is presented in Fig. 1.
We have taken the exact approach [11,12] to investigate five different exact potentials, not taken the approximate parametrization approach for w [13,14]. We have also mainly concerned the terminal boundary condition at z ≃ 0 [11][12][13][14][15][16][17][18][19][21][22][23] and not z ≫ 1 [4,11,15,17]. Adding to the above-mentioned points, there are two more new aspects of this paper, 1) We have imposed a new constraint on w(z), which was discovered in Ref. [20] (Planck 2015). The main point of the new constraint is w(z) ≤ −0.91 at z = 0.65 ∼ 2.0, which has not been previously taken into consideration in the literature to our best knowledge. We consider that it is a handhold constraint to investigate the quintessence models.
2) We noticed that many investigated exact potentials, perhaps for simplicity, have two parameters and to determine these parameters from observations, we need dw/da and d 2 w/da 2 , adding w 0 . By assuming two parameters (dw/da, and d 2 w/da 2 ), the potential parameters are fixed and we can calculate the evolution of scalar field Q and variation of w(z). It becomes possible to compare the calculated variation of w(z) with the observed constraint on w(z). From this procedure, we can derive the allowed parameter space in dw/da − d 2 w/da 2 . The point is that the allowed region for each potential is different from each other and a coordinate point is a one-to-one correspondence to a set of potential parameters. In future, we hope that the more detailed direct observations of dw/da − d 2 w/da 2 space become possible and they could distinguish the quintessence different potentials and the potential parameters. So it must be the most fundamental step to understand the detailed dark energy in cosmology.
In this paper, the allowed dw/da−d 2 w/da 2 space was not obtained by direct observations of dw/da nor d 2 w/da 2 but was obtained from the comparison of the observational constraint of w(z ≥ 0) and the numerical simulations via scalar potentials, because the potential parameters were associated with dw/da and d 2 w/da 2 .
In §2 the derivatives of w with a and the basic parameter ∆ are presented. The diagram and meaning of dw/da−d 2 w/da 2 are introduced and the numerical procedure with boundary conditions are explained. The main results are presented in Fig. 1 and the followings are mainly the explanation for them. In §3 the forbidden region for each potential is described and the allowed region is derived through the comparison with the observations [20]. The detailed changes of w(z) with the redshift z for potentials and parameters in the dw/da − d 2 w/da 2 space are calculated. In §4 the parameter ∆ has been changed from 0.1 and the results of the allowed region for each potential are presented. We investigate the effect of another observational constraint [20] for the allowed regions for the potentials in dw/da − d 2 w/da 2 space in §5. In §6 the allowed regions for the investigated potentials are summarized and the possibility of the thawing and freezing type potentials are discussed. Although it is calculated exactly in the text, it is approximated by Taylor expansion for the above potentials in Appendix A.

A. Scalar field
For the dark energy, we consider the scalar field Q, where the action for this field in the gravitational field is described by where S M is the action of the matter field and G is the gravitational constant, usually putting G = 1 [6]. In the expanding universe, the equation for the scalar field Q(t) becomes where H is the Hubble parameter, overdot is the derivative with time, and prime of V ′ is the derivative with Q. This is the equation to calculate the development of the field Q, if the potential is fixed with the estimated parameters.
The energy density and pressure for the scalar field are written by and respectively. Then the parameter w Q for the equation of state is described by It is assumed that the current value of w Q is slightly different from a negative unity by By using Eq. (5),Q 2 is written asQ which becomes,Q where it is used the density parameter Ω Q = ρ Q /ρ c =0.68, being ρ c = 3H 2 /(8π) the critical density of the universe. Combining Eqs. (7) and (8), V is given by From Eqs. (8) and (9),Q is expresseḋ Since ρ c is given by the observation through the Hubble parameter H,Q is determined by Ω Q and ∆, which also determine the value of V . If we adopt the form and parameters of each potential, the value of V could be used to estimate the value of Q. Actually, the evolution of H in Eq. (2) depends on the background densities which include radiation density. The effect of radiation density can be ignored in the near past (z ≤ 10 3 ) and so is not considered in this work. The numerical calculations must be done simultaneously with Eq. (2) and Friedmann Eq. forä as where ρ m is the matter density, including cold dark matter and baryonic matter. We take ρ m (z) = ρ m0 (1 + z) 3 , Ω m0 = ρ m0 /ρ c = 0.32 and Ω Q + Ω m = 1, taking ρ m0 = ρ m (z = 0) and assuming flat universe. It must be noted that if time t is normalized by t n = ( 4πGρ c /3) −1 as τ = t/t n , it becomes (da/dτ )/a = √ 2 at a = 1 for H 2 = ((da/dt)/a) 2 = 8πGρ c /3. :Inverse power law type (brown-curve like region), Exponential type (orange), Mixed type (yellow), cosine type (green), and Gaussian type (blue). In Fig. 1, d 2 w/da 2 are plotted against dw/da for each scalar potential. A series of colored regions are the boundary and the allowed regions of d 2 w/da 2 versus dw/da which were obtained a posteriori, so as to make the numerical calculation w(z) satisfy the observational constraint w obs (z) < −0.91 at z = 0.65−2.0 [20], shown in Fig. 2 (b) by blue step lines. In the previous papers [5,6], dw/da and d 2 w/da 2 were associated with two parameters of each scalar potential. When the fixed values of w(z = 0),Q(z = 0), and V (Q)(z = 0) are taken as the boundary conditions, each potential shape was found to be restricted in a some extent so as to reproduce the w obs (z) band data [20].
C. Implications Figure 1 tells us that the allowed regions are almost mutually exclusive among the several potentials and that a coordinate point of (dw/da, d 2 w/da 2 ) is nearly one-to-one correspondence to a set of the potential parameters. According to the diagram of Fig. 1, the direct observation of (dw/da, d 2 w/da 2 ) enables us to identify not only the type of the scalar potential but also the potential parameters. The allowed regions in the d 2 w/da 2 versus dw/da diagram may be convenient and helpful to explore a likely scalar potential in quintessence models.
In the region dw/da > 0 of Fig. 1, w is increasing at z = 0 which shows a what is called "thawing feature". On the other hand, in dw/da < 0 , w is decreasing at z = 0 which shows a "freezing feature". In d 2 w/da 2 > 0, dw/da is increasing and w(z) indicates a concave curve at z = 0 which shows a "thawing feature" and a "freezing feature" [18]. In d 2 w/da 2 < 0, dw/da is decreasing and w(z) indicates a convex curve at z = 0 which shows a "freezing feature".
So "freezing feature" appears in the lower-left and upper-left regions on dw/da − d 2 w/da 2 space and "thawing feature" appears in the upper-right region on dw/da − d 2 w/da 2 space.
In the other region, w(z) shows a complex feature such as, for example in the right-lower region (dw/da > 0 and d 2 w/da 2 < 0) where the initial conditions must be unnatural (to be discussed later).

D. Numerical procedure
The procedure how to test the potential parameters and obtain the allowed regions in Fig. 1 is as follows. First, a set of values of (dw/da, d 2 w/da 2 ) was chosen with a full plain 6 scan, for example, of −0.4 < dw/da < 1.2 and −1.5 < d 2 w/da 2 < 1.0 by a step of 0.01.
For each scalar potential, the two parameters were calculated as in the previous paper [5] and touched in this paper. Second, a set of the terminal boundary conditions (given below) was imposed to solve the simultaneous differential equations of Eqs. (2) and (11) [20]. The values of dw/da and d 2 w/da 2 were associated with two parameters of each scalar potential as implicit parameters. One should note that the allowed regions are almost mutually exclusive and that a coordinate point is one-to-one correspondence to a set of potential parameters.
reproduces within the observational constraint of w obs (z) [20]. If yes, then the values of (dw/da, d 2 w/da 2 ) are the allowed ones. Thus, the diagram was obtained.

E. Boundary conditions
The first set of the terminal boundary conditions was imposed to solve the simultaneous differential equations of Eqs. (2) and (11) for the investigated potentials, where ∆ = 0.1 and time normalization t n = (3/(4πρ c )) 1/2 = 1 is taken. The initial value of Q(z = 0) is usually derived from V (Q)(z = 0), where the form parameters of each potential are already determined by adopted dw/da and d 2 w/da 2 .
Under the boundary conditions, the potential parameters were restricted to reproduce w(z) within the observational constraint. A non-monotonic solution of Q with respect to z was found to be a key to cause the meandering behavior in w obs . In §4, the effects of variation of ∆ in w = −1 + ∆ were investigated.

III. FIRST AND SECOND DERIVATIVES OF w Q
To fix the characteristic parameters of each potential, we use dw Q /da and d 2 w Q /da 2 [2,5,6] as follows. The first derivative of w Q is given by where M pl is the Planck mass.
From the form of dw Q /da, the second derivative of w Q is given by After the calculation in the paper [6], d 2 w Q /da 2 becomes From this equation, we estimate d 2 w Q /da 2 for each potential [6].
Thus dw Q /da and d 2 w Q /da 2 are associated with the parameters of potential form.
A. Inverse power law potential α is given by Using Eqs. (12), (15), and (16), α is described by d 2 w Q /da 2 and dw Q /da as Then the constraint α > 0 separates the region by the parabolic curve in Fig. 1 as This parabolic criterion for ∆ = 0.1 is shown by red curve in Fig. 1 and the downside region indicates α > 0.
From the previous study [5] This potential has two parameters of α and M. α is estimated by dw/da and d 2 w/da 2 through Eq. (17). V ′ /V is estimated by dw/da through Eq. (12). Using these α and V ′ /V , By takingQ estimated by Eq.  Table I. The observational constraints by blue step lines and numerical calculations of w as a functionof z for the cases a) ∼ e) are shown in Fig. 2 (b). Q as a functionof z is presented in Fig. 2 (c). In Fig. 2 (d), V as a functionof Q is shown where z is an implicit parameter.
In  Fig. 4 (b) by red curves the w evolution of almost 100 cases evenly adopted from the brown region in Fig. 2 (a).  Table I  Confront color sorting with Table I In this left-hand part of the brown region, the signature of d 2 w/dz 2 is negative and they are the curves of convex upward which are represented by blue, and sky-blue color curves. Especially the case a) represented by blue curve shows that w has decreased from upper value to −0.9 at present (z = 0) which meansQ 2 has decreased in the potential V . The freezing potential has the feature that w will decrease and approach to −1.  of which slopes are negative at z = 0 (dw/dz = −(1 + z) −2 dw/da < 0 ). In these cases w has increased from −1 to −0.9 which means thatQ 2 has increased in the potential. The thawing type potential has this feature. The second derivative signature of w with z is positive and they are the curves of downward convex (concave).
In Fig. 3 (a), it can be seen that all cases except b) in Table I have climbed up the potential at first and then passed down. All cases except b) have begun from a negative value of dQ/dt and then increased in Fig. 3 (b). When dQ/dt passes through 0, w takes the value of −1 which is noticed in Fig. 4 (a). The blue a) and sky-blue b) cases have decreased the value of dQ/dt after z ∼ 0.5 which can be seen in Fig. 3 (b). It may be due to the expansion drag effect of the second term in Eq. (2). (c) Relation between dQ/dt and V (Q) for each case. In Fig. 3 (c), it can be seen that dQ/dt has increased for almost all cases, however it has decreased (after z ∼ 0.5) for the blue a) and sky-blue b) cases. (d) Relation between d 2 Q/dt 2 and V (Q) for each case. It can be noticed in Fig. 3 (d) that d 2 Q/dt 2 is positive in almost all situation, however it becomes negative (after z ∼ 0.5) for the blue a) and sky-blue b) cases.
In Fig. 4 (a), under the observational constraint, the w evolution of cases a) ∼ e) are shown, which are the same as Fig. 2 (b), except the part of −1.05 ≤ w ≤ −0.5 in the vertical axis is enlarged. The a) and b) cases in the left-hand part of the brown region (dw/da < 0) in Fig. 2 (a) are represented in Fig. 4 a) by blue and sky-blue curves of which slopes are positive at z=0 (dw/dz = −(1 + z) −2 dw/da > 0). In Fig. 4 (b), almost 100 cases are presented by red curves which are evenly adopted from the brown region in Fig. 2 (a).
Evolutions of d 2 a/dt 2 for cases a) ∼ e) are presented in Fig. 4 (c) where the signature has changed to positive around z = z a ∼ 0.6, meaning to become accelerating universe. The fact that z a is so close to zero is the coincidence problem. In Fig. 4 (d), the evolution of da/dt are shown. Expansion velocity was decreasing at first and then changed to increasing at around z ∼ 0.6. These changes of Figs. 4 (c) and (d) are almost the same for the following cases in Tables II ∼ V, so they are not shown anymore.
Even for the freezing type potential, there are some cases which have the thawing features that w has increased from lower value (≃ −1) to −0.9. These characteristics are common for the following exponential and mixed type potentials.

B. exponential potential
The constraint β > 0 is given from the relation Taking Q > 0, the criterion of β > 0 and α > 0 is the same as the positive signature of the square brackets [ ] and the same region in dw Q /da − d 2 w Q /da 2 plane of Fig. 1.
The cases a) ∼ e) with values of dw/da, d 2 w/da 2 , βM and Q 0 are presented in Table II.
In Fig.5 (a), it is the same as Fig. 1 for the allowed orange colored region for the exponential potential which is enlarged. The cases a) ∼ e) in Table II  for each case is shown in Fig. 5 (c) Three cases ( a), c) and e)) have climbed the potential at first (z ∼ 2) and then gone down. When dQ/dt passes through 0, w takes the value of −1 which is noticed in Fig. 6 (a). In Fig. (d), the evolution of the velocity dQ/dt for each case is presented. Three cases ( a), c) and e)) have started from dQ/dt < 0 at first (climbed the potential) and then increased to become positive.
The observational constraints by blue step lines and numerical calculations of w as a function of z for the cases a) ∼ e are shown in Fig. 6 (a). Many cases in the orange region in Fig. 5 (a) are shown by red curves in Fig. 6 (b) For this potential, the allowed parameter region in dw Q /da−d 2 w Q /da 2 plane is presented by orange color in Figs. 1 and 5 (a). The example curves in Table II    The constraint γ > 0 is given by This relation is reduced to which is almost to the same of the α > 0 constraint (Eq. (19)), however the last constant . This criterion is shown in Fig. 1 by dark-blue parabolic curve.
For this potential, the allowed parameter region in dw Q /da − d 2 w Q /da 2 plane is shown in yellow in Figs. 1 and 7 (a). The typical cases a) ∼ e) with values of dw/da, d 2 w/da 2 , γ, and Q 0 are presented in Table III   All cases except c) have started from dQ/dt < 0 at z = 2 and then become positive. When dQ/dt passes through 0, w takes the value of −1 which is noticed in Fig. 8 (a).
The minimum point of this potential is realized at Q min = (γ/(8π)) 1/2 . The values of  For the cosine type potential, there is the constraint that the following Y must be negative.  Fig. 1 for the allowed green colored region for the cosine potential which must be noticed in different scales. The same with Fig. 1 for sky-blue, red and dark-blue curves. The cases a) ∼ e) in Table IV are  and then slowed down. (d) Evolution of the velocity dQ/dt for each case is presented. Three cases (c) d) and e)) have started from dQ/dt < 0 at and then become positive. When dQ/dt passes through 0, w takes the value of −1 which is noticed in Fig. 10 (a). Especially case e) of red curve seems to be stagnating around some part (top?) of the potential.  It must be noted that, as the signature d 2 w/da 2 is positive, the allowed curves are the downward convex (concave).
The above constraint reduced to the following equation which means that there is an allowed region, given by The critical boundary of the parabolic equation is shown by the sky-blue curve in Fig.   1. For this potential, the allowed parameter region in dw Q /da − d 2 w Q /da 2 plane is shown by green in Figs. 1 and 9 (a). The typical cases a) ∼ e) with values of dw/da, d 2 w/da 2 , f, and Q 0 are presented in Table IV. The example curves in Table IV  The scales must be noticed to be different in the vertical axes of Figs. 1, 2 (a) and 9 (a). Then the allowed region is large for this potential compared with the previous freezing type potential. If we estimate the probability of the potential from the regional surface of dw/da − d 2 w/da 2 space, the probability of this potential is large.

E. Gaussian type potential
For the Gaussian type potential, it is interesting that the constraint σ 2 > 0 is the opposite region of α > 0 in dw Q /da − d 2 w Q /da 2 plane of Fig. 1 [5]. For this potential, the allowed parameter region in dw Q /da − d 2 w Q /da 2 plane is presented by blue in Figs. 1 and 11 (a).
The typical cases a) ∼ e) with values of dw/da, d 2 w/da 2 , σ, and Q 0 are presented in Table   V. The example curves in Table V  One should notice that the scale height of d 2 w/da 2 in Fig. 11 (a) is much higher than the similar figures. If one estimate the probability of the potential from the allowed region surface in dw/da − d 2 w/da 2 space, the allowed surface is the vastest compared with other potential regions. The height of the allowed region becomes much higher than 4000 in numerical simulations. The large value of d 2 w/da 2 means the rapid increase of w around z ≃ 0 (see Fig. 12 (a) especially red case e)).  (9)). If the time normalization Evolution of the potential height V (Q) for each case is shown. Two cases (c) and d) have climbed at first (z ∼ 2) and then slowed down. (d) Evolution of the velocity dQ/dt for each case is presented. Two cases c) and d) have started from dQ/dt < 0 at z = 2 and then become positive. When dQ/dt passes through 0, w takes the value of −1 which is noticed in Fig. 12 (a). Case e) seems to be stagnating around some part (top?) of the potential.   Figs. 11 and 12 (a). The cases of a) to e) are adopted along the blue region in Fig. 11 (a).
It should be noted that the parameter σ and the value of Q 0 at present are different for each case. Although it is difficult to deny the possibility of the large value of d 2 w/da 2 , it seems to be low for its rapid increase of w around z ≃ 0 (see Fig. 12 (a) especially red case e)). At least such a rapid change of w will be checked in near future. For the inverse power potential, the results are presented in Fig. 13 (a). The enlarged allowed region is presented in Fig. 13 (b). The two cases are presented in Figs. 14 (a) and (b). The evolution w with time z is shown in Fig. 14 which is the same with Fig. 4 except for the value of ∆ = 0.0001.
To analyze the pink region, we have presented two cases in Fig. 14 (a) and (b). The red curve is within an allowed region (pink) in Fig. 13 and green curve is out of the allowed region (right-hand side). The value of d 2 w/da 2 is about −10 for both cases.
In Fig. 14 (a), the w value of green one at z = 2 is greater than −0.91 which is not allowed by the observation (blue line). The detailed features around z ≃ 0.05 is shown in Fig. 14 (b). In Fig. 14 (b), it is almost the same with Fig 14 (a), except for the horizontal range where the range 0 ≤ z ≤ 0.1 is enlarged. Both curves show the similar features. Toward the past, the field Q decreases (not shown here) which means that the field Q climbs up the potential V = M 4+α /Q α and then returns back down the potential. At the turnover the value of dQ/dt becomes null and w becomes −1. Then the field Q comes down the potential which shows the bouncing feature of w with time (z).
If the second derivative of w with z is negative, it is the curve of convex upward which is observed far from the bouncing region. Around the bouncing region, it is the curve of convex downward (concave) where the second derivative of w with z is positive, even though the second derivative w with a is negative at . The left-hand side of the allowed pink region is forbidden by the constraint adw/da + 6∆(1 − ∆/2) ≥ 0 which could be understood in Eq. (12).
For the exponential potential, the results are presented in Fig. 15. For the mix type potential, the results are presented in Fig. 16. These freezing type potentials have almost to z ≃ 0, then the extreme convex feature appeared in the redshift evolution of w(z). For the Gaussian type potential, the results are presented in Fig. 18. If one estimate the probability of the potential from the surface of dw/da − d 2 w/da 2 space, the Gaussian potential has the largest possibility to represent the quintessence field, because the region bounded by the red and pink curves is very vast compared with other potentials. The second large possibility is the cosine type potential bounded by red and pink curves in Fig. 17.
In the limit ∆ → 0 the allowed region for every potential includes the part of dw/da ≃ 0 and d 2 w/da 2 ≃ 0. In this limit, it is hard to discern the potential type in the quintessence scenario as well as the cosmological constant model (w = −1).  Fig. 20 by blue lines (z ≤ 2) is the one in the paper (Fig. 7: Planck+BSH) [20]. To some extent, both constraints are similar around z ≃ 0.6, however the adopted constraint of green curve in this section has the long constraint until z = 5. The allowed regions for the investigated potentials shrink as presented in Fig. 19 where the allowed regions are much smaller than those in Fig. 1. The example curves for the case of the inverse power potential are shown in Fig. 20 by red curves. The constraint for the large redshift (z ≤ 5) has the effect on the decrease of the allowed region for each potential in dw/da − d 2 w/da 2 space.
It is noticed that the constraint of z ≃ 5 has the effect to shrink the allowed region, as the brown colored region in Fig. 1 has decreased a little bit in Fig. 19. The allowed regions for the adopted potentials under another observational constraint which is the one in the paper (Fig. 5) [20]. Compare with Fig. 1 where the observational constraint is the one in the paper (Fig. 7) [20]. The parabolic curves are the same with Fig. 1. Each shrunk colored region corresponds to the same potential as Fig.   1, respectively. The brown curve-like region is for the inverse power potential. The orange region is for the exponential potential. The yellow colored region shows the allowed region for the mixed type potential. The blue region is for the Gaussian potential and the green region is the allowed region for the cosine potential. To some extent, both constraints are similar, however the adopted in this section has the long constraint until z = 5 and the allowed regions have shrunk. In the simulation, we take ∆ = 0.1 and Ω Q = 0.68.
It is apparent that other colored region in Fig. 1 has decreased in Fig. 19. On the other hand, the allowed region exists for any potential, even though the region has shrunk.

VI. RESULTS AND DISCUSSION
We study the accelerating universe in terms of the time evolution of the equation of state w(z) for the thawing and freezing potentials in the quintessence scenario.
The investigated exact potentials have two parameters for their forms. To determine the parameters from the observation, we need dw/da and d 2 w/da 2 and we could not find out the papers to refer this point except [5,6]. If we know the exact potential with two parameters, it becomes possible to simulate the evolution of the scalar field with the two cases of the inverse power potential are shown by red color curves. It is noticed that the constraint of z ≤ 5 has the effect to shrink the allowed region, as the brown colored region in Fig. 1 has decreased a little bit in Fig. 19.
boundary conditions Q and dQ/da (it is equivalent to have w 0 and Ω Q (Eqs. (9) and (10)), being the second differential equation, which determines the evolution of Q. In the end, it is necessary to observe four values w 0 , Ω Q , dw/da 0 , and d 2 w/da 2 , which determine the potential parameters and evolution of the scalar field.
The main procedure of this paper is the following: By adopting the values of the derivatives dw/da and d 2 w/da 2 , the characteristic parameters of each potential can be fixed and then the theoretical time evolution of the scalar field Q in the potential is numerically calculated from the present to the past which is time reversible. The time evolutions of w(z) which pass through the observational constraint are selected and those values of dw/da and d 2 w/da 2 are presented as an allowed region for each scalar potential which are summarized in Fig. 1.
By comparing the observations of w(z) variation with the simulations, we can discriminate the quintessence potentials through Fig. 1. For the moment as there is a lot of uncertainty in the observations, every investigated potential has the possibility to explain them, however the form parameters of each potential have been limited which are related to dw/da and d 2 w/da 2 as stated in the paper and especially Appendix B. The parameter space of dw/da − d 2 w/da 2 are limited as shown in Fig. 1. If future observations become improved and more precise in dw/da − d 2 w/da 2 space, it will be possible to discriminate the potential types and potential parameters.
There are some cases which show the thawing feature in the freezing type potentials. It has been assumed that the scalar field flows down the freezing type potential. However, under the observational constraint of z(≤ 2), there is a possibility that scalar field climbs up the potential and stop there, where w becomes −1, and then flows down where w increases from −1. If we consider that such initial conditions are not natural, there will be no allowed dw/dz − d 2 w/da 2 space for the freezing potentials. It is due to the observational constraint that w must be lower −0.91 around 0.65 ≤ z ≤ 2 and w is −0.9 at present (z = 0). If we accept the rather unnatural conditions that Q climbs up the potential, situated there and roles down for the freezing potentials, there is a possibility that freezing type potentials could survive.
Even for thawing potential, there are rather unnatural solutions that Q climbs up the potential, situated there and roles down. These solutions are investigated by Swaney and Scherrer [19]. We accept those solutions from the phenomenological point that the observations do not exclude such solutions for the moment. It will be a future problem why such initial conditions are realized if such solutions are found in observations.
If one estimate the probability of the potential from the surface region of dw/da−d 2 w/da 2 space, the Gaussian potential is the first, because the region bounded by the red and pink curves in Fig. 18 (a) is the largest compared with other potentials. The second possibility is the cosine type potential bounded by red and pink curves in Fig. 17. Under the observations [20], there are two ample space regions in dw Q /da − d 2 w Q /da 2 space for thawing type potentials (Gaussian type and cosine type). Only the relatively small regions are allowed for freezing type potentials.
Although there is a lot of uncertainty at the moment, it seems to be preferable for the thawing model against the freezing model under the comparison with the numerical results and the observations. It should be stressed that the detailed observation of dw/da, and d 2 w/da 2 will determine the potential type of the quintessence scenario as shown in Figs. 1, 9 (a), 11 (a), and 19, respectively.
About the other observational constraint (z ≤ 5) in §5, we hope the constraint becomes much severe and extends to higher z, however the observational constraint must be very hard where the ratio of dark energy to matter is smaller than 10 −2 ∼ 1% at z = 5 as There is a lot of many works related dark energy under a quintessence scenario and this work concentrates on this model. Although one parameter (∆) formula is investigated in the reference [16] and some interesting features are predicted, our work introduces two more parameters and various other possibilities are predicted under the observational constraint [20]. There are many other trials to investigate w(z) of dark energy and they have the original perspective in the quintessence scenario [10][11][12][13][14][15][16][17][18][19][20][21][22][23]. We believe that it is urgent to determine ∆ > 0 for quintessence and other models. Moreover, if ∆ ≤ 0 (w ≤ −1), we have to consider utterly different scenarios, [1,[24][25][26][27][28].
We hope that it will be improved in detail the observation of w(z) variation as presented by Planck (2016) [20] in the following projects; BOSS: Baryon Oscillation Spectroscopic Survey, DES: Dark Energy Survey, WFIRST: Wide-Field Infrared Survey Telescope, Euclid (ESA satellite), and CMB (COrE) (ESA satellite) projects [28].
(d 2 w/da 2 ). As we have expanded w(a) in Taylor expansion [20] given by where a, w 0 , w a , w a2 and w a3 are the scale factor (a = 1 at current), the current value of w(a), the first, second and third derivatives of w(a) as w a = −dw/da w a2 = d 2 w/da 2 and w a3 = −d 3 w/da 3 , respectively. We take it as the following relation of z = (1 − a)/a, being 1 + z = a 0 /a and a 0 = 1.
w(z) = w 0 − w 0 (z/(1 + z)) + 1 2 w a2 (z/(1 + z)) 2 − 1 3! w a3 (z/(1 + z)) 3 , where we neglect the higher expansion term. It is different from the treatment in the text where the exact evolution of w under each potential is calculated. with the first 4 terms (including until the third derivative). The allowed regions for the inverse power law, exponential, cosine, and Gaussian potentials are colored by brown, orange, blue and green, respectively. There is no region for the mixed type potential. In the simulation, we take ∆ = 0.1 and Ω Q = 0.68.
The results are presented in Fig. 21 which is analogous in Fig. 1. In Fig. 21, the allowed regions for inverse power, exponential, cosine, and Gaussian potentials are displayed by brown, orange, blue and green color, respectively. There is no allowed region for the mixed type potential in this Taylor expansion formalism. In the text, the evolution of Q is exactly calculated. Although Taylor expansion is the approximation, however there are some similarities between Fig.21 and Fig. 1, such as the allowed region for the inverse power law potential is almost on the one-dimensional curve.

Appendix B : Derivation of parameters
In this appendix, it is outlined the procedure for the derivation of parameters for each potential. The equation number Eq. (x) in the reference (Muromachi et al. [5]) is referred as Eq. (M x) in this appendix.
Taking this and d 2 w/da 2 , Y = V ′′ /V is derived. Using the relation Y = X 2 + γ/Q 2 + 8π, γ/Q 2 is obtained as in Eq. (M 53). Dividing Eq. (M 50) by Q, X/Q is estimated. As we know X, Q is derived. Then it is easy to get γ, through X or Y .

B. 3 Cosine Potential
By adopting dw/da in Eq.