1 Introduction

Scalar field models have played a vital role in cosmological theoretical studies in nearly half a century. Those assumed scalar fields appeared in different cosmological research aspects to settle different cosmological problems [1], such as to drive inflation, to explain a time variable cosmological \('\)constant\('\) and so on. After the discovery of the accelerating expansion of universe, scalar fields have played another important essential role as a candidate of dark energy. There are so many phenomenological dark energy models of scalar fields, such as quintessence, phantom, quintom and the scalar fields with non-canonical kinetic energy term (for a review, see [2, 3]).

Phase-plane analysis is a very useful and common method (see Ref. [4] and recent papers, e.g. [58]) to study the dynamical evolution of those scalar fields models and their cosmological implications. However, most of those works only focus on the quintessence models (including phantom quintessence and quintom) with unique exponential potential and tachyon models (including phantom tachyon) with inverse square potential, and correspondingly, the dynamical systems are two-dimensional autonomous systems (see the references cited in [11, 27]). Using a method which considers the potential related variable \(\Gamma \) as a function of another potential related variable \(\lambda \) (see Eq. (10) for the definition of \(\Gamma \) and \(\lambda \)) [11, 27, 46], we are able to analyze the phase plane of the dynamical systems of the quintessence and tachyon models with many different potentials. When the potentials are beyond the special type such as exponential or inverse square potentials, the dynamical systems consequently become a three-dimensional autonomous systems. This method is quite effective and powerful, it therefore has been generalized to several other cosmological contexts [1220, 45]. However, there is very few work focusing on the dynamical behavior of the scalar field with a general modified kinetic term, such as K-essence (\(L=V(\phi )F(X)\)) and general non-canonical scalar field (\(L=F(X)-V(\phi )\)). Recently, Josue De-Santiago et al. analyzed the dynamical system of general non-canonical scalar field with the lagrangian \(L=F(X)-V(\phi )\) and studied the phase plane after a suitable choice of variables [40]. They obtained the three-dimensional autonomous system of this non-canonical scalar field after specifying the kinetic term as \(F(X)=AX^{\eta }\) and choosing the potential as \(V(\phi )=V_0(\phi -\phi _0)^{1/(1-\Gamma )}\) (i.e., a special case that \(\Gamma =\mathrm{const}\)) and studied the critical points as well as their stability.

Motivated by Ref. [40], we try to extend our work [11, 27] in this paper to give the three-dimensional autonomous dynamical systems for most of the popular scalar field dark energy models including (phantom) quintessence, (phantom) tachyon, K-essence, and general non-canonical scalar field models in current work. We will show that the three-dimensional autonomous systems of general non-canonical scalar field and K-essence will reduce to the quintessence and tachyon scalar field, respectively. Unlike most of the previous works, here we express the three-dimensional autonomous systems from variables \((x, y, \lambda )\) to the observable related variables \((w_{\phi }, \Omega _{\phi }, \lambda )\). It will be very convenient to investigate the dynamical properties of the autonomous system based on the observable related variables \(w_{\phi }\) and \(\Omega _{\phi }\) (see [2128] and a recent paper about the general property of dynamical quintessence field [30]). Since the definition of the variables x and y could vary with different scalar field models, while the meaning of \(w_{\phi }\) and \(\Omega _{\phi }\) are the same for different dark energy models and therefore the differential equations of \(w_{\phi }\) and \(\Omega _{\phi }\) are model independent. The paper is organized as follows. We firstly present the basic theoretical framework for (phantom) quintessence, (phantom) tachyon, K-essence, and general non-canonical scalar field models in Sect. 2, and we try to give the relationships between those different scalar fields in this section. We then give the three-dimensional autonomous dynamical systems for those scalar fields and switch the dynamical variables from (xy) to \((w_{\phi }, \Omega _{\phi })\) in each subsection of Sect. 3. Additionally, using the dynamical systems, we give the exact solution of \(w_{\phi }\) and \(\Omega _{\phi }\) for a special case of tachyon model when the potential is chosen to be a constant in Sect. 3.2. We show that the dynamical autonomous system of K-essence can reduce to tachyon model, and investigate another special case called kinetically driven quintessence with the lagrangian \(p(X, \phi )=f(\phi )(-X+X^2)\) detailedly in Sect. 3.3. In Sect. 3.4, we show that the dynamical autonomous system of general non-canonical scalar field can reduce to quintessence and tachyon model, respectively, for some special cases. We also studied two special cases of the so-called purely kinetic united model \(L = F(X)\) in detailed in this subsection. We try to give the cosmological implications of the three-dimensional dynamical autonomous system and present the conclusion in Sect. 4. We also raise a question about the possibility of chaotic behavior in the spatially flat single scalar field FRW cosmological models in the presence of the ordinary matter.

2 Basic framework for various scalar fields

Let us restrict ourselves to a flat universe described by the FRW metric and consider a spatially homogeneous real scalar field \(\phi \) with non-canonical kinetic energy term. The lagrangian density is given as

$$\begin{aligned} L=L(X,V) \end{aligned}$$
(1)

where L is a function of X and potential \(V(\phi )\), \(X=\frac{1}{2}\nabla _{\mu }\phi \nabla ^{\mu }\phi =\frac{1}{2}{\dot{\phi }}^2 \) for a spatially homogeneous scalar field. The pressure, energy density, and the Friedmann equations of the scalar field could easily be obtained as follows:

$$\begin{aligned} p_{\phi }=L(X, V), \end{aligned}$$
(2)
$$\begin{aligned} \rho _{\phi }=2X L_X -L, \end{aligned}$$
(3)
$$\begin{aligned} H^2=\left( \frac{\dot{a}}{a}\right) ^2=\frac{1}{3M^2_{pl}}[\rho _{\phi }+\rho _b], \end{aligned}$$
(4)
$$\begin{aligned} \dot{H}=-\frac{1}{2M^2_{pl}}[2X L_X+\gamma _b\rho _b], \end{aligned}$$
(5)

where \(8\pi G=\kappa ^2=1/M^2_{pl}\), \(\rho _b\) is the density of a barotropic fluid component with the equation of state \(p_b=w_b \rho _b=(\gamma _b-1) \rho _b\). \(\gamma _b=1\) for matter and \(\gamma _b=4/3\) for radiation. \(L_X\) is the derivative of \(L(X,\phi )\) with respect to X.

For the quintessence, general non-canonical scalar field, tachyon, and K-essence model, the pressure p and energy density \(\rho \) are, respectively,

$$\begin{aligned} p_q=\varsigma \frac{1}{2}\dot{\phi }^2-V(\phi ),\quad \rho _q=\varsigma \frac{1}{2}\dot{\phi }^2+V(\phi ), \end{aligned}$$
(6)
$$\begin{aligned} p_g=F(X)-V(\phi ),\quad \rho _g=2XF_X-F(X)+V(\phi ), \end{aligned}$$
(7)
$$\begin{aligned} p_t=-V(\phi )\sqrt{1-\varsigma \dot{\phi }^2},\quad \rho _t=\frac{V(\phi )}{\sqrt{1-\varsigma \dot{\phi }^2}}, \end{aligned}$$
(8)
$$\begin{aligned} p_K= & {} -V(\phi )F(X),\nonumber \\ \rho _K= & {} 2L_XX-L(X, \phi )=V(\phi )(F-2XF_X). \end{aligned}$$
(9)

If \(\varsigma =1\), Eqs. (6) and (8) correspond to the quintessence and tachyon scalar field. If \(\varsigma =-1\), Eqs. (6) and (8) correspond to the phantom quintessence and phantom tachyon scalar field. General non-canonical scalar field Eq. (7) can recover to (phantom) quintessence Eq. (6) if \(F(X)=\varsigma {\dot{\phi }}^2=2\varsigma X\). K-essence model Eq. (9) can recover to (phantom) tachyon Eq. (8) if \(F(X)=\sqrt{1-\varsigma \dot{\phi }^2}=\sqrt{1-2\varsigma X}\). Moreover, if the scalar field is redefined, it is demonstrated that the K-essence model described by Eq. (9) with a linear kinetic function \(F(X)=X+1\) can reduce to any quintessence model as described by Eq. (6). It means that any quintessence can be contained into K-essence frame, so each quintessence model is kinematically equivalent to a K-essence model. The authors also give the relationship between the potentials of the two models [9]. For example, the exponential potential \(V(\phi )=V_0 e^{-\lambda \phi }\) in quintessence model plays the similar role as the inverse square potential \(V(\phi )=(\frac{1}{2}\kappa \lambda \phi -c_1)^{-2}\). We will also show that the role of inverse square potential in K-essence model is very similar with exponential potential in quintessence in the next section.

3 Dynamical system of various scalar fields

In this section, we will give the dynamical system for the quintessence, tachyon, K-essence, and general non-canonical scalar field model. We will summarize the dynamical system analysis and give our comments.

3.1 Dynamical system for quintessence and phantom quintessence scalar field

For the (phantom) quintessence scalar field with lagrangian \(L=\varsigma \frac{1}{2}\dot{\phi }^2-V(\phi )\), we can define the following dimensionless variables:

$$\begin{aligned} x=\frac{\kappa \dot{\phi }}{\sqrt{6} H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3} H},\quad \lambda (\phi )=-\frac{V'}{V},\quad \Gamma (\phi )=\frac{V V''}{V'^2} \end{aligned}$$
(10)

where \(V'=\mathrm{d}V(\phi )/\mathrm{d}\phi , V''=\mathrm{d}^2V(\phi )/\mathrm{d}\phi ^2 \). The parameters \(\Gamma (\phi )\) and \(\lambda (\phi )\) of the potentials can be related with the famous slow roll parameters \(\epsilon _V\) and \(\eta _V\) (e.g., see [10]):

$$\begin{aligned} \epsilon _V(\phi )=\frac{M^2_{pl}}{2}\left( \frac{V'}{V}\right) ^2=\frac{M^2_{pl}}{2}\lambda ^2(\phi ), \end{aligned}$$
(11)
$$\begin{aligned} \eta _V(\phi )=M^2_{pl}\frac{V''}{V}=M^2_{pl}\frac{V}{V'^2}\frac{V''}{V}\frac{V'^2}{V}=M^2_{pl}\Gamma (\phi )\lambda ^2(\phi ). \end{aligned}$$
(12)

Using Eqs. (4), (5), (6), and (7), We can write down the following equations for the evolution of the (phantom) quintessence:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}N}= & {} -3x+\frac{\sqrt{6}}{2}\varsigma \lambda y^2\nonumber \\&+\frac{3}{2} x\left[ \varsigma (1-w_b)x^2+(1+w_b)(1-y^2)\right] , \end{aligned}$$
(13)
$$\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}N}=-\frac{\sqrt{6}}{2}\lambda xy+\frac{3}{2}y\left[ \varsigma (1-w_b)x^2+(1+w_b)(1-y^2)\right] , \end{aligned}$$
(14)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=-\sqrt{6} \lambda ^2(\Gamma -1)x, \end{aligned}$$
(15)

where \(N=\mathrm{ln}(a)\), a is the scale factor. \(\varsigma =1\) or \(-1\) for the quintessence and phantom quintessence model. Here we should emphasize that the system in Eqs. (13)–(15) is not a dynamical autonomous system since the parameter \(\Gamma (\phi )\) is unknown. The energy density fraction of the dark energy scalar field is

$$\begin{aligned} \Omega _{\phi }=\frac{\rho _{\phi }}{3M^2_{pl}H^2}=\varsigma x^2+y^2, \end{aligned}$$
(16)

while the equation of state of the dark energy scalar field is

$$\begin{aligned} \gamma _{\phi }=1+w_{\phi }=1+\frac{\varsigma x^2-y^2}{\varsigma x^2+y^2}=\frac{2\varsigma x^2}{\varsigma x^2+y^2}=\frac{2\varsigma x^2}{\Omega _{\phi }}. \end{aligned}$$
(17)

On the other hand, it is more convenient to rewrite the dynamical system Eqs. (13) and (14) from the dependent variables \((x, y, \lambda )\) directly to the observable quantities \((\Omega _{\phi }, \gamma _{\phi }, \lambda )\) [24]:

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }), \end{aligned}$$
(18)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}=(2-\gamma _{\phi })\left( \lambda \sqrt{3\varsigma \gamma _{\phi }\Omega _{\phi }}-3\gamma _{\phi }\right) , \end{aligned}$$
(19)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=-\sqrt{3} \lambda ^2(\Gamma -1)\sqrt{\varsigma \gamma _{\phi }\Omega _{\phi }}. \end{aligned}$$
(20)

Equations (18)–(20) can reduce to quintessence when \(\varsigma =1\) (Eqs. (3)–(5) in Ref. [25] or Eqs. (17)–(18) in Ref. [24]) and reduce to phantom quintessence when \(\varsigma =-1\) [23]. Equations (18)–(20) are very useful to study the cosmological implication of the evolution behavior of the dynamical system because the dynamical variables \(\Omega _{\phi }\) and \(\gamma _{\phi }\) are the observable quantities. For example, \(\mathrm{d} \gamma _{\phi }/\mathrm{d}N>0\) corresponds to the thawing model and \(\mathrm{d} \gamma _{\phi }/\mathrm{d}N<0\) corresponds to the freezing model of the evolution of equation of state of dark energy [25, 26].

If \(\Gamma =1\), Eq. (15) (or Eq. (20)) will disappear, Eqs. (13) and (14) (or Eqs. (18) and (19)) will become a two-dimensional dynamical autonomous system. \(\Gamma =1\) corresponds to the exponential potential \(V_0 e^{-\lambda \phi }\), which has been studied in much of the literature (e.g., see Ref. [11]). However, the system described by Eqs. (13)–(15) (or Eqs. (18)–(20)) is not a dynamical autonomous system for other potentials because the potential related parameter \(\Gamma \) is unknown. Since \(\lambda \) is a function of quintessence scalar field \(\phi \) and \(\Gamma \) is also a function of \(\phi \), then \(\Gamma \) can generally be expressed as a function of \(\lambda \). So if we consider \(\Gamma \) as a function of \(\lambda \), namely \(\Gamma =\Gamma (\lambda )\), then Eqs. (13)–(15) (or Eqs. (18)–(20)) are definitely a dynamical autonomous system, we therefore can study its properties and dynamical evolution using a phase-plane and critical points analysis. Moreover, considering \(\Gamma \) as a function of \(\lambda \) can cover many potentials beyond the exponential potential [11].

3.2 Dynamical system for tachyon and phantom tachyon scalar field

For the tachyon and phantom tachyon scalar field with lagrangian \(L=-V(\phi )\sqrt{1-\varsigma \dot{\phi }^2}\), we can define the dimensionless variables as follows:

$$\begin{aligned} x=\dot{\phi },\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3} H},\quad \lambda =\frac{V'}{\kappa V^{3/2}},\quad \Gamma =\frac{V V''}{V'^2} \end{aligned}$$
(21)

where \(V'=\mathrm{d}V(\phi )/\mathrm{d}\phi , V''=\mathrm{d}^2V(\phi )/\mathrm{d}\phi ^2 \). Using Eqs. (4), (5), (6), and (21), the evolution of (phantom) tachyon can be described in the following dynamical form [27]:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}N}=-\sqrt{3}(1-\varsigma x^2)(\sqrt{3}x+\varsigma \lambda y), \end{aligned}$$
(22)
$$\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}N}=\frac{\sqrt{3}}{2}y(\lambda xy+\frac{\sqrt{3} y^2(\varsigma x^2-\gamma _b)}{\sqrt{1-\varsigma x^2}}+\sqrt{3}\gamma _b), \end{aligned}$$
(23)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=\sqrt{3} x y\lambda ^2(\Gamma -\frac{3}{2}). \end{aligned}$$
(24)

The density parameter of tachyon field \(\Omega _{\phi }\) and the equation of state \(w_{\phi }\) are

$$\begin{aligned} \Omega _{\phi }=\frac{y^2}{\sqrt{1-\varsigma x^2}}, \end{aligned}$$
(25)
$$\begin{aligned} \gamma _{\phi }=1+w_{\phi }=\varsigma x^2. \end{aligned}$$
(26)

We can also rewrite the dynamical system Eqs. (22)–(24) from the dependent variables \((x, y, \lambda )\) directly to the observable quantities \((\Omega _{\phi }, \gamma _{\phi }, \lambda )\):

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }), \end{aligned}$$
(27)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}=-2(1-\gamma _{\phi })\left[ 3\gamma _{\phi }+\lambda \sqrt{3\varsigma \gamma _{\phi } \Omega _{\phi }}(1-\gamma _{\phi })^{\frac{1}{4}}\right] , \end{aligned}$$
(28)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=\sqrt{3\varsigma \gamma _{\phi }\Omega _{\phi }}(1-\gamma _{\phi })^{\frac{1}{4}}\lambda ^2\left( \Gamma -\frac{3}{2}\right) . \end{aligned}$$
(29)

Equations (27)–(29) are also obtained in [28, 29]. The critical points and the dynamical evolution of this system describing by Eqs. (22)–(24) or Eqs. (27)–(29) had been studied in [27].

If \(\Gamma =3/2\), Eq. (24) (or Eq. (29)) will disappear, Eqs. (22)–(23) (or Eqs. (27)–(28)) will become a two-dimensional dynamical autonomous system. For (phantom) tachyon scalar field, \(\Gamma =3/2\) corresponds to the case that the form of potential is inverse square potential. But for other potentials, the system described by Eqs. (22)–(24) (or Eqs. (27)–(29)) will not be a dynamical autonomous system any more since the potential related parameter \(\Gamma \) is unknown, and we cannot exactly analyze the evolution of universe like the inverse square potential any more. However, since \(\lambda \) is the function of tachyon field \(\phi \) and \(\Gamma \) is also the function of \(\phi \), \(\Gamma \) can be generally expressed as a function of \(\lambda \). So as the method used in [11], we can consider \(\Gamma \) as a function of \(\lambda \), \(\Gamma =\Gamma (\lambda )\), then Eqs. (22)–(24) (or Eqs. (27)–(29)) will become a dynamical autonomous system, we therefore can study its properties and dynamical evolution using the phase plane and critical points analysis. For each form of the function \(\Gamma (\lambda )\), we can figure out the detailed form of potential, so this method can cover many potentials beyond the inverse square potential [27].

For the most simple case of potential \(V(\phi )=V_0\), \(\lambda =0\), Eqs. (27)–(29) become a very simple differential equation:

$$\begin{aligned} \frac{\mathrm{d}\Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }),\quad \frac{\mathrm{d}\gamma _{\phi }}{\mathrm{d}N}=-6\gamma _{\phi }(1-\gamma _{\phi }). \end{aligned}$$
(30)

We can get the exact solution for the above equations:

$$\begin{aligned} \gamma _{\phi }=\frac{1}{1+c_1 e^{6N}},\quad \Omega _{\phi }=\frac{\sqrt{1+c_1 e^{6N}}}{\sqrt{1+c_1 e^{6N}}+c_2 (e^{6N})^{(1-\gamma _b)/2}}. \end{aligned}$$
(31)

The tachyon scalar field with this constant potential has been studied in [34]. According to the best-fit values of the parameters one obtained, we can obtain the value of the integral constant \(c_1=0.0082\) and \(c_2=0.59\) here. We know from Eq. (31), when \(N\rightarrow +\infty \), that we have \(\gamma _{\phi }\rightarrow 0\) and \(\Omega _{\phi }\rightarrow 1\); the universe will be a de Sitter-like universe filled with the tachyon scalar field.

3.3 Dynamical system for K-essence scalar field

For the K-essence scalar field with lagrangian \(L=-V(\phi )F(X)\), we define the following dimensionless variables to be the same as in Eq. (21):

$$\begin{aligned} x=\dot{\phi },\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3} H}, \quad \lambda =\frac{V'}{\kappa V^{3/2}},\quad \Gamma =\frac{V V''}{V'^2}. \end{aligned}$$
(32)

Using Eqs. (4), (5), (6), and (32), we get the following dynamical system:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}N}=\frac{-\sqrt{3}}{F_X+2XF_{XX}}\left[ \sqrt{3}xF_X-y\lambda (F-2X F_X)\right] , \end{aligned}$$
(33)
$$\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}N}=\frac{y}{2}\left[ \sqrt{3}\lambda x y+6XF_X(\gamma _b-1)y^2-3\gamma _b F y^2+3\gamma _b\right] , \end{aligned}$$
(34)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=\sqrt{3} x y\lambda ^2\left( \Gamma -\frac{3}{2}\right) . \end{aligned}$$
(35)

The density parameter of tachyon scalar field \(\Omega _{\phi }\), the equation of state \(w_{\phi }\) are

$$\begin{aligned} \Omega _{\phi }=(F-2XF_X)y^2=(F-xF')y^2, \end{aligned}$$
(36)
$$\begin{aligned} \gamma _{\phi }=\frac{xF'}{xF'-F}=-\frac{xy^2F'}{\Omega _{\phi }}, \end{aligned}$$
(37)

where \(F'=\mathrm{d} F(X)/d x\) and \(F''=\mathrm{d}^2 F(X)/dx^2\), they both are the functions of x. Using Eqs. (36) and (37), we can also rewrite the dynamical system Eqs. (33), (34), and (35) from the dependent variables \((x, y, \lambda )\) directly to the observable quantities \((\Omega _{\phi }, \gamma _{\phi }, \lambda )\):

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }), \end{aligned}$$
(38)
$$\begin{aligned} \frac{d \gamma _{\phi }}{dN}= & {} -\frac{\sqrt{3}(\lambda x y+\sqrt{3}\gamma _{\phi })[xF''(1-\gamma _{\phi })+F']}{xF''}\nonumber \\= & {} \sqrt{3}(\lambda x y+\sqrt{3}\gamma _{\phi })\left( \gamma _{\phi }-1-\frac{1}{2\Xi +1}\right) , \end{aligned}$$
(39)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=\sqrt{3}x y \lambda ^2\left( \Gamma -\frac{3}{2}\right) , \end{aligned}$$
(40)

where \(\Xi =XF_{XX}/F_X=(xF''-F')/2F'\). We have pointed out the relationship between quintessence and K-essence in the previous section. K-essence model described by Eq. (9) with a linear kinetic function \(F(X)=X+1\) can reduce to any quintessence model described by Eq. (6) [9]. So all quintessence models with any potentials can be considered as special cases of the K-essence model. Moreover, it proved that the correspondence of the exponential potential \(V(\phi )=V_0 e^{-\lambda \phi }\) in quintessence model is exactly the inverse square potential \(V(\phi )=(\frac{1}{2}\kappa \lambda \phi -c_1)^{-2}\) in K-essence model. This can explain why the dynamical system Eqs. (13)–(15) (or Eqs. (18)–(20)) of quintessence reduces to two-dimensional autonomous system for an exponential potential (\(\Gamma =1\)), while the same situation occurs for the inverse square potential (\(\Gamma =3/2\)) in tachyon and K-essence model. The author had considered another non-canonical scalar field lagrangian defined as \(L_{\phi }=Vf(B)\) with \(B=X/V\) [35]. It also includes the canonical quintessence if we choose \(f(B)=B-1\).

Let us discuss whether Eqs. (38)–(40) could be considered as an autonomous system. Firstly, we realized that the variables x and y still appear in Eq. (39) because we do not know the detailed form of the function F(X) and cannot figure out the solution of x and y. However, from Eqs. (36)–(37), we are sure that the variables x and y are functions of the variables \(\Omega _{\phi }\) and \(\gamma _{\phi }\). Since \(\Xi (=\,XF_{XX}/F_X=(xF''-F')/2F')\) is a function of x, it is also a function of \(\Omega _{\phi }\) and \(\gamma _{\phi }\). Secondly for the potential related parameter \(\Gamma \), if we consider \(\Gamma \) as a function of \(\lambda \) just similar with tachyon scalar field in Sect. 3.2, Eqs. (38)–(40) can eventually become a three-dimensional dynamical autonomous system for any K-essence models, and then we can easily study the critical points and the dynamical evolution beyond the inverse square potential.

We can take the tachyon model as a very simple example of K-essence model. For the (phantom) tachyon scalar field described in Sect. 3.2, the function F(X) should take the form of \(\sqrt{1-2\varsigma X}=\sqrt{1-\varsigma x^2}\), and then Eqs. (36)–(37) will become Eqs. (25)–(26). We can obtain the relationship as follows:

$$\begin{aligned} x=\sqrt{\varsigma \gamma _{\phi }},\quad y=\sqrt{\Omega _{\phi }}(1-\gamma _{\phi })^{\frac{1}{4}},\quad 1+\frac{1}{2\Xi +1}=2-\gamma _{\phi }. \end{aligned}$$
(41)

Putting Eqs. (41) into (38)–(40), these equations will reduce to Eqs. (27)–(29). That means the dynamical system Eqs. (38)–(40) we derived for the K-essence model are correct.

Another example is the lagrangian \(p(\phi , X)=f(\phi )(-X+X^2)\), which is proposed as a kinetically driven quintessence model [36]. The pressure and energy density are given by

$$\begin{aligned}&p(\phi , X)=f(\phi )(-X+X^2);\nonumber \\&\rho =2X\frac{\partial p}{\partial X}-p=f(\phi )(-X+3X^2). \end{aligned}$$
(42)

Using Eqs. (36) and (37) and (42), we get the following equation from Eqs. (38)–(40):

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }), \end{aligned}$$
(43)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}=\frac{(\lambda \sqrt{3(4-3\gamma _{\phi })\Omega _{\phi }}+3\gamma _{\phi })(\gamma _{\phi }-2)(3\gamma _{\phi }-4)}{3\gamma _{\phi }-8}, \end{aligned}$$
(44)
$$\begin{aligned} \frac{\mathrm{d}\lambda }{\mathrm{d}N}=\lambda ^2\sqrt{3(4-3\gamma _{\phi })\Omega _{\phi }}\left( \Gamma -\frac{3}{2}\right) . \end{aligned}$$
(45)
Fig. 1
figure 1

The evolution of the equation of state \(w_{\phi }\) with respect to N when \(w_{\phi }\) cross the Phantom Line. Figures 1 and 2 are plotted using differential equations (43)–(45)

Equations (43)–(45) completely describe the dynamical evolution of the kinetically driven quintessence. Eq. (45) will vanish when \(\Gamma =3/2\), then Eqs. (43)–(44) are a two-dimensional autonomous system which describes the dynamical evolution of the kinetically driven quintessence with the lagrangian \(p(\phi , X)=f(\phi )(-X+X^2)\) and potential \(f(\phi )=(\frac{1}{2}\kappa \lambda \phi -c_1)^{-2}\). In [37], authors obtained the two-dimensional dynamical autonomous system with the dimensionless variables (xy) for this type of K-essence model. They studied the phase-space properties and the cosmological implications of the critical points in detail. However, here we give the two-dimensional autonomous system Eqs. (43) and (44) with the observational quantities \((\Omega _{\phi }, \gamma _{\phi })\) instead of the variables (xy). We can obtain the critical points of the observational quantities \((\Omega _{\phi }, \gamma _{\phi })\) directly, so it will be more convenient to study the properties of the critical points and their cosmological implication with these observational quantities. Furthermore, if we consider \(\Gamma \) as a function of \(\lambda \), we can study the critical points and the evolution of the universe beyond the inverse square potential, just like the method used in [27].

The equation of state \(w_{\phi }(=\gamma _{\phi }-1)\) and the effective sound speed of perturbation \(c_s^2\) can be obtained from Eq. (42):

$$\begin{aligned} w_{\phi }= & {} \frac{X-1}{3X-1},\quad c_s^2=\frac{p_{,X}}{\rho _{,X}}=\frac{2X-1}{6X-1}\nonumber \\= & {} \frac{-(w_{\phi }+1)}{3(w_{\phi }+1)-8}=\frac{-\gamma _{\phi }}{3\gamma _{\phi }-8}. \end{aligned}$$
(46)

It is interesting that \(w_{\phi }\) could be larger or less than \(-1\) for this kind of K-essence model (Fig. 1). \(w_{\phi }\) is larger than \(-1\) when \(X<1/3\) (note that \(X\ge 0\) since \(X=\dot{\phi }^2/2\)) or \(X>1/2\), and less than \(-1\) when \(1/2>X>1/3\). So it is possible for the equation of state \(w_{\phi }\) crossing the Phantom Line. We have solved numerically Eqs. (43)–(45) for several different potentials and plot the evolution of \(w_{\phi }\) crossing the Phantom Line in Fig. 1. However, when the equation of state \(w_{\phi }\) crosses the Phantom Line, the effective sound speed of perturbations \(c_s^2\) will change its sign from positive to negative simultaneously (Fig. 2). For the stability with respect to the general metric and matter perturbation, the condition \(c_s^2\ge 0\) is necessary, so the background models with \(c_s^2<0\) are violently unstable and do not have any physical significance. Therefore this transition model is not realistic [38, 39]. However, we point out that not all the kinetically driven quintessence we discussed here will cross the phantom divide in the future (see Figs. 3 and 4, and their discussion). It is determined by the potentials and initial conditions. Here we just show the possibility that the kinetically driven quintessence can cross the phantom divide.

Fig. 2
figure 2

The evolution of \(c_s^2\) with respect to N. All the initial conditions are the same as Fig. 1. Figures 1 and 2 are plotted under the model of kinetically driven quintessence. The curves in yellow, blue, green, and red in Figs. 1 and 2 are for \(\Gamma =0, ~1, ~3/2, ~2\), corresponding to the potentials being of the form of \(V_0(\phi +c), ~V_0 e^{c\phi }, ~V_0(\phi +c)^{-2}, ~V_0(\phi +c)^{-1}\), respectively. The initial conditions for these four potentials are the same, \(\gamma _{\phi 0}=0.1,~\Omega _{\phi 0}=0.7\) and \(\lambda =0.1\) when \(N=0\) (at present time)

Fig. 3
figure 3

The evolution behavior of \(\gamma _{\phi }\) with respect to \(\Omega _{\phi }\). \(\lambda _0=-0.2, -0.4\) for red and black solid curve, respectively. The initial condition when \(\Omega _{\phi }=0\) is (\(\gamma _{\phi }=0\), \(\lambda =-0.2\)) and (\(\gamma _{\phi }=0\), \(\lambda =-0.4\)) for the three dashed blue, green and yellow curves around the red and black solid line, respectively

If the equation of state \(w_{\phi }\) is near \(-1\) (so \(\gamma _{\phi } \sim 0\)) for the kinetically driven quintessence model, we can drop terms of higher order in \(\gamma _{\phi }\), and further get a simple differential equation for \(\gamma _{\phi }\) with the dependent variable from N to \(\Omega _{\phi }\) from Eqs. (43) and (44):

$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d} \Omega _{\phi }}=-\frac{2\lambda }{\sqrt{3\Omega _{\phi }}(1-\Omega _{\phi })}-\frac{\gamma _{\phi }}{\Omega _{\phi }(1-\Omega _{\phi })}. \end{aligned}$$
(47)

We assume that \(\lambda \) is approximately constant (\(=\lambda _0\)) when \(\gamma _{\phi }\) is near 0, so that the above equation can be solved exactly:

$$\begin{aligned} \gamma _{\phi }=1+w_{\phi }= & {} -\frac{\sqrt{3}}{3}\lambda _0\frac{1-\Omega _{\phi }}{\Omega _{\phi }}\nonumber \\&\times \left( ln\left[ \frac{1-\sqrt{\Omega _{\phi }}}{1+\sqrt{\Omega _{\phi }}}\right] +\frac{2\sqrt{\Omega _{\phi }}}{1-\Omega _{\phi }}\right) ; \end{aligned}$$
(48)

we have chosen the boundary condition that \(\gamma _{\phi }=0\) at \(\Omega _{\phi }=0\). Equation (48) gives the general behaviors of all the kinetic driven quintessence model for all sufficiently flat potentials.

Figures 3 and 4 show how accurate the analytic result (Eq. (48)) is and they show the tiny differences among different potentials. Solid black and red curves are plotted using Eq. (48) directly, and all other curves are numerically plotted from Eqs. (43)–(45).

Fig. 4
figure 4

The evolution behavior of \(\gamma _{\phi }\) with respect to \(\Omega _{\phi }\). All the initial conditions are the same as Fig. 3 except \(\lambda _0=0.2, 0.4\). Solid curve in Figs. 3 and 4 are plotted using Eq. (48), all other curves plotted using differential equations (43)–(45) with different potentials (namely, \(\Gamma =1, 3/2, 2\))

The solid curves (red and black) in Figs. 3 and 4 show the general relationship (Eq. (48)) of \(\gamma _{\phi }\), \(\Omega _{\phi }\), and \(\lambda _0\). Around each solid curve are the dashed color curves (blue, green, and yellow), which are the numerical results of Eqs. (43)–(45) with different values of \(\Gamma \) since different values of \(\Gamma \) corresponding to different forms of the potentials: a dashed blue curve for \(\Gamma =1\) (\(~V_0 e^{c\phi }\)), a dashed green curve for \(\Gamma =3/2\) (\(~V_0(\phi +c)^{-2}\)), and a dashed yellow curve for \(\Gamma =2\) (\(V_0(\phi +c)^{-1}\)). For the dashed green curve, \(\Gamma =3/2\) and then \(\lambda \) is a constant, so it is easily understood why the difference between dashed green curve and red(black) solid curve is so small. It also demonstrated that Eq. (48) is valid when the value of \(\gamma _{\phi }\) is small. The only difference between Figs. 3 and 4 is in the opposite values of \(\lambda _0\). The initial values of \(\lambda \) when \(\Omega _{\phi }=0\) are \(-0.2\) (red) and \(-0.4\) (black) in Fig. 3 while 0.2 (red) and 0.4 (black) in Fig. 4. We find that the red solid curve and the dashed color curves around it in both Figs. 3 and 4 are more close to the Phantom Line (\(w_{\phi }=-1\)). So the smaller the initial value of \(|\lambda |\) is (i.e., the flatter the potential is), the less deviation the equation of state \(w_{\phi }\) has from \(-1\). It is interesting that, for different potentials \(~V_0 e^{c\phi }\), \(V_0(\phi +c)^{-2}\) and \(V_0(\phi +c)^{-1}\), the equation of state \(w_{\phi }\) can be larger or less than \(-1\), and also it can increase or decrease with respect to \(\Omega _{\phi }\), only depending on the initial value of \(\lambda \) (determined by the initial value of \(\phi \)). Unlike the evolution in Figs. 1 and 2, \(w_{\phi }\) plotted in Figs. 3 and 4 does not cross the Phantom Line due to the different choice of the initial conditions.

3.4 Dynamical system for general non-canonical scalar field

In the last subsection of Sect. 3, we consider the general non-canonical scalar field with a lagrangian \(L=F(X)-V(\phi )\) in a spatially flat Fridmann–Lemaitre–Robertson–Walker (FLRW) cosmology. The motion of this scalar field and the evolution of the universe are described by Eqs. (4), (5) and (7).

In order to obtain the autonomous system we define the variables as follows [40]:Footnote 1

$$\begin{aligned} x= & {} \frac{\sqrt{(2XF_X-F)\varsigma }}{\sqrt{3}M_{pl}H},\quad y=\frac{\sqrt{V}}{\sqrt{3}M_{pl}H},\sigma \nonumber \\= & {} -\frac{M_{pl}V_{\phi }}{V}\sqrt{\frac{2X}{3(2XF_X-F)\varsigma }}\mathrm{sign}(\dot{\phi }). \end{aligned}$$
(49)

Using Eqs. (7) and (49), Eqs. (4) and (5) can be expressed as the following dynamical system:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}N}= & {} \frac{3}{2}[\sigma \varsigma y^2-x(w_k+1)]\nonumber \\&+\frac{3}{2}x[(1+w_b)(1-y^2)+(w_k-w_b)\varsigma x^2], \end{aligned}$$
(50)
$$\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}N}=-\frac{3}{2}\sigma xy+\frac{3}{2} y[(1+w_b)(1-y^2)+(w_k-w_b)\varsigma x^2], \end{aligned}$$
(51)
$$\begin{aligned} \frac{\mathrm{d}\sigma }{\mathrm{d}N}= & {} -3 \sigma ^2 x(\Gamma -1)+\frac{3\sigma [2\Xi (w_k+1)+w_k-1]}{2(2\Xi +1)(w_k+1)}\nonumber \\&\times \left( w_k+1-\frac{\sigma y^2}{\varsigma x}\right) , \end{aligned}$$
(52)

where

$$\begin{aligned} w_k= & {} \gamma _k-1=\frac{F}{2XF_X-F},\quad w_{\phi }=\frac{p_{\phi }}{\rho _{\phi }}=\frac{w_k\varsigma x^2-y^2}{\varsigma x^2+y^2},\nonumber \\ \Omega _{\phi }= & {} \varsigma x^2+y^2,\quad \Xi =XF_{XX}/F_X. \end{aligned}$$
(53)

From Eq. (53), we get the following relationships:

$$\begin{aligned} x^2=\frac{\gamma _{\phi }\varsigma }{\gamma _k}\Omega _{\phi },\quad y^2=\frac{\gamma _k-\gamma _{\phi }}{\gamma _k}\Omega _{\phi }. \end{aligned}$$
(54)

Using Eqs. (53)–(54), the dynamical system Eqs. (50)–(52) can be rewritten with the variables \((\Omega _{\phi }, \gamma _{\phi }, \sigma )\) related with the observable quantities:

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }\left( 1-\Omega _{\phi }\right) , \end{aligned}$$
(55)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}= & {} 3\left[ \varsigma \gamma _{\phi }\left( \gamma _{\phi }-1-\frac{1}{2\Xi +1}\right) \right. \nonumber \\&+\left. 2\sigma \sqrt{\frac{\varsigma \Omega _{\phi }\gamma _{\phi }}{\gamma _k}}\frac{(\gamma _k-\gamma _{\phi })(\Xi +1)}{(2\Xi +1)\gamma _k}\right] , \end{aligned}$$
(56)
$$\begin{aligned} \frac{\mathrm{d}\sigma }{\mathrm{d}N}= & {} -3\sigma \left( \sigma (\Gamma -1)\sqrt{\frac{\varsigma \Omega _{\phi }\gamma _{\phi }}{\gamma _k}}+\left( \frac{1}{2\Xi +1}-\frac{\gamma _k}{2}\right) \right. \nonumber \\&\times \left. \left[ \sigma \sqrt{\frac{\varsigma \Omega _{\phi }\gamma _{\phi }}{\gamma _k}}(\frac{1}{\gamma _k}-\frac{1}{\gamma _{\phi }})+1\right] \right) . \end{aligned}$$
(57)

Beside the three variables \((\Omega _{\phi }, \gamma _{\phi }, \sigma )\), there are other quantities in Eqs. (56) and (57) which are not constant: \(\Xi , \gamma _k\) and \(\Gamma \). We know that \(\gamma _{\phi }=p_{\phi }/\rho _{\phi }+1=\gamma _{\phi }(X,\phi )=\gamma _{\phi }(\dot{\phi }, \phi )\) and \(\sigma =\sigma (\dot{\phi }, \phi )\). So generally speaking, we can obtain the expressions of \(\phi =\phi (\gamma _{\phi }, \sigma )\) and \(\dot{\phi }=\dot{\phi }(\gamma _{\phi }, \sigma )\). For the parameters \(\Xi , \gamma _k\), and \(\Gamma \), we know \(\Xi =\Xi (\dot{\phi }), \gamma _k=\gamma _k(\dot{\phi })\), and \(\Gamma =\Gamma (\phi )\), so they can be expressed as \(\Xi (\gamma _{\phi }, \sigma ), \gamma _k(\gamma _{\phi }, \sigma )\), and \(\Gamma (\gamma _{\phi }, \sigma )\). Therefore, Eqs. (55)–(57) could in principle be a dynamical autonomous system though it is actually quite complicated.

We take two simple examples to check the correctness of Eqs. (55)–(57) and study its properties. The first example: we know that the general non-canonical scalar field with lagrangian \(F(X)-V(\phi )\) will reduce to (phantom) quintessence if \(F(X)=\varsigma X\), and Eqs. (56)–(57) should reduce to Eqs. (19)–(20) (Eq. (55) has the same form as Eq. (18)). \(F(X)=\varsigma X\) makes \(\Xi =0\), \(\gamma _k=2\), and \(\sigma =\sqrt{\frac{2}{3}}\lambda \), then we found Eqs. (56)–(57) have the same form as Eqs. (19)–(20).

Another example is that, if for the potential \(V(\phi )=-1\) in K-essence and the potential \(V(\phi )=0\) in a general non-canonical scalar field, these two scalar field models will have the same form of lagrangian \(L=F(X)\), then both will reduce to the purely kinetic united model [41, 42]. In this case, \(\mathrm{d}\lambda /\mathrm{d}N=0\) and \(\mathrm{d}\sigma /\mathrm{d}N=0\), then both Eqs. (38)–(40) and (55)–(57) will reduce to a two-dimensional dynamical system as follows:

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }), \end{aligned}$$
(58)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}=3\gamma _{\phi }\left( \gamma _{\phi }-1-\frac{1}{2\Xi +1}\right) . \end{aligned}$$
(59)

For the case of the purely kinetic united model \(L=F(X)\), we know from Eqs. (49) and (53) that \(\gamma _{\phi }\) is a function of X, \(\gamma _{\phi }=\gamma _k=\frac{F}{2XF_X-F}\). In the meantime, \(\Xi \) is also a function of X because \(\Xi =XF_{XX}/F_X=(xF''-F')/2F'\). So generally speaking, \(\Xi \) can be expressed as a function of \(\gamma _{\phi }\). Then the system of Eqs. (58) and (59) could be a two-dimensional autonomous dynamical system. For some special cases of F(X), we can even get the exact solution for \(\Omega _{\phi }\) and \(\gamma _{\phi }\). We take two examples to illustrate our viewpoints.

The most simple is the case that \(\Xi \) is a constant. We know that \(\Xi =\frac{2-\alpha }{2(\alpha -1)}=\frac{XF_{XX}}{F_X}\), we can integrate and get the form for F(X):

$$\begin{aligned} F(X) =C X^{\frac{\alpha }{2(\alpha -1)}}+X_0. \end{aligned}$$
(60)

We can also get the exact solution for \(\Omega _{\phi }\) and \(\gamma _{\phi }\) from Eqs. (58)–(59):

$$\begin{aligned} \gamma _{\phi }(N)= & {} \frac{1}{c_3 e^{3\alpha N}+1/\alpha },\nonumber \\ \Omega _{\phi }(N)= & {} \frac{1}{c_4(c_3\alpha e^{3\alpha N}+1)^{-1}e^{3(\alpha -\gamma _b)N}+1}, \end{aligned}$$
(61)

where \(\alpha =1+\frac{1}{2\Xi +1}\), C, \(X_0\), \(c_3\), and \(c_4\) are the integral constants. If we set the equation of state of dark energy \(\gamma _{\phi }(0)=\gamma _0 \) and the energy density of dark energy \( \Omega _{\phi }(0)= \Omega _0\) at present \(N=0\), we then get \(c_3=\frac{\alpha -\gamma _0}{\alpha \gamma _0}\) and \(c_4= \frac{\alpha (1-\Omega _0)}{\Omega _0\gamma _0}\).

Fig. 5
figure 5

The evolution of \(\Omega _{\phi }\) with respect to N when \(F(X) =C X^{\frac{\alpha }{2(\alpha -1)}}+X_0\) and \(\gamma _b=1\), the initial condition is chosen as \(\gamma _{\phi }=0.1\) and \(\Omega _{\phi }=0.7\) at present (\(N=0\)). The color lines from the top to the bottom are plotted with \(\alpha =0.3, ~0.9, ~1.0, ~4/3, ~2.0\)

Fig. 6
figure 6

The evolution of \(\Omega _{\phi }\) with respect to N when \(F(X) =C X^{\frac{\alpha }{2(\alpha -1)}}+X_0\). The initial condition and the values of the parameter \(\alpha \) are the same as Fig. 5 except \(\gamma _b=4/3\). Figures 5 and 6 are plotted using Eq. (61)

For the cosmic evolution in very early time, N is negative and |N| is very large, we get \(\gamma _{\phi }\approx \alpha \) from Eq. (61), and then the energy density of scalar field will behave as \(\rho _{\phi }\sim a^{-3\alpha }\). For the cosmic evolution in late time, N is positive and very large, \(\gamma _{\phi }\approx 0\) and \(\Omega _{\phi }\approx 1\), scalar field behaves as the cosmological constant. Noted that there is only a kinetic term in the lagrangian which gives the cosmological constant solution. Moreover, we know from Eq. (61) that \(\Omega _{\phi }\) might not be 0 in the very early time. Its value depends on the value of \(\alpha \) comparing with the value of \(\gamma _b\) (\(\gamma _b=1\) for matter and \(\gamma _b=4/3\) for radiation). We have plotted the evolution of \(\Omega _{\phi }\) with respect to N for different \(\alpha \) and different \(\gamma _b\) in Figs. 5 and 6, and it is shown that the value of \(\Omega _{\phi }\) in the very early time could be 1, 0 or a positive constant which is less than 1. When \(N \rightarrow -\infty \), we can get \(\Omega _{\phi }(N) \approx 1/(c_4 e^{3(\alpha -\gamma _b)N}+1)\) from Eq. (61). Then \(\Omega _{\phi }\rightarrow 1\) for \(\alpha >\gamma _b\), \(\Omega _{\phi }\rightarrow 1/(c_4+1)\) for \(\alpha =\gamma _ b\) and \(\Omega _{\phi }\rightarrow 0\) for \(\alpha <\gamma _b\). If \(\Omega _{\phi }>0\) in very early time, it would be very interesting to investigate its impact on the evolution history of the early universe.

The second case is the following lagrangian:

$$\begin{aligned} F(X)=A_1\sqrt{X}-A_2 X^{\beta } \end{aligned}$$
(62)

where \(A_1\), \(A_2\), and \(\beta \) are constants. This form of F(X) was proposed in [41, 43]. From Eq. (53), we get \(\Xi =(2\beta -\gamma _{\phi })/2\gamma _{\phi }\), then the dynamical system of Eqs. (58)–(59) becomes the following equations:

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}= & {} 3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }),\end{aligned}$$
(63)
$$\begin{aligned} \frac{\mathrm{d} \gamma _{\phi }}{\mathrm{d}N}= & {} 3\gamma _{\phi }\left( \frac{2\beta -1}{2\beta } \gamma _{\phi }-1\right) . \end{aligned}$$
(64)

Solving the above differential equations, we obtain the following exact solution for \(\gamma _{\phi }\) and \(\Omega _{\phi }\):

$$\begin{aligned} \gamma _{\phi }(N)= & {} \frac{1}{c_5 e^{3N}+1/\eta },\nonumber \\ \Omega _{\phi }(N)= & {} \frac{1}{c_6(c_5\eta e^{3N}+1)^{-\eta }e^{3(\eta -\gamma _b)N}+1} \end{aligned}$$
(65)

where \(c_5\) and \(c_6\) are the integral constants, \(\eta =\frac{2\beta }{2\beta -1}\) is also a constant. Equation (65) is very similar but a little different from the evolution described in Eq. (61):

$$\begin{aligned} \gamma _{\phi }(N)= & {} \frac{1}{c_3 e^{3\alpha N}+1/\alpha },\\ \Omega _{\phi }(N)= & {} \frac{1}{c_4(c_3\alpha e^{3\alpha N}+1)^{-1}e^{3(\alpha -\gamma _b)N}+1}. \end{aligned}$$

For the cosmic evolution in very early time, N is negative and |N| is very large, we get \(\gamma _{\phi }\approx \eta =2\beta /(2\beta -1)\) from Eq. (65). The scalar field will mimic the evolution of matter with zero pressure \(\rho \sim a^{-3}\) in the limit of \(\beta \rightarrow \infty \). For the cosmic evolution in late time, N is very large, \(\gamma _{\phi }\approx 0\) and \(\Omega _{\phi }\approx 1\), the scalar field behaves as the cosmological constant. So it is the second case that a lagrangian without a potential term gives the cosmological constant solution. Moreover, we know from Eq. (65) that \(\Omega _{\phi }\) might not be 0 in the very early time. It depends on the value of \(\eta \) comparing with the value of \(\gamma _b\) (\(\gamma _b=1\) for matter and \(\gamma _b=4/3\) for radiation). When \(N \rightarrow -\infty \), we can get \(\Omega _{\phi }(N) \approx 1/(c_6 e^{3(\eta -\gamma _b)N}+1)\) from Eq. (65). Then \(\Omega _{\phi }\rightarrow 1\) for \(\eta >\gamma _b\), \(\Omega _{\phi }\rightarrow 1/(c_6+1)\) for \(\eta =\gamma _b\) and \(\Omega _{\phi }\rightarrow 0\) for \(\eta <\gamma _b\). In this case, it will be interesting to investigate the impact of the scalar field on the evolution history of the early universe. We have plotted the evolution of \(\Omega _{\phi }\) with respect to N for different \(\eta \) and different \(\gamma _b\) in Figs. 7 and 8, and it is shown that the value of \(\Omega _{\phi }\) in the very early time can actually be 1, 0 or an arbitrary positive constant which is less than 1.

4 Cosmological implications and conclusion

The main purpose of the paper is not to analyze the dynamical behavior about different scalar fields in detail, so we will not investigate the detailed critical points and their stable properties for each dynamical system. What we want to focus on is the dynamical system itself.

4.1 Variables (xy) vs. observable quantities \((\gamma _{\phi }, \Omega _{\phi })\)

The dynamical variables (xy) in the previous papers concern the scalar field and its first derivative. Though for some scalar field models (quintessence or phantom quintessence), the combination of variables (xy) has a certain meaning, namely, \(y^2 \pm x^2=\Omega _{\phi }\), generally speaking, these dynamical variables have no direct cosmological meaning, and the autonomous system for \(\mathrm{d}x/\mathrm{d}N\) and \(\mathrm{d}y/\mathrm{d}N\) also varies from different models. The authors studied the dynamical properties of the scalar field based on the quantities \((\gamma _{\phi }, \Omega _{\phi })\) instead of the variables (xy) [2128, 30]. It is more convenient if we change the variables from (xy) to the observable quantities \((\gamma _{\phi }, \Omega _{\phi })\). Firstly, \((\gamma _{\phi }, \Omega _{\phi })\) are directly related to the observable quantities and also have to do with the properties of dark energy. Analyzing the system based on \((\gamma _{\phi }, \Omega _{\phi })\), we can figure out how the equation of state of dark energy \(w_{\phi }\) and the density parameter \(\Omega _{\phi }\) evolve. Secondly, though the forms of the autonomous system for \(\mathrm{d}x/\mathrm{d}N\) and \(\mathrm{d}y/\mathrm{d}N\) are completely different for different models, it is quite interesting that the function for \(\mathrm{d}\Omega _{\phi }/\mathrm{d}N\) has the same expression in the quintessence, tachyon, K-essence, and general non-canonical scalar field model (i.e., Eqs. (18), (27), (38), (55)) as follows:

$$\begin{aligned} \frac{\mathrm{d} \Omega _{\phi }}{\mathrm{d}N}=3(\gamma _b-\gamma _{\phi })\Omega _{\phi }(1-\Omega _{\phi }). \end{aligned}$$
(66)
Fig. 7
figure 7

The evolution of \(\Omega _{\phi }\) with respect to N when \(F(X)=A_1\sqrt{X}-A_2 X^{\beta }\) and \(\gamma _b=1\), initial condition is chosen as \(\gamma _{\phi }=0.1\) and \(\Omega _{\phi }=0.7\) at present  (\(N=0\)). The color lines from the top to the bottom are plotted with \(\eta =0.3, ~0.9, ~1.0, ~4/3, ~2.0\)

Fig. 8
figure 8

The evolution of \(\Omega _{\phi }\) with respect to N when \(F(X)=A_1\sqrt{X}-A_2 X^{\beta }\). The initial condition and the values of the parameter \(\eta \) are the same as Fig. 7 except \(\gamma _b=4/3\). Figures 7 and 8 are plotted using Eq. (65)

The only difference is the form of the function \(d\gamma _{\phi }/dN\). In fact, it is well known that Eq. (66) holds for all non-coupled dark energy models as long as they satisfy the following equations:

$$\begin{aligned} H^2=\frac{1}{3M^2_{pl}}(\rho _b+\rho _{de}),\quad \dot{\rho }_i+3H(p_i+\rho _i)=0, \end{aligned}$$
(67)

where the subscript i denotes each energy component such as dark energy, matter or radiation. If we set \(\gamma _{de}=w_{de}+1=p_{de}/\rho _{de}+1\) and \(\Omega _{de}=\rho _{de}/3M^2_{pl}H^2\), Eq. (66) can be derived from Eq. (67). So the dark energy density parameter \(\Omega _{\phi }\) obeys the same differential equation (namely Eq. (66)) independent on the scalar field models considered whenever the dark energy is uncoupled in GR frame. For example, the authors got the same equation even for the purely kinetic coupled gravity model which modified the standard general relativity action through the addition of a coupling between functions of the metric and kinetic terms of a free scalar field [44] (Eq. (29) in this paper and \(\gamma _b\) is taken as 1).

We can conclude from Eq. (66) that there are only three possible destinies (three types of critical points) for \(\Omega _{de}\) in these models, namely \(\Omega _{de}=0\), \(\Omega _{de}=1\), and the case \(\gamma _{de}=\gamma _b\) where the value of \(\Omega _{\phi }\) is determined by the other equation in the dynamical system. The cases of \(\Omega _{de}=0\) or \(\Omega _{de}=1\) are completely opposite destinies, corresponding to the universe completely dominated only by the scalar field or by the barotropic fluid. Generally speaking, we can obtain \(0<\Omega _{de}<1\) for the case of \(\gamma _{de}=\gamma _b\). However, for this scaling solution, the equation of state of the dark energy \(w_{de}\) is the same as the equation of state of the barotropic fluid \(w_b\), so there is no accelerating expansion. Since the observation suggested that we are living in an accelerated expanding universe with \(\Omega _{de}\sim 0.7\), none of these three destinies corresponds to the present universe we observed. This could be considered as a clue for the possibility of the interaction between dark energy and other barotropic fluids (see [45] for such model) if we want to solve or at least alleviate the cosmological coincidence problem without fine-tunings. This result is valid for not only all the non-coupled dark energy models, but also for many modified gravity models as long as the energy density and the pressure of dark energy or effective dark energy satisfy the continuity equation Eq. (67). However, we should emphasize that our result is not new, there is much work on the study of the interacting dark energy model [21, 22, 3133].

4.2 Two-dimensional vs. three-dimensional dynamical autonomous system

Another important thing we want to emphasize is that it is more reasonable and more scientific to investigate the dynamical behaviors of a dark energy model under the three-dimensional autonomous system rather than the two-dimensional system. Firstly, the two-dimensional dynamical autonomous system is just a specific case when the potential takes a special form. If we want to completely study the general dynamical properties of a dark energy model, we need to study the system beyond a special potential. Then we can find more critical points than the ones found in a two-dimensional system. We therefore are able to analyze which critical points are possessed by a class of dark energy models and which ones exist only due to the concrete potentials. The method of studying the three-dimensional dynamical autonomous system beyond one special potential originally was meant for quintessence [11, 46] and then developed for other dark energy models [1220, 45]. Here we extend this method to the more general scalar field models in Sect. 3. Secondly, more stable attractors can be found in terms of a three-dimensional autonomous system. For example, a new critical point is found only in a three-dimensional dynamical system of power-law kinetic quintessence, which corresponds to the dark energy dominated universe (\(\Omega _\phi =1\)) where power-law kinetic quintessence behaves as a cosmological constant with the sound speed \(c_s^2\) being 0 [65]. Thirdly, from the viewpoint of chaos theory, the dynamical properties of the three-dimensional autonomous system is more fruitful than the two-dimensional system. According to the Poincaré–Bendixson theorem, chaos does not exist in any two-dimensional autonomous dynamical system [47, 48] but could be possible in a three-dimensional autonomous dynamical system. For a number of three-dimensional systems, such as the famous three-dimensional Lorenz equations which is a model describing the atmospheric convection [49], there exists chaos for certain values of the parameters.

4.3 Stable attractors vs. chaotic behaviors

The study of chaotic dynamics in cosmological models has a long history. Chaotic properties have been reported in spatially closed scalar field FRW cosmological models [5055], spatially flat FRW cosmological model with two or more scalar fields [56, 57], Bianchi IX universe [58, 59], Bianchi I universe [60], and the mixmaster universe [61]. It would be very interesting and also a big challenge for the theoretical study of dark energy if the dynamical systems we consider here (i.e., Eqs. (18)–(20), (27)–(29), (38)–(40), and (55)–(57)) have chaotic properties. Then the evolution of \(\Omega _{\phi }\) and \(\gamma _{\phi }\) will be very sensitive to the initial condition, and therefore predicting their evolution in the future becomes totally impossible. However, it is proved that there is no chaotic behavior in spatially flat single scalar field FRW cosmological models [62, 63]. Since for the spatially flat case with \(k=0\), the dynamical system can be described by a three-dimensional autonomous system with a set of variables \((H, \phi , \dot{\phi })\) under a Hamiltonian constraint, the dynamical system is actually a two-dimensional autonomous system (a and \(\dot{a}\) appear only in the combination \(H=\dot{a}/ a\)) [64]. We know that for the two-dimensional autonomous systems, there are not enough degrees of freedom for chaos to exist, so this proved the non-chaotic dynamics in the spatially flat scalar field FRW cosmological model. However, we note that this result is obtained in the absence of matter and radiation. This may be the case in the very early time when our universe is undergoing an inflation era and is completely dominated by the scalar field. However, for the study of dark energy of late-time cosmic acceleration, the component of matter is comparable with the density of dark energy and should not be ignored when we investigate the dynamical behavior of the scalar field. In the presence of matter, the scale factor a will reappear in the dynamical system beside the variables \((H, \phi , \dot{\phi })\), and then the system cannot be reduced to two-dimensional autonomous dynamical system any more. So here we argue that, beside the ordinary attractors (such as a dark energy dominated solution, a de Sitter like solution, and a scaling solution), it is still possible for the chaotic behavior in spatially flat single scalar field FRW cosmological models in the presence of matter. It is very like the case of spatially non-flat (\(k\ne 0\)) single scalar field FRW cosmological models where the dynamical system cannot be reduced to two-dimensional autonomous dynamical system too. What we argued here is also supported by the equations in Sect. 3 (i.e., Eqs. (18)–(20), (27)–(29), (38)–(40) and (55)–(57)), which described three-dimensional autonomous dynamical systems. However, we are not sure whether there truly exists chaotic behavior in spatially flat scalar field FRW cosmological models now; to find the chaotic behavior is beyond the scope of this paper, and it should be investigated in the future.