Dynamics of the universe with disformal coupling between the dark sectors

We use the dynamical analysis to study the evolution of the universe at late time for the model in which the interaction between dark energy and dark matter is inspired by disformal transformation. We extend the analysis in the existing literature by supposing that the disformal coefficient depends both on the scalar field and its kinetic terms. We find that the dependence of the disformal coefficient on the kinetic term of scalar field leads to two classes of the scaling fixed points that can describe the acceleration of the universe at late time. The first class exists only for the case where the disformal coefficient depends on the kinetic terms. The fixed points in this class are saddle points unless the slope of the conformal coefficient is sufficiently large. The second class can be viewed as the generalization of the fixed points studied in the literature. According to the stability analysis of these fixed points, we find that the stable fixed point can take two different physically relevant values for the same value of the parameters of the model. These different values of the fixed points can be reached for different initial conditions for the equation of state parameter of dark energy. We also discuss the situations in which this feature disappears.


Introduction
The observed cosmic acceleration at late time is the one of the most important mystery in the universe [1,2]. This phenomena may be explained by introducing unknown form of energy to govern the dynamics of the late-time universe [3,4]. For the simplest case, this unknown form of energy, dubbed dark energy, is supposed to be in the form of evolving scalar fields. In general, viable dark energy models should have a mechanism to alleviate the coincidence problem, which is the problem why the energy density of dark energy and a e-mail: khampheek@nu.ac.th b e-mail: stharporns57@email.nu.ac.th matter are comparable in magnitude at present although they evolve independently throughout the whole evolution of the universe [5][6][7]. A possible assumption for alleviating the coincidence problem is based on the introduction of the interaction between dark energy and dark matter. Various phenomenological form of the interaction between dark energy and dark matter have been proposed and investigated in literature [8][9][10][11]. Interestingly, it has recently been shown that models of dark energy in which the interaction between dark energy and dark matter is assumed can satisfy the bound on the Hubble parameter at redshift 2.34 from BOSS data, while the Λ CDM model predicts too large Hubble parameter at this redshift [12].
In addition to the phenomenological models of the interaction between dark energy and dark matter, the models of the interaction between dark energy and dark matter can be constructed from the frame transformation of the theory of gravity. Applying the conformal transformation to some classes of scalar-tensor theories, one obtains the coupling terms between dark energy and dark matter in the Einstein frame in which the gravity action takes the Einstein-Hilbert form [13][14][15]. The cosmological consequences of the the interaction between dark energy and dark matter due to the conformal transformation have been investigated in [16,17]. However, in order to transform general scalar-tensor theories to Einstein frame, we need transformations that are more general than the conformal transformation. It has been shown that subclasses of the GLPV theory which is the generalization of the Horndeski theory can be transformed to the Einstein frame using the disformal transformation defined as [18][19][20] g µν = C(φ )g µν + D(φ , X)φ ,µ φ ,ν , where subscript , denotes partial derivatives, and X = − 1 2 ∂ µ φ ∂ µ φ is the kinetic energy of scalar field. Here, the conformal coefficient C depend only on scalar field φ , while the disformal coefficient D can depend both on φ and its kinetic term X. In the case where D depends only on the scalar field, the above transformation provides relations among some pieces of Lagrangian in the Horndeski theory, but cannot generate a piece of the GLPV Lagrangian from the Horndeski theory [21]. However, it has been shown that if C also depends on X, the application of the transformation to GLPV action can generate terms that do not belong to the GLPV theory [19] and therefore these terms might be a cause of Ostrogradski's instibilities in the theories. Nevertheless, according to the discussion in [22], the Ostrogradski's instabilities can be eliminated by hidden constraint in some cases.
The cosmological consequences of the interaction between dark energy and dark matter due to the disformal transformation, called disformal coupling between dark energy and dark matter, has been studied in various aspects for the case where the conformal and disformal coefficients depend on the field φ only. It has been shown in [23] that the disformal coupling between dark energy and dark matter leads to a new stable fixed point compared with the case of conformal coupling, and the cosmological parameters at this fixed point can satisfy the observational bounds. The metric singularity of the new fixed point found in [23] presents the phantom behaviour in the Jordan frame [24]. The influences of the disformal coupling on the observational quantities such as the CMB and matter power spectra have been investigated in [25][26][27]. In this work, we study the disformal coupling between dark energy and dark matter for more general case where the disformal coefficients depends on both φ and its kinetic terms X. The physical motivation for such disformal coupling is related to the frame transformation among the general scalar-tensor theories presented above. Our aims are to study how the kinetic-dependent disformal coupling influences the evolution of the universe at late time by finding and analyzing the physically relevant fixed points of the model, rather than search for all possible fixed points of the model. We will show in the following sections that there are features arising only for the case where the disformal coefficient depends on both φ and X.
In section 2, the evolution equations for dark energy and dark matter with disformal coupling are presented in the covariant form. The autonomous equations for this model of dark energy are computed in section 3, and the evolution of the late-time universe is studies using the dynamical analysis in section 4. The conclusions are given in section 5.

Disformal coupling between dark energy and dark matter
In this section, we derive the disformal coupling between scalar field and matter arising from the general disformal transformation defined in eq. (1). From the metric in eq. (1), we havē In order to study coupling between dark energy and dark matter due to disformal transformation, we suppose that the field φ in the disformal transformation plays a role of dark energy, and therefore the interaction between dark energy and dark matter can occur when the Lagrangian of dark matter depends on metricḡ µν defined in eq. (1). Thus we write the action for gravity in terms of metric g µν and write the action for the dark matter in terms ofḡ µν as is the Lagrangian of the scalar field, V (φ ) is the potential of the scalar field, L c is the Lagrangian of dark matter and L M is the Lagrangian of ordinary matter including baryon and radiation. Varying this action with respect to g αβ , we get where G αβ is the Einstein tensor computed from g µν , and the energy momentum tensor for scalar field and matter are defined in unbarred frame as Using these definitions of the energy momentum tensor and the conservation of the energy momentum tensor of ordinary matter, we have ∇ α (T αβ φ + T αβ c ) = 0 due to ∇ α G αβ = 0. However, we see that the energy momentum tensors of dark energy and dark matter do not separately conserve because the Lagrangian of dark matter depends on field φ . On the other hand, sinceḡ αβ does not depend on ψ, the energy momentum tensor of dark matter is conserved in the barred frame such as where∇ α is defined from barred metric, and the energy momentum tensor in the barred frame is related to that in the unbarred frame defined in eq. (6) as To compute the interaction terms between scalar field and dark matter, we vary the action with respect to φ as One can show that where ; denotes the covariant derivative and V ,φ = ∂V /∂ φ , and Combining eq. (10) with the above equation, we obtain the evolution equation for scalar field φ , Multiplying the above equation by φ ,λ , we get In order to write the barred quantities in the interaction term Q in terms of unbarred quantities, we write eq. (8) as where J ≡ √ −ḡ/ √ −g. Hence, we get and therefore These relations yield Inserting eqs. (17), (18) and (19) into eq. (12) we can write Q as where The form of this interaction terms can reduce to that in [25,27] when λ 3 = 0.

Evolution equations for the FLRW universe
We now compute the evolution equations for all matter components in the spatially flat FLRW universe. Using the perfect fluid model for radiation and matter as well as the FLRW line element in the form the component (00) of eq. (4) yields where a prime denotes derivative with respect to conformal time τ, ρ r , ρ b and ρ c are the energy density of radiation, baryon and dark matter respectively. Furthermore, the interaction terms Q in eq. (20) becomes From now on, we will use X ≡ (φ ′ ) 2 /(2a 2 ). Inserting this expression for the interaction terms into eq. (13), we can where Substituting this expression for ρ ′ c into eq. (25) and using eq. (12), we get where and The evolution equations for ρ r and ρ b can be computed from the conservation of their energy yielding

Autonomous equations
In order to study the evolution of the universe for the disformal coupled model of dark energy, we analyze solutions for the evolution equations presented in the previous section using dynamical analysis. For concreteness, we derive the autonomous equations using the conformal coupling, disformal coupling and the scalar field potential of the form where λ 1 , λ 2 , λ 3 , λ 4 and C 0 are the dimensionless constant parameters, while M and M v are the constant parameters with dimension of mass. Here, we extend the analysis in the literature by supposing that the disformal coefficient D also depends on the kinetic term X through X λ 3 which is the simplest extension. Using the following dimensionless variables, we can write eqs. (27) and (30) in the form of autonomous equations as where N ≡ ln a and The evolution of the universe is completely described by these autonomous equations and the constraint equation which is obtained from eq. (24) as In order to derive the above autonomous equations, we also use which can be obtained by differentiating the constraint equation with respect to N. From the above equation, we see that where t is the cosmic time.

Dynamical analysis
Here, we concentrate on dynamics of the universe at late time, so that we ignore the contribution from radiation density in the autonomous equations. Moreover, we are mainly interested in the physical fixed points that correspond to the acceleration of the universe at late time.
The fixed points of the autonomous equations can be obtained by setting the LHS of eqs. (33) -(36) to zero and solving the resulting polynomial equations for x 1 , x 2 , x 3 and Ω b . The obtained solutions are the fixed points of the system denoted by variables with subscript f , e.g., x 1 f , x 2 f and x 3 f are the fixed points for x 1 , x 2 and x 3 respectively. It follows from eqs. (41) and (36) that the fixed points which correspond to the acceleration of the universe can exist if x 2 f = 0 and Ω b f = 0. When x 2 f = 0, eq. (34) gives the following relation for the fixed points: Inserting this relation into eq. (35), supposing that x 1 f = 0 and using the fact that both dx 3 /dN and dx 1 /dN vanish at fixed point, we get Applying eqs. (42) and (43) to eq. (33), we can compute the fixed points for the case where both x 1 f and x 2 f do not vanish. This case correspond to the scaling solution which will be analyzed in detail in section 4.2.
We note that the relation in eq. (43) is derived by supposing x 1 f = 0. In the case x 1 f = 0, eq. (42) gives x 2 f = 1, i.e., this fixed point is the potential dominated solution. We will consider this case in detail in the next section.

Potential dominated solution
The potential dominated solution corresponds to the fixed point (x 1 f , x 2 f ) = (0, 1). Substituting this fixed point into eq. (33) and setting Ω b = 0, we get λ 4 = 0 at the fixed point. Setting λ 4 = 0, eq. (33) yields lim Substituting this relation into eq. (35), we obtain This implies that the fixed point that corresponds to the potential dominated solution occurs in two situations. The first is the situation where the disformal coefficient is much smaller than the conformal coefficient, i.e., x 3 f = 0 The second is the situation where the disformal coefficient do not depend on X, i.e., λ 3 = 0. This result can be easily understood by noting that the field φ is nearly constant in time within the potential dominated regime, so that the ratio of disformal coefficient to conformal coefficient nearly vanish if the disformal coefficient depends on kinetic term of the scalar field. Performing the usual stability analysis, one can show that the fixed point for the potential dominated solution is stable for λ 3 ≥ 0.

Scaling and field dominated solutions
We now consider the case where both x 1 f and x 2 f do not vanish. In our consideration, x 2 1 f + x 2 f ≤ 1, i.e., Ω c ≥ 0 at the fixed points, so that these fixed points correspond to scaling and field dominated solutions.
According to eq. (43), the existence of the fixed points requires x 3 f = 0 or Since x 3 f = 0 implies that disformal coefficient vanishes at the fixed point, this fixed point is the conformal scaling solution. Hence, the case x 3 f = 0 corresponds to the disformal scaling solutions, in which the condition given in eq. (46) is required for the existance of fixed points. We will consider these fixed point in the following sections.

Conformal scaling solutions
Substituting eq. (42) into eq. (33) and then setting x 3 f = 0, we obtain the following fixed points: The first fixed point is actually the scalar field dominated solution because x 2 1 f + x 2 f = 1, while the second fixed point is the scaling solution. Since Ω b always vanishes at the fixed point in our consideration, we will perform stability analysis by linearizing only eqs. (33) -(35). The inclusion of the evolution equation for the baryon density given in eq. (36) will give rise to additional eigenvalue, µ = 3x 2 1 f − 3x 2 f , which is negative for the fixed points corresponding to accelerating universe. To simplify our analysis, we will write the eigenvalues of the fixed points in terms of the density parameter Ω d and equation of state parameter w d of dark energy at fixed point. In terms of the dynamical variables defined in eq. (32), the quantities Ω d and w d at fixed point can be written as where Ω d f and w d f are the value of Ω d and w d at fixed points respectively. Inserting the fixed points from eqs. (47) and (48) into the above equations, we respectively get , for the fixed points in eq. (47) , (50) for the fixed points in eq. (48) .
The eigenvalues µ ≡ (µ 1 , µ 2 , µ 3 ) for the fixed points in eqs. (47) and (48) are respectively given by where γ f ≡ 1 + w d f and The stabilities of the fixed points in eqs. (47) and (48) have been already discussed in literature [16,28], so that we will not consider in detail here. However, we would like to check whether the values for parameters λ 4 , λ 3 , λ 2 , λ 1 can imply the evolution of the universe at late time using the dynamical analysis. Before considering more complicate fixed points in the next sections, let us start with the potential dominated solution discussed in section 4.1 and conformal scaling solutions given in eqs. (47) and (48). In the calculation for the conformal scaling solutions, we suppose that x 1 f = 0, so that w d f > −1 and therefore λ 4 = 0 is not allowed in eqs. (50) and (51). This implies that the universe will evolve towards the potential dominated solution, corresponding to the De Sitter expansion, at late time if λ 4 = 0 and λ 3 ≥ 0.
In the case λ 4 < 0, the universe evolves towards the stable fixed point (w d , Ω d ) = (w d f , 1) if λ 1 < −6w d f / 3γ f and λ 3 as well as λ 2 are suitably chosen, e.g. λ 3 ∼ λ 2 ∼ O(1). Similarly, for the case λ 4 > 0, the universe will evolve towards the stable fixed point (w d , Ω d ) = (w d f , 1) if λ 1 > 6w d f / 3γ f . Here, w d f can be specified by λ 4 through eq. (50). However, if λ 4 and λ 1 satisfy eq. (51) and λ 3 as well as λ 2 are suitably chosen, the universe will reach the stable fixed point (w d , Ω d ) = (w d f , Ω d f ) at late time, where w d f and Ω d f are related to λ 1 and λ 4 through eq. (51). We note that if we set λ 1 , λ 2 , λ 3 and λ 4 such that eq. (46) is satisfied, the first eigenvalue in eqs. (52) and (53) will be zero. Consequently, one can show that these fixed points are saddle points, and therefore the universe will finally evolve towards the disformal scaling solutions discussed in the next section. Hence, in this section, we will consider only the cases where the relation in eq. (46) is not satisfied.
It is interesting to check whether the fixed points in eqs. (47) and (48) can be stable for the same value of λ 1 , λ 2 , λ 3 and λ 4 , i.e., the universe has two possible stable fixed points for the same value of parameters. Let us suppose that w d f for the fixed point in eq. (47) is w d f 1 , so that eq. (50) gives λ 4 = ∓ 3(1 + w d f 1 ). We then set (w d f , Ω d f ) = (w d f 2 , Ω d f 2 ) for the fixed point in eq. (48), and use eq. (51) to show that for this fixed point. In the case where λ 4 computed from both fixed points are the same, we can write Ω d f 2 in terms of w d f 1 and w d f 2 as Applying a simple numerical calculation to the above equation, it can be checked that both values of Ω d f 2 are unphysical unless −0.99 ≤ w d f 2 < w d f 1 ≤ 1. Hence, if this condition on w d f 1 and w d f 2 is satisfied and λ 1 as well as λ 4 satisfy eq. (51), the universe will evolve towards the fixed point (w d f , Ω d f ) = (w d f 2 , Ω d f 2 ) because this fixed point is stable. The universe will not evolve to fixed point given in eq. (47) because this fixed point becomes unstable in this situation. To show that the fixed point in eq. (47) is unstable in this case, we compute λ 1 from eq. (51) by setting , Ω d f 2 ) and then insert the result into the third eigenvalue in eq. (52). Therefore we get for the case −0.99 ≤ w d f 2 < w d f 1 ≤ 1, Ω d f 2 satisfies eq. (55) and 0 < Ω d f 2 < 1. Nevertheless, if λ 1 and λ 4 do not satisfy eq. (51), eq. (56) is not valid and hence the universe can evolve to the stable fixed point (w d , Ω d ) = (w d f , 1) associated with eq. (47).

Disformal scaling solutions
In the case where the relation in eq. (46) is satisfied, the fixed points at which the disformal coefficient does not vanish, i.e.,  Table 1 The fixed points for the disformal scaling solutions. For these fixed points, the parameters λ 1 and λ 3 are arbitrary, while λ 2 and λ 4 are replaced by eqs. (46) and (57) respectively.
and Ω d f by inserting eq. (42) into eq. (49) and solving the resulting equations as This implies that we can choose x 1 f and λ 4 by fixing w d f and Ω d f . Substituting these relations into eq. (42), we obtain Since the above equation is computed from eq. (42), this equation is a consequence of the vanishing of the LHS of eq. (34). Substituting relations in eqs. (46) and (42) into eq. (33), we obtain the polynomial equation which has degree 2 in x 3 f and degree 6 in x 1 f . Since it is not easy to solve this equation for x 1 f that lies inside the physical phase space, we instead solve this equation for x 3 f , and write x 1 f and λ 4 in the solutions in terms of w d f , Ω d f using eq. (57) as where x 3 f 1 and x 3 f 2 are the solutions for the polynomial equation of the fixed points. From the above calculations, we conclude that there are two classes of the fixed points for the disformal scaling solutions shown in table (1). It follows from eq. (59) that x 3 f 1 → ∞ when λ 3 = 0, implying that the class I of fixed point does not exist for the case where the disformal coefficient does not depend on the kinetic terms of scalar field. It can be checked that for the case λ 3 = 0, the fixed points belonging to the class II are the fixed points discussed in [23]. According the result in [23], one of the eigenvalues for each fixed points in the class II is zero for the case where λ 3 = 0. Hence, to simplify our analysis, we check whether the fixed points in the both classes have zero eigenvalue. We compute the metric M i j ≡ ∂ E i /∂ x j , where the indices i and j run from 1 to 3 and E i is the RHS of eqs. (33) -(35) respectively. Evaluating this matrix at the fixed points, we get Inserting the fixed points in table (1) into this matrix, we get det M = 0 for all fixed points which suggests that one of the eigenvalues for each fixed points vanishes. Therefore, to determine the stability of the fixed points, we have to go beyond linear analysis. However, as discussed above for the case of disformal solution, we have only two relations among x 1 f , x 2 f and x 3 f which are not sufficient for writing the fixed points completely in terms of the parameters of the model. Hence, is interesting to check whether the constraint equation given in eq. (46) can imply the constraint equation for the dynamical variables x 1 , x 2 and x 3 . Using the relation in eq. (46) and the definitions of x 1 , x 2 and x 3 , we get where r 0 is constant which controls magnitude of x 3 for given x 1 and x 2 . Hence, the third relation among x 1 f , x 2 f and x 3 f can be constructed from constraint given in eq. (46) by introducing other constant parameter r 0 . The above constraint equation suggests that for a given r 0 , x 3 can be computed from x 1 and x 2 , so that the dimension of the phase space is reduced [23]. Substituting this relation and the relation in eq. (42) into eq. (33), we obtain polynomial equation degree 6 + 4λ 3 of x 1 f . Similar to the previous analysis, instead of solve this equation for x 1 f , we solve this equation for r 0 , so that we get two expressions for r 0 : Moreover, we also use the relations in eq. (57) to write the above expressions in terms of w d f and Ω d f . The above relations suggest that For any given w d f , Ω d f and the parameters of the model, the LHS of eq. (35) vanishes (fixed point exists) if r 0 satisfies the above expressions.
Here, x 1 f and λ 4 are presented in terms of w d f and Ω d f , therefore the different values of w d f and Ω d f corresponds to different values of x 1 f , λ 4 as well as x 2 f (according to eq. (58)). In terms of r 0 , the fixed points in table (1) now can be presented in table (2). class I: Table 2 The fixed points for the disformal scaling solutions. For these fixed points, the parameters λ 1 and λ 3 are arbitrary, while λ 2 , λ 4 and  Fig. 1 The regions I, II and III represent the regions in which the fixed point II + is a saddle point for the cases where the values of (λ 3 , λ 1 ) are (0, 2), (0, 5) and (0, 10) respectively. The fixed point is stable outside these regions. Note that the region II totally overlaps with a part of region III, and the region I totally overlaps with a part of region II.
It can be checked that the fixed points in the class II reduce to the fixed points in [23] when λ 3 = 0. Since the analytic expressions for the above fixed points are complicate, we perform the stability analysis using numerical calculation and show the region of the cosmological parameters at fixed point (w d f , Ω d f ) in which the fixed point is stable in figures (1) -(4). In our consideration, λ 4 , x 1 f and x 2 f can be computed from w d f and Ω d f , while λ 2 and r 0 can be computed from w d f , Ω d f , λ 1 and λ 3 . Therefore, the stability of the fixed points for the disformal scaling cases can be explored by plotting the stability regions of the fixed points in the w d f -Ω d f plane for various values of λ 1 and λ 3 .
Let us first consider the case λ 3 = 0. It is clear that this case is not allowed for the fixed points in the class I. The numerical investigation shows that the fixed point II − is stable for a wide range of λ 1 when −1 < w d f < 0 and 0 < Ω d f < 1. However, it follows from figure (1) that the parameters region in which the fixed point II + is saddle increases area with increasing λ 1 . The fixed point II + can become a saddle point within the region w d f ∈ (−0.99, −0.97) and Ω d f ∈ [0.7, 1) until λ 1 ≫ 1 and λ 3 ≫ 1.
We The fixed points in the class I are saddle points for the whole region of −1 < w d f < 0 and 0 < Ω d f < 1 when λ 1 = 1 and λ 3 = 1. However, according to figure (3), these fixed points become stable within some regions in the w d f -Ω d f plane when λ 1 or λ 3 increases. For this class, the fixed points can be stable within the region w d f ∈ (−0.99, −0.97) and Ω d f ∈ [0.7, 1) if λ 1 is sufficiently large, i.e., the slope of the conformal coefficient is large.
We now consider whether the same value of parameters λ 1 , λ 2 , λ 3 and λ 4 can lead to different stable fixed point at late time. Since the fixed points in the class I do not support λ 3 = 0, we consider the case where λ 3 = 0 for the fixed points in the class II. In the case where λ 3 = 0 , eq. (64) becomes This shows that for a given value of r 0 , or equivalently increases, and shifts to the larger value when λ 3 increases. According to our numerical investigation, we also find that w s d f gets closer to 0 as λ 3 gets larger, but w s d f will not exist when λ 1 30. This means that the fixed point II + can take only one physically relevant value for a given value of the parameters when λ 1 30. From figure (4), we see that for the fixed point II − , the solution w s d f shifts towards 0 when λ 3 increases. From the detail of the numerical analysis for fixed point II − , we find that Ω d f associated with w s d f becomes larger than unity when λ 1 1. Nevertheless, the value of Ω d f can be reduce by enhancing the value of λ 3 , e.g., Ω d f < 1 for λ 1 = λ 3 = 5. For both II + and II − fixed points, the solution w s d f does not exist if Ω * d f 0.9. Based on the above analysis, we conclude that the fixed point associated with w s d f will not exist if λ 1 is sufficiently larger than unity or the values of r 0 and λ 4 correspond to Ω * d f 0.9.
We now study the situation in which the different fixed points with the same value of the parameters of the model can be reached. In order to perform, we solve eqs.  such that the initial value of w d is significantly larger than w * d f , the universe will evolve towards the fixed point associated to the solution w s d f in figure (4). Unfortunately, the present value of w d may lie outside the observational bounds if the universe evolves towards this fixed point. Hence, the existence of the solution w s d f seems to be a problem, which can be avoided by setting λ 1 to be sufficiently larger than unity or setting the value of r 0 and λ 4 to be matched with Ω * d f 0.9. We stress that this conclusions are based on the situation where λ 3 > 0 and λ 1 > 0.

Conclusions
In this work, we study dynamics of the universe at late time for the model in which dark energy directly interacts with dark matter through disformal coupling. When the disformal coefficient depends on the kinetic terms of scalar field, there exist two classes of fixed point which can describe the acceleration of the universe in addition to that found in literature. The fixed points in the first class are saddle points unless λ 1 is sufficiently larger than unity, and exist only for the case where the disformal coefficient depends on the kinetic terms of scalar field. The fixed points in the second class can be stable within the parameters ranges that correspond to the accelerating universe.
In the case where the disformal coefficient depends only on the scalar field, the fixed points in the second class become the fixed points that found in the literature. Interestingly, in the case where the disformal coefficient also depends on the kinetic terms of scalar field, the stable fixed points in the second class can take different physically relevant values for the same value of the parameters of the model. For the case where λ 1 ∼ λ 3 ∼ 1 and the value of r 0 and λ 4 are set such that the fixed point can occur at 0.9 Ω d f 0.7 and w d f ∼ −0.99, the universe will evolve towards the fixed point (w d , Ω d ) = (w d f , Ω d f ) if the initial value of w d is close to w d f . Nevertheless, if the initial value of w d is sufficiently larger than w d f , the universe will evolve towards another value of the fixed point at which the present value of w d may not be in agreement with the observational bounds. The existence of two different values of the fixed point for the same value of the parameters can be avoided if λ 1 is sufficiently larger than unity, or the values of r 0 and λ 4 are set from the fixed point w d f ∼ −0.99 and Ω d f 0.9.