Dynamical system analysis of three-form field dark energy model with baryonic matter

A cosmological model having matter field as (non) interacting dark energy (DE) and baryonic matter and minimally coupled to gravity is considered in the present work with flat FLRW space time. The DE is chosen in the form of a three-form field while radiation and dust (i.e; cold dark matter) are the baryonic part. The cosmic evolution is studied through dynamical system analysis of the autonomous system so formed from the evolution equations by suitable choice of the dimensionless variables. The stability of the non-hyperbolic critical points are examined by Center manifold theory and possible bifurcation scenarios have been examined.


Introduction
The cosmologists for the last 2 decades are tirelessly searching for a theory which shows the present accelerated expansion of the universe as predicted by a series of observational results [1]. Cosmologists are not unanimous over this issue rather they have two distinct opinions -one group is in favour of standard cosmology and they introduce some exotic matter (known as DE having large negative pressure) to explain this accelerated expansion while the other group prefers some modified gravity theory instead of Einstein gravity and they speculate that the extra geometric terms may responsible for this accelerated expansion [2].
In the context of the opinion of the first group the best candidate so far for DE is the cosmological constant which is simple in nature and observationally most favourable. But unfortunately, it has two severe drawbacks namely cosmological constant problem and cosmic coincidence problem [3]. As a result, several dynamic DE models in the form of perfect fluid with variable equation of state parameter (of various forms) come into picture (namely quintessence, a e-mail: soumyachakraborty150@gmail.com b e-mail: sudipcmiiitmath@gmail.com (corresponding author) c e-mail: schakraborty.math@gmail.com tachyon, phantom [4] and chameleon [5,6] etc). Also some complicated fields namely spinors [7], vectors [8], higher order spin field and three-form field [9,10] are considered as DE (also one may refer [11][12][13][14] on the p-form and their stability properties, the three-form inflation and the center manifold). Also other positive aspects of the three-form fields in the context of cosmology are to obtain dynamically the latetime acceleration, to describe the phantom like behaviour [15] and to produce non-Gaussianities [16].
In the present work the three-form field is chosen as DE candidate. It has been shown that the dual of the threeform field is a scalar field which has an equivalence with K -inflation model [17] if the potential is non-quadratic in form. The present cosmological model contains this threeform field with (interacting, non-interacting) baryonic matter in the form of dust and radiation. Due to complicated form of the Einstein field Equations, the evolution equations are converted into an autonomous system by suitable choice of the dimensionless variables and non-hyperbolic equilibrium points are analyzed by Center manifold theory on the Hilbert space [18][19][20]. Bifurcation analysis [21] has been done to find the qualitative changes of global behavior due to different parameter values. The plan of the paper is as follows: Sect. 2 deals with basic equations for the three form field, Einstein field equations and energy conservation equations of both types of fluids. Autonomous system is formed and critical points are determined in Sect. 3. Also stability analysis and possible bifurcation scenarios have been examined in this section. Finally, the summary of the present work is proposed in Sect. 4.

The model and basic equations
For a three-form field A μνρ , the field strength tensor is given by where square bracket indicates anti-symmetrization of the indices. The action of this three-form field is given by where V (A 2 ) is the potential of the field with A 2 = A μνρ A μνρ . The energy-momentum tensor corresponding to this action is given by and the evolution equation of the three form field is obtained by varying the above action (2) with respect to A μνρ as α F αμνρ = 12V (A 2 )A μνρ (4) where an over dash denotes differentiation with respect to the argument. The present cosmological model is considered in the background of homogeneous and isotropic flat Friedmann-Lemaître-Robertson-Walker (FLRW) space-time manifold having line element (choosing c = 1) ds 2 = −dt 2 + a 2 (t) dr 2 + r 2 dθ 2 + sin 2 θ dφ 2 . (5) So the three form field has only time dependence and the above evolution Eq. (4) results A 0μν = 0. Thus the spatial components of the three-form field can be written as where i jk is the three dimensional Levi-Civita symbol (with 123 =1). So the scalars A 2 and X are related as A 2 = 6X 2 and the equation of motion (4) of the three-form field gives the evolution of the scalar field X (t) as where differentiation with respect to the cosmic time 't' is denoted by over dot in the above equation. Now from the expression (3) of the energy-momentum tensor for the threeform field, the explicit form of the energy-density and the pressure of the scalar field X are given by having equation of state parameter Thus if V, X (= dV d X ) is positive then the three-form field behaves as DE (i.e; non-phantom) while the field behaves as phantom field if V, X < 0. Now for the present cosmological model the action integral for the whole system (assuming the three-form field to be minimally coupled to gravity) is given by where κ 2 = 8π G, R is the usual Ricci scalar and S B , the standard action integral for the baryonic matter (in the form of dust and radiation) is given by The last equality in the above expression for S B holds upto a total derivative [22]. Here ρ B and p B are the energy density and thermodynamic pressure of the baryonic matter respectively. Now varying the action with respect to the metric tensor gives the usual Friedmann equations as where ρ r , ρ d are respectively the energy density of radiation and dust part of the matter and ω r = 1 3 is the equation of state parameter for radiation. Also the energy conservation equations are [assuming three form field interacts with dust, the cold dark matter (CDM)] Now due to this interaction, the evolution of the scalar field X (i.e., Eq. (7)) modifies tö In the present work, two possible choices for Q are taken as (i) Q = αρ d (Ẋ + 3H X) and (ii) Q = αρ d H , where α is the coupling parameter. Note that the above choices of Q are purely phenomenological and the only motivation of choosing such Q is the formation of the autonomous system (in Sect. 3 below) with the choices of the dimensionless variables [in Eq. (20)]. Further, in order to derive the evolution Eqs. (15) and (17) for the three-form field and dust, the interaction Lagrangian is of the form [22] where ' f ' is an arbitrary function which will specify the particular model. Also ' f ' depends only on the dynamical degrees of freedom of the fluids.

Formation of autonomous system: critical point and stability analysis
We now define a set of dimensionless variables as [23] x ≡ κ X, y ≡ κ the equation of state parameter can be written in the form and the cosmic evolution equations (in the last section) can be written in an autonomous system as where 'dash' over a variable denotes differentiation with respect to N = ln a. Note that all the above dimensionless variables are not independent, rather due to first Friedmann equation (i.e; Eq. (13)) they are constrained by the relation and we have 4D phase space for the dynamical system. In the present work for simplicity of calculation the potential of the scalar field X (i.e; V (X )) is chosen as where V 0 (> 0) and μ are constant parameters, so that λ(X ) is a non-zero constant. Also μ can be interpreted as the rate of decrease of the logarithm of the potential function. The properties of the critical points of the above autonomous system (22)(23)(24)(25) are represented in the following lemmas.

Lemma 1 Critical points can be located either on unullcline or on v-nullcline.
Proof Let us prove this by contradiction. We assume that both v and u = 0. Then from Eqs. (24) and (25) we have 1 + R = 0 and But it is not possible to exist any such R. 2 and λ(x)y < 0.

Non-interacting three-form field
The vector fields of the autonomous system (22-25) for noninteracting three-form field for exponential potential can be analyzed as follows. By taking exponential potential we will get λ(x) = μ. The set of critical points, existence of critical points and the value of cosmological parameters are shown in Table 1 and the eigenvalues and the nature of critical points are shown in Table 2.
The first two critical points (i.e., C 0 , C 1 ) represent the evolution of the universe with barotropic matter and vanishing DE. In fact, the critical point C 0 describes the relativistic radiation era of evolution while C 1 represents the dust phase of evolution. The remaining three critical points namely C 2 −C 4 are fully dominated by DE and they correspond to de-Sitter phase.

Stability analysis
Critical point C 0 The Jacobian matrix at the critical point C 0 can be put as The eigenvalues of J (C 0 ) are −3, 2, 1 2 and 4 (hyperbolic in character) with eigenvectors [1, 0, 0, 0] T , can use Hartman-Grobman theorem for analyzing the stability of this critical point. As three eigenvalues are positive and one is negative, so the critical point C 0 is a saddle node and unstable in nature. For μ > 0 the phase portrait in all possible 3D coordinate system near the origin are shown as in Fig. 1.
The Jacobian matrix at the critical point C 1 can be put as The eigenvalues of J (C 1 ) are −3, 3 2 , 3 and − 1 2 (hyperbolic in character) with eigenvectors [1, 0, 0, 0] T , 2 we can use Hartman-Grobman theorem for analyzing the stability of this critical point. Since two eigenvalues are positive and another two eigenvalues are negative, hence the critical point C 1 is unstable due to its saddle nature.
Critical point C 2 The Jacobian matrix at the critical point C 2 can be put as (33) Fig. 1 All possible 3 dimensional phase plots corresponding to the critical point C 0 (0, 0, 0, 1) for μ > 0; a represents phase plot in xyv-coordinate system (saddle node), b represents phase plot in yvucoordinate system (unstable node), c represents phase plot in xyu-coordinate system (saddle node), d represents phase plot in xvucoordinate system (saddle node). We refer the interested readers to Appendix A for corresponding Mathematica code respectively. To apply Center Manifold Theory, we first trans- so that C 2 moves to the origin. As a result, the autonomous system (22)(23)(24)(25) changes to Now, we can see that there is a linear term of Y in the R.H.S. of (34), so we have to introduce another set of new coordinates so that the Jacobian matrix (33) transforms into the diagonal form. By using the eigenvectors of the Jacobian matrix (33), we introduce the following coordinate system ⎡ and in these new coordinates the equations are transformed into ⎡ So by center manifold theory there exists a continuously dif- Differentiating both sides with respect to N , we get where We only concern about the non-zero coefficients of the lowest power terms in Center Manifold Theory as we analyze arbitrary small neighborhood of the origin. Comparing coefficients of the corresponding power of Y T , we get the following center manifold This implies that one dimensional center manifold lies on the (X T Y T )-plane and tangent to the center subspace (Y T axis) at the origin. The flow on the center manifold near the origin is determined by The flow on the center manifold depends on the sign of μ.
For both the cases μ > 0 and μ < 0, the origin is a saddle node and unstable in nature (Fig. 2). As the new coordinate system (X T , Y T , V T , U T ) is topologically equivalent to the old one, hence the origin in the new coordinate system, i.e., the critical point C 2 in the old coordinate system (x, y, u, v) is a saddle node and unstable in nature.
The Jacobian matrix at the critical point C 3 and the corresponding eigenvalues and eigenvectors are same as for the critical point C 2 . If we put forward similar argument as we have mentioned for the analysis of C 2 then we get the following center manifold and the flow on the center manifold is determined by Since our interest only on the coefficient of the lowest power term of Y T in the expression of the center manifold and also the flow on the center manifold, so the vector field near the origin is same as for the critical point C 2 (Fig. 2). Hence, the critical point C 3 is unstable due to its saddle nature in the old coordinate system.
Critical point C 4 The critical point C 4 exists only for μ = 0. The Jacobian matrix at the critical point C 4 is same as (33), so the eigenvalues and corresponding eigenvectors are same as for (33). Since, the critical point C 4 is non-hyperbolic in nature so we can use center manifold theory. But here the center manifold is determined by and the flow on the center manifold is determined by Since, from this information we can not determine the nature of the vector field near the origin. So, we only try to define stability of the vector field near the origin on each plane. The stability of the vector field on each plane are shown in Table 3. Coordinate plane Stability Vector field near the origin is stable star Vector field near the origin is stable star Vector field near the origin is stable star 3.2 Coupling three-form field with dark matter The existence of the coupling can be represented by the modified continuity equations (15) and (17), where ρ d stands for energy density of dust, ρ A is the energy density of three-form field and Q is the energy transfer between dark energy and dark matter. Then the autonomous system changes due to Eqs. (15) and (17) as follows where I = κ Q √ 6(Ẋ +3H X)H 2 and 'prime' denotes derivative w.r.t N = ln a.
The properties of the critical points of the above autonomous system (depending on the choices of I ) are presented in the form of lemmas as follows (we shall consider the relation x = 2 3 y due to Eq. (56)).

Corollary 1
If only u = 0 and I = v y α, then α = − 3v 2 and critical points contain any real y,v and z with y 2 +v 2 + z 2 = 1.

Lemma 7
If only u and z = 0, then I = − 3 2 v 2 y and critical points contain any real v and y with v 2 + y 2 = 1.

Corollary 2
If only u, z = 0 and I = v y α, then α = − 3 2 vy 2 and critical points contain any real v and y with v 2 + y 2 = 1.

Lemma 8
If v = 0 and I = vα, then y must be 0 and u = ±1.

Lemma 9
Consider I = v 2 α. Then if v = 0, then z = 0 and critical points contain any real u and y with u 2 + y 2 = 1.
So eliminating √ 6λ(x)yz 2 from (66) and (67)   (2) If both u and v are non-zero, then y has to be non zero. Otherwise, Theorem 1 contradicts. Now by Lemma 5 z turns out to be 0 and I = v 2 2y . On the other hand, if we put the expression of I and 4u 2 + 3v 2 = 4 in (57) we derive that v 2 + 4y 2 = 0 which implies both v and y are 0. But this contradicts our hypothesis. 2 v 2 , αv 2 2y respectively. Now for these two choices of Q we analyze the stability of the vector field near the origin for the autonomous system (56-59).
For this choice of interaction the autonomous system (56-59) changes to The critical points, existence of critical points and the value of cosmological parameters are shown in Table 4 and the eigenvalues and the nature of critical points are shown in Table 5.
In Table 4 there are six critical points of which the first four namely P 0 -P 3 are equivalent to the critical pointsC 0 , C 2 -C 4 respectively in Table 1 from cosmological view point. The critical point P 4 describes the cosmic evolution due to interacting DE and CDM for α = ±3. For α = ±3; the matter is purely in the form of cosmological constant. So this critical point corresponds to scaling solution in cosmology. As long as the coupling parameter 'α' is restricted as α 2 < 3 then CDM dominates over DE otherwise there is accelerated expansion. For critical point P 5 the DE in the form of a scalar field interacting with CDM. But here cosmic evolution is dominated by DE and is in the de-Sitter era of evolution.
For all μ and f or all Always non-hyperbolic for any α Always non-hyperbolic for any α Always non-hyperbolic for any α Non-hyperbolic for α = ±3 and hyperbolic for all α (−3, 0) ∪ (0, 3)

Stability analysis
Critical point P 0 The Jacobian matrix at P 0 is same as (31). So the eigenvalues and the corresponding eigenvectors are also same. Since, the critical point is hyperbolic, we can analyze the stability of this critical point by Hartman-Grobman theorem. The stability of this critical point is same as the stability for C 0 (Fig. 1).

Critical point P 1
The Jacobian matrix at the critical point P 1 can be put as Here we will take four choices of α for analyzing the stability of this critical point : the corresponding eigenvectors. We put forward similar argument as we have mentioned for the analysis of C 2 then the center manifold is given by (44-46) and the flow on the center manifold is determined by (47). So, the origin is a saddle node and unstable in nature (Fig. 2). Hence, in this case the origin in the new coordinate system, i.e., the critical point P 1 in the old coordinate system is a saddle node and unstable in nature.

Case (ii): α = −3
In this case the Jacobian matrix at the critical point P 1 can be put as The eigenvalues of the above Jacobian matrix are −3, 0, 0, −2 (non-hyperbolic in nature). Since the algebraic multiplicity and the geometric multiplicity corresponding to each eigenvalues are equal so we can transform (73) to the diagonal form. Similarly as above we will take the same transformations so that P 1 moves to the origin and then introduce another transformation (38) (by using the eigenvectors of (73)) to make the Jacobian matrix (73) into the diagonal form. In these coordinates, our system of equations So by center manifold theory there exists two continuously differentiable functions χ :R 2 → R and φ:R 2 → R such that Now differentiating both sides with respect to N , we get Comparing L.H.S. and R.H.S. of (77) and (78), we get a 1 = − 4λ 3 , a 2 = 0, a 3 = 0 and b i = 0 for all i. Then the center manifold is given by The flow on the center manifold is determined by Then the coefficient of lowest power terms of (81) and (82) both are positive. Now divide both sides of (81) and (82) by 2 √ 6 and ( 3 2 + √ 6μ) respectively and since by dividing both sides any positive number there will no effect on the flow of the vector field, we take r 2 = Y 2 T + V 2 T , in arbitrary small neighborhood of the origin. Differentiating both sides with respect to N , we get r = Y T r . For Y T < 0 one has r < 0 while r > 0 for Y T > 0. So, for μ > 0 the origin is a saddle node and unstable in nature (vector field near the origin is same as Fig. 2b).
Then the coefficient of lowest power terms of (81) and (82) both are negative. We assume 2 √ 6μ=-σ 2 and ( 3 2 + √ 6μ) = −γ 2 and divide (81) by σ 2 and (82) by γ 2 and we take r 2 = Y 2 T + V 2 T , in arbitrary small neighborhood of the origin. Differentiating both sides with respect to N , we get r = −Y T r . For Y T < 0 one has r > 0 while r < 0 for Y T > 0.
So, for μ < − 1 2 3 2 the origin is a saddle node and unstable in nature (vector field near the origin is same as Fig. 2a).
As for both of the subcases in this case the origin is unstable due to its saddle nature. Hence, in the old coordinate system (x, y, v, u) the critical point P 1 is a saddle node, i.e., unstable in nature. Critical point P 2 The Jacobian matrix at the critical point P 2 can be put as For avoiding similar calculation, we only state the expression of the center manifold and the flow on the center manifold.
Here we also take four choices of α for analyzing the stability of this critical point : are the corresponding eigenvectors. In this case the center manifold is same as (48-50) and the flow on the center manifold is determined by (51). So, the origin is a saddle node and unstable in nature (Fig. 2). Hence in the new coordinate system due to saddle nature of the origin, the critical point P 2 is a saddle node and unstable in nature.
Case (ii): α = 3 In this case, the center manifold is given by The flow on the center manifold is determined by Case (ii)(a): μ > 1 2 3 2 In this case as above, divide both sides of (86) and (87) by 2 √ 6μ and (− 3 2 + √ 6μ) respectively and by taking r 2 = Y 2 T + V 2 T , we can get r = Y T r . If Y T > 0 then r > 0 and r < 0 while Y T < 0. So, the origin is a saddle node and unstable in nature (same as Fig. 2b).
Case (ii)(b): μ < 0 In this case also by taking r 2 = Y 2 T + V 2 T , we can get r = −Y T r and as above we can determine that the origin is a saddle node and unstable in nature (same as Fig. 2a).
As for both of the subcases in this case also the origin is a saddle node, hence in the old coordinate system (x, y, v, u) the critical point P 2 is a saddle node and unstable in nature.
In this case we will get the same eigenvalues and eigenvectors as in Case (iii) of the critical point P 1 and the center manifold is given by (48-50) and the flow on the vector field near the origin is determined by (51).

Case (iv) α = −3
In this case we will get the same eigenvalues and eigenvectors as in Case (iv) of the critical point P 1 and the center manifold is given by (48-50) and the flow on the vector field near the origin is determined by (51).

Critical point P 3
The Jacobian matrix at the critical point P 3 can be put as and the center manifold is same as (52-54) and the flow on the center manifold is determined by (55). So, here we only define the stability near the origin on each coordinate plane which is shown in Table 6.

X T Y T -plane
Vector field is stable about Y T axis X T V T -plane Vector field near the origin is stable star (if y c = 0 or α > − 3 yc ), vector field near the origin is a saddle node (if α < − 3 yc ), vector field is stable about V T axis (α = − 3 yc ) X T U T -plane Vector field near the origin is stable star yc ), vector field near the origin for α = − 3 yc is parallel to V T axis and the direction of the vector field is from negative V T axis to positive V T axis Vector field is stable about Y T axis U T V T -plane Vector field near the origin is stable star (if y c = 0 or α > − 3 yc ), vector field near the origin is a saddle node (if α < − 3 yc ), vector field is stable about V T axis (for The Jacobian matrix at the critical point P 4 can be put as Table 7 The set of critical points, existence of critical points and the value of cosmological parameters corresponding to the critical points for the autonomous system (90-93) For all μ and α 2 3 For all μ and f or all  , λ 2 , λ 3 , λ 4 ) of the Jacobian matrix for the autonomous system (90-93) corresponding to the above critical points and the nature of the critical points: Critical points λ 1 λ 2 λ 3 λ 4 Nature of critical point The critical points, existence of critical points and the value of cosmological parameters are shown in Table 7 and the eigenvalues and the nature of critical points are shown in Table 8. The first three critical points B 1 -B 3 in Table 7 are cosmologically equivalent to the critical points C 2 -C 4 in Table 1. The critical point B 4 represents a scaling cosmological solution where DE is interacting with CDM. The CDM dominates the evolution for −1 < α < 0, otherwise there is accelerated expansion due to dominance of DE. For α = −3 the matter is purely in the form of cosmological constant.

Stability analysis
The Jacobian matrix at the critical point B 1 is same as (72). So here we also take four choices of α (same as P 1 ) for analyzing the stability of this critical point.
In this case we have the same eigenvalues and corresponding eigenvectors of the Jacobian matrix. We put forward similar argument as we have mentioned for the stability analysis of P 1 in case (i), then the center manifold is given by (44)(45)(46) and the flow on the vector field near the origin is determined by (47) (Fig. 2).
In this case we also get the same Jacobian matrix (73), so have the same eigenvalues and corresponding eigenvectors. We put forward similar argument as we have mentioned for the analysis of P 1 in case (ii), then we get the center manifold same as (79-80) and the flow on the center manifold is determined by Case (ii)(a): μ > 0 Divide (94) and (95) by 2 √ 6μ and √ 6μ respectively and as the flow on the vector field is not changed due to divide by positive number on the flow equation, we take r 2 = Y 2 T + V 2 T , then we get r = Y T r . If Y T > 0 then r > 0 and r < 0 while Y T < 0. So, the origin is a saddle node and unstable in nature (Fig. 2b).
We put forward similar argument as we have mentioned above then we also obtain that the origin is a saddle node and unstable in nature (same as Fig. 2a). Table 9 Center manifold, flow on the CM and the stability near the origin corresponding to the critical point B 2

Case
Center manifold Flow on the center manifold Stability α = −3, 1, 3 Origin is a saddle node for both of the cases μ > 0 and μ < 0 6μV T Y T + higher or der ter ms Origin is a saddle node for both of the cases μ > 0 and μ < 0.
Origin is a saddle node for both of the cases λ > 0 and λ < 0

X T U T -plane
Vector field near the origin stable star , vector field for α = −3 is shown in Fig. 3

Y T U T -plane
Vector field is stable about Y T axis U T V T -plane Vector field near the origin is stable star (if α > −3), vector field near the origin is a saddle node (if α < −3 ) As for both of the subcases in this case the origin is unstable due to its saddle nature. Hence, in the old coordinate system (x, y, v, u) the critical point B 1 is a saddle node, i.e., unstable in nature. The Jacobian matrix at this critical points is also same as (72).
Here also arises four cases. The center manifold and the flow on the center manifold are shown in the Table 9. As for all possible cases in the new coordinate system (X T , Y T , V T , U T ) the origin is a saddle node, i.e., unstable in nature. Hence, in the old coordinate system (x, y, v, u) the critical point B 2 is a saddle node and unstable in nature.

Critical point B 3
The Jacobian matrix at the critical point B 3 is same as (72).
Here, we will get the center manifold same as (52-54) and the flow on the center manifold is same as (55). So, here we only define the stability of the vector field near the origin on each coordinate plane, shown in Table 10.
We clearly see that this critical point is non-hyperbolic only when α = −3 and the stability analysis for this case already studied for B 1 . So here we analyze the stability of this critical point only when α (−3, 0). For this case the two eigenvalues are positive and another two are negative. So by Hartman-Grobman theorem we can say that the critical point B 4 is unstable due to its saddle nature.

Global behavior and bifurcation analysis
For non-interacting model with exponential potential, we have five critical points (C 0 -C 4 ) for various choices of μ. It is to be noted from matrix (32), that C 0 is saddle with one stable and three unstable eigen-directions. On the other hand, C 1 is saddle with two stable and two unstable eigen-directions. The stability of C 0 and C 1 do not depend on the parameter μ. At C 2 and C 3 the effective potential [10] reaches its extreme point in the dark energy dominated period. We can have global behavior of phase space and the main property of bifurcations on different nullclines. On the Y -nullcline, the vector field does not depend on μ. But on the center manifold of C 2 and C 3 the flow reverses its direction when μ passes through the value μ = 0. But the vector fields remain topologically equivalent even after μ changes its sign. At μ = 0, the stability is classified by considering the flow of the remaining eigen-directions. At μ = 0, new non-isolated critical point C 4 appears to exist. C 4 is normally hyperbolic and stable in nature. It is to be noted that, for μ = 0, C 2 and C 3 are saddle-node in nature. So there exists a non-generic de-Sitter evolution of the universe through only one trajectory (i.e; invariant manifold) from past (C 2 or C 3 ) to future (C 3 or C 2 ). For μ > 0, if the trajectory starts near C 2 and ends asymptotically to C 3 , then the three form field transit from phantom field to non-phantom field and reverse transition happens for μ < 0. As C 4 is stable with de-Sitter phase, so there is a generic evolution near C 2 or C 3 towards C 4 at μ = 0.
Next, we have considered the interaction Q = α(ρ d )(Ẋ + 3H X). For exponential potential there exists three critical points (P 0 -P 5 ) for various choices to μ. At α = ±3, P 1 and P 2 change its stability. On the other hand, at μ = 0, new non-isolated critical point P 3 appears to exist. P 3 is normally hyperbolic and on the Y -nullcline P 3 becomes stable node to saddle node when α becomes smaller than − 3 y c , (y c = 0) for any fixed y c (0 < y c 1). When α touches ±3 and lies in [−3, 0) ∪ (0, 3], a new critical point P 4 appears to exist. P 1 and P 2 are the special case of P 4 at α = ±3. For interaction Q = α(ρ d )y H, at α = 3, −3, both P 1 and P 2 are saddle node in nature. So there exists a non-generic de-Sitter evolution of the universe through the invariant manifold (86) from past P 2 to future P 3 .
For the interaction αρ d H , a new critical point B 4 appears to exist for α ∈ [−3, 0). This critical point changes its stability at −3 and μ = 0. For μ = 0, B 4 is saddle node in nature. So the evolution of the universe through the invariant manifold is performed as in the previous case. Now, it is to be noted that at μ = 0, the system associated with non-interacting model becomes structurally unstable. On the other hand, two new bifurcation values appear for the interacting model, namely α = −3, 3. Moreover, at μ = 0, the system associated with interacting model becomes structurally unstable. Also one may note that the potential becomes runaway to non-runaway through bifurcation value μ = 0 [23]. On the other hand, if we neglect the contribution of matter and radiation, the fixed points at the extremum point of the effective potential [10] are same as in our non-interacting dynamical equation i.e., C 2 and C 3 which trigger non-generic evolution of the universe.

Summary
The present work is an example where without making any attempt of solving the complicated coupled cosmic evolution equations (Einstein field equations), the cosmological interferences have been done using the technique of dynamical system analysis. The cosmic matter is chosen as three-form field (DE) interacting/non-interacting with baryonic matter (radiation with CDM). By suitable choice of the dimensionless variables the evolution equations are converted into an autonomous system for non-interacting case and for two suitable choices of the interaction. Out of total fifteen equilibrium points (presented in Tables 1, 4 and 7) the following sets of critical points (C 0 , P 0 ), (C 2 , P 1 , B 1 ), (C 3 , P 2 , B 2 ), (C 4 , P 3 , B 3 ) are equivalent both cosmologically as well as from dynamical system analysis. It is to be noted that the set of critical points (C 4 , P 3 .B 3 ) are essentially a line of critical points. The first set describes the radiation era of evolution while the remaining three sets correspond to accelerated expansion at the de-Sitter phase. The critical point C 1 describes the CDM dominated decelerated expansion. Although the critical point P 5 represents interacting DE and CDM but it also corresponds to de-Sitter phase. The most interesting critical points are P 4 and B 4 which are representing scaling cosmological solutions. Also the cosmic evolution representing these two critical points are observationally important as by choosing the coupling parameter α to be closed to ±3 (for critical point P 4 ) / −3 (for B 4 ) the equation of state parameter agrees (within confidence limit) with the latest plank data [24].
The stability of the critical points are analysed by studying the eigenvalues of the Jacobian matrix for hyperbolic critical points while center manifold theory has been employed for non-hyperbolic critical points. The two parameters μ (in the potential) and α (coupling) on the system are found to be important for bifurcation analysis. Usually, it is speculated that at the bifurcation point there may be a cosmic phase transition.
Finally, the general nature of the critical points both for non-interacting and interacting cases have been presented in the form of lemmas and theorems. Lastly, one may note that the present work does not depend effectively on the choice of V (X ). Fig. 3 Vector field near the origin in Y T V T -plane for α = −3. After shifting and matrix transformation, the line of critical points B 3 (in the old coordinate system (x, y, v, u)) is completely lying on the Y T axis in the new coordinate system. For determining this vector field we have taken y c = 1 2 . So, in the new coordinate system the origin in Y T V Tplane is corresponding to 0, 1 2 , 0, 0 in old coordinate system. We refer the interested readers to Appendix B for corresponding Matlab code