Phase space of multi-fluid universe in $F(T)$-gravity and some enhancements for the oscillating interaction model

Recently, a Friedmann-Robertson-Walker universe filled with various cosmological fluids has been considered by S.D. Odintsov et al. in [30] from phase space vantage point where various expressions for the Equation-of-State (EoS) parameter were studied. Since these types of EoS parameters are generative of appreciable results in the Hilbert-Einstein model, hence we intend to investigate all the cases in a homogeneous $F(T)$-gravity ($T$ is the torsion) through phase space analysis in precise detail. In short, three viable models of interaction between dark matter and dark energy, including usual-type, power-law type, and oscillating type, are investigated comprehensively. It is indicated that the power-law interaction in the related dynamical systems should be of increasing nature with time to get more critical points. Due to the failure of the oscillating model of ref. [30] in $F(T)$-gravity, four modified models are suggested and examined in both $F(T)$ and Hilbert-Einstein models. As to be seen, the modified models not only are generative of critical points equivalent to that of ref. [30], but also give rise to further critical points covering crucial stages of the evolution of the universe. In the context of these four models, such as the old one, at early times the interactions are negligible and they commence to grow as the cosmic time approaches the late-time in which the unification of early inflation and late acceleration is obtained. Using an indirect method, it is shown that the oscillating models have substantial roles in transitions between eras.


Introduction
According to several observational data including supernova type Ia [1,2], weak lensing [3], CMB studies [4], largescale structure [5], and baryon acoustic oscillations [6], the current status of the expansion of the universe is on the accelerating phase. The explanation of this acceleration is one of the major challenges for cosmologists because it is inconsistent with the standard Einstein's general relativity. In order to explain the late-time-accelerated expansion two main classes of approaches have been suggested up to now. In the first approach, the universe is assumed to be made up of an exotic liquid, the so-called 'dark energy' which comprises 70% of the universe. Einstein's cosmological constant was historically the first model of the accelerated expansion of the Universe. To avoid its 'fine-tuning' and 'cosmic coincidence' problems [7] and unclear physical nature the alternative models were proposed [8], called as 'dark energy' models: the tachyon field [9], quintessence [10,11], quintom [12,13], phantom field [14,15], etc. The second approach is to modify Einstein's general relativity by making the action of the theory dependent upon a function of the curvature scalar R . As one expects, the theory reduces to general relativity in certain limits of the parameters. Besides f (R) -gravity, various theories such as scalar-tensor theories, F (T ) -gravity (T is the torsion), F (T ) -gravity with an unusual term [16], and etcetera have recently been proposed with fruitful results.
As is known, there are several tools for studying cosmological models amongst which the Noether symmetry approach, reconstruction methods, and phase space analysis are mentionable. In the Noether symmetry approach, one finds the exact solutions of the model, which carry some conserved currents. Besides conserved currents, presenting the form of unknown functions is the most important feature of this tool. Indeed, the Noether symmetry approach is a powerful standard way for acquiring the forms of unknown functions of alternative theories of gravity. It is important to mention that some hidden conserved currents may not be obtained by this approach [17,18]. Carrying more conserved currents as well as new forms of unknown functions are the main motivations for enhancing the method to the 'Beyond Noether Symmetry approach' (B.N.S.-approach) [19] and 'Combination of Sub-symmetries through Special Selections' (CSSS-approach) [20] that have recently been proposed with noteworthy outcomes. There are several ways to reconstruct a model of alternative theory such as reconstruction by parameterizing the Hubble parameter in terms of the redshift [21], reconstruction using the scalar field [22,23], reconstruction using the Raychaudhuri equation [24,25], and etcetera. Some new collections of the forms of unknown functions may also be achieved through the reconstruction method especially when one does the reconstruction through the scalar field. A general prescription to reconstruction based on the scalar field has been suggested in ref. [23] using the approach proposed in ref. [22]. A general powerful approach, the 'Boundary function method' ( B-function method), to analyze the results of this approach has recently been suggested [26]. The B-function method can be applied to phase space analysis and exact solutions of alternative theories of gravity. It also can be regarded as a new reconstruction method. The phase space analysis is also a powerful tool for investigating a model. In this approach, the equations of our system are converted into first-order differential equations through defining some dimensionless variables. Although this process could be challenging in some cases, several approaches have been proposed for those challenges. In phase space analysis, we find Critical Points (CPs) of a model along with its eigenvalues and possibly some further helpful parameters. Indeed, one can use these CPs to analyze the model. In general, the nature of a CP can be attractor, saddle, repeller, and null according to its eigenvalues. Needless to say, this approach has attracted cosmologists to take it as a prospective in analyzing gravitational theories.
The phase space of a Friedmann-Robertson-Walker universe which is filled with various cosmological fluids (with and without interactions) has recently been considered in detail in ref. [30]. Several forms of the equation of state parameter, especially arisen from viscous fluid and a new class of interaction between dark energy and dark matter namely cosine type, have been studied in ref. [30] successfully. In this paper, we examine all the EoS suggested in [30] and some further improved interaction types in an F (T ) -model in detail via phase space analysis.
The paper is organized as follows: The model is presented in sect. 2. In sect. 3, some well-known features of the dynamical system approach in cosmological context with generalized fluids suggested in [30], are considered. In sect. 4, we add dark matter, which interacts with dark energy candidate, to our assumed components. The generalized formalism of sect. 4 is carefully considered in sect. 5 in two main modes. Section 6 is devoted to studying some interacting models between dark matter and dark energy. Furthermore, some enhancements to cosine-type of interaction are suggested and examined successfully. Finally, the conclusions are provided at the end of the paper.
where e = det(e i ν ) = √ −g with e i ν being a vierbein (tetrad) basis, κ 2 = 8πG, T is the scalar torsion, and L M is the matter Lagrangian.
Adopting flat FRW metric with the signature (+, −, −, −) , the Friedmann equations for (1) would be 3 where ρ and p are the total density and pressure, respectively, H is the Hubble parameter, the dot denotes differentiation with respect to the cosmic time, F ,T = dF/dT , and F ,T T = d 2 F/dT 2 . These equations describe the dynamical evolution of the model. The total energy conservation equation for the flat FRW universe, namelẏ is taken to be true for different groups of components that are studied in this paper. However, F (T ) is an unknown function of the scalar torsion, but in order to avoid computation problems, let us proceed with a well-known special form as [27] F (T ) = c 1 + 2c 2 √ −T + αT , where c 1 , c 2 , and α are arbitrary constants. Throughout this paper, we assume that α = −1 , for the left-hand-side of Einstein's equations cannot be produced when one takes α = −1 . In F (T ) -gravity, the constant and the non-linear terms have correspondence with the cosmological constant equation of state [28]. Furthermore, as it has been discussed in ref. [29], the power-law form of F (T ) (i.e. F (T ) ∼ T n (n > 1) ) such as T 2 may remove the finite-time future singularity. When n = 0 , F (T ) may be regarded as a cosmological constant. The case n = 1/2 is helpful in understanding power-law inflation and also in elucidating the little-rip and pseudo-rip cosmology [29], hence making our form more appealing than other options.

Standard approach on dynamical systems and cosmological dynamics
In this section, a simple case in which the universe is filled with radiation, ρ r , and a perfect fluid with a non-trivial EoS of the form 1 p v = −ρ v + f (ρ v ) + G(H) [32], is considered. These types of non-trivial EoS parameters may be deemed as some sort of viscous fluid or generalized EoS fluid [33] as well as an effective fluid presentation of some modified gravity theory [34,35]. The motivation of the appearance of the term G(H) which might seem unconventional, has clearly been described in ref. [30].
The corresponding continuity equations are as follows: where the prime denotes differentiation with respect to the e-folding number N = ln(a) in which a is the scale factor of the FRW universe. In order to rewrite the cosmological equations in terms of dimensionless variables, we define dimensionless density parameters as Therefore, dynamical equations turn out to be Like ref [30], a simple case of EoS is considered in this section. To this end, we assume f = ρ v (1 + w 0 ) and G = w 1 H 2 . Therefore, equations (8)-(9) take the following forms: Equating the left-hand sides of (10)- (11) to zero, three critical points are found: CP2: CP3: CP1: CP1 is in the physical region when α > −1 . It corresponds to a universe filled with a viscous fluid only (i.e. without radiation).
CP2: If one sets (w 1 /(1 + w 0 )) ≤ 0 , then CP2 will be in the physical era. There is a singularity at w 0 = −1 besides w 1 = ∞ = w 0 . Apart from the special case w 1 = 0 , that refers to the vacuum case, for w 1 = 0 this point such as CP1 shows a universe filled with a viscous fluid without radiation.
CP3: If the condition (w 1 /(3w 0 − 1)) ≤ 0 is adopted, then according to the relation y 3 = −x 3 + 1 + α , it is sufficient to take α ≥ −(1 − x 3 ) for putting CP3 in the physical region. The sum of the densities of radiation and viscous fluid is controlled by the value of α . There is a singularity at w 0 = 1/3 . The special cases w 1 = 0 and {x 3 = 1 + α & α = −1} correspond to pure radiation and pure viscous cases, respectively. And when {w 1 = 0 & w 0 = 1/3} , the background of the model will be the vacuum.
The eigenvalues of the Jacobian matrix for these CPs are as follows: CP2: CP3: where γ 0 = κ 2 w 1 /(1 + α) . CP1: Because λ 2−cp1 = λ 1−cp1 + 4 , hence CP1 can be stable or unstable depending upon the values of constant parameters: CP2: CP2 is an attractor along y -direction, while at the next direction it can be stable or unstable based on the chosen values of parameters. Pursuant to the condition of staying CP2 in the physical region, w 0 and w 1 must belong to one of the following sets: In the case that w 0 and w 1 belong to 'Set-1', if κ 2 |w1| 1+α < 3(1 + w 0 ) is satisfied, then CP2 will be an attractor otherwise it will be a saddle point. In the case that w 0 and w 1 belong to 'Set-2', if κ 2 |w1| 1+α > 3|(1 + w 0 )| is satisfied, then CP2 will be an attractor, otherwise, it will be a saddle point.
CP3: In general, CP3 is called unstable. More precisely, along x -direction, it has stability while along y -direction its nature depends upon the amounts of constant parameters. According to the physical conditions of CP3, besides 0 ≤ x 3 ≤ (1 + α) , w 0 and w 1 must belong to one of the following sets:  (11). Plot 'a' indicates the phase diagram when the values of the constant parameters are selected as: κ 2 = 1 , w0 = +1 , w1 = −1 , and α = +1 (i.e. the chosen values belong to 'Set-1' and 'Set-3'), while plot 'b' shows the phase plane when we pick up the values of constant parameters from 'Set-2' and 'Set-4' as w0 = −8 , w1 = +8 , α = +1 , and κ 2 = 1 . The pink points refer to the critical points.
If w 0 and w 1 belong to 'Set-3', then CP3 will be a saddle point, while if they belong to 'Set-4', it will be a completely unstable point (i.e. all trajectories will be repelled).
As an example, two phase portraits have been illustrated in fig. 1. The above-mentioned explanations are clear in this figure. 4 Models with dark matter interacting to dark energy In this section, beside the radiation, ρ r , and a viscous fluid, ρ v , the dark matter, ρ dm , is also added so that it interacts with dark energy via a usual 2 interaction model. The dynamical system for this case may be taken as where Q is the coupling parameter that describes the interaction and −3Hζ corresponds to the bulk viscous pressure of the dark-matter fluid. Let us assume that the bulk viscous coefficient is of the form [36][37][38] The relation between ζ and total density is estimated astronomically; see ref. [39]. However, the relation which we use here is ζ ∼ √ ρ tot. instead of the general relation ζ ∼ ρ λ tot. , which is due to the fact that this ansatz yields a good fit between astronomical constraints and the fluid approach [30]. The dimensionless variables for this case are introduced as follows: Therefore, by taking the previous EoS, f (ρ v ) = ρ v (1 + w 0 ) and G(H) = w 1 H 2 , (18)- (21) in terms of the independent dynamical variables turn out to be Although this system is analytically solvable, the results are unpresentable because they are very long excluding one CP. Thanks to the abundance of such solutions, let us use the acronym 'VLT' for such Very Long Terms, for the sake of brevity. It means that the related term is very long so that it is somewhat unpresentable, or even if it is, it does not lead to any clear outcome. For the above system, we may present only one CP as: Discussing this CP is not so interesting as the parameters cause confusion. Indeed, all the related eigenvalues are 'VLT'. About other CPs, we can say that all of them have a common property: z| all other CPs = 0 which means that for the remained CPs, the density of radiation is zero and they are on the interaction plane, namely x − y plane. In short, they are not on the radiation-dominated era. Therefore, the coordinates of these remained CPs can be used to reach some fixed relations between the densities of dark matter and dark energy.
To avoid 'VLT'-type solutions, it seems that there are two approaches: 1. Assuming some relations among constant parameters; 2. Specifying the values of all or some constant parameters. Examples clarifying these approaches will be presented in the next sections.

Generalized form of cosmological fluid
In order to include a generalized form of the equation-of-state, the case studied in Sect. 3 is extended in this section. A universe filled with radiation and a perfect fluid with a non-trivial EoS, [32], is considered here. To convert to an autonomous system of first-order differential equations, the set of dimensionless variables are defined as from which we obtain According to the definitions, one has ρ v = 3x/(κ 4 E) and H 2 = 1/(κ 2 E) . Obviously, there are only two independent dimensionless variables, but introducing a further variable, E , leads to obtaining more information about the system. This trick was firstly proposed in ref. [40]. It would be also helpful for interpreting the evolution of the system that we introduce the effective EoS. This parameter which is defined as 3 takes the following form for this case of study:

A simplified form of EoS
Let us first consider the simplest case of EoS namely Therefore, the equations (28)-(30) would be The EoS (32) also turns out to be For the above system, we obtain three critical points: CP2: CP3: where γ 1 = κ 2 w 1 /(3w 0 − 1) and m 3 is an arbitrary value. The eigenvalues of the Jacobian matrix for these CPs are found as follows: where γ 2 = κ 2 w 1 /(1 + α) . The corresponding values of the EoS parameter for these CPs take the following forms: CP1: CP1 is in the physical region if 0 < −γ 1 < 1 + α . Generally, it is an unstable point because λ 1 = λ 2 = +4 -if (γ 2 + 3w 0 − 1) < 0 , then it is completely unstable and it means that all trajectories of phase-portrait are repelled away from it, while if (γ 2 + 3w 0 − 1) > 0 , then it will be a saddle point (i.e. the trajectories of y -direction are attracted while the trajectories of x and E directions are repelled; indeed, the eigenvalues of x and E are +4 ). The positive eigenvalue of E -direction indicates the expansion of the universe and decreasing nature of the Hubble parameter. Under every condition, this point belongs to the radiation-dominated era (see (43)) and for this reason, obviously, the saddle option should be ascribed to this point in the numerical/qualitative study of phase portrait. In the special case namely (γ 2 + 3w 0 − 1) = 0 , however, the point is still unstable but let us study it deeply as it should be a saddle point. Under the condition (γ 2 + 3w 0 − 1) = 0 , CP1 has a one-dimensional center manifold. First of all, it is useful to consider the system at CP1-origin. It means that we perform a transformation as (x, y, E)| CP1 −→ (X, Y, R) = (0, 0, 0) . In order to convert the system to a standard form, three new variables (u, v, s) that are connected with (X, Y, R) via are introduced. The center manifold is represented in the form where we assume that the mappings h 1 and h 2 are of the forms where k i s are constant parameters, and they satisfy the well-known relations of the center manifold theorem. Computations culminate in the following dynamical system which is topologically equivalent to the previous system It is clearly observed that the point is completely unstable. Therefore, under the condition (γ 2 + 3w 0 − 1) = 0 , CP1 is not a good option for the period dominated by radiation. It is concluded that CP1 is a suitable candidate for the radiation-dominated era only under the following conditions: CP2: According to (38), CP2 is called a physical point if it belongs to one of the following sets: Indeed CP2 is not a critical point, but it is a Critical Line (CL) because E = an arbitrary value ≡ m 3 . This CL is so interesting as according to its effective EoS, w eff. = −1 , it demonstrates the current status of the universe with acceptable accuracy. Specifically speaking, E = an arbitrary value ≡ m 3 can be fixed exactly to the current time by tuning m 3 = (κ 2 H 2 0 ) −1 in which H 0 is the present value of the Hubble parameter. Indeed, the current status of the universe is a point of this CL. Therefore, if one desires to suppose this CP as the current time, then all of its eigenvalues should be negative 4 . According to (41), even by putting (γ 2 + 3w 0 + 3) > 0 , we still cannot distinguish its stability status and it needs to be studied by the center manifold theorem. If one assumes (γ 2 + 3w 0 + 3) > 0 , then the related center manifold is one-dimensional. We first consider this case: Like the previous case, first, the system is translated to the origin, (x, y, E)| CP2 −→ (X, Y, R) = (0, 0, 0) , and then utilizing the following relations, we get the new coordinate (u, v, s) : Using two mappings, where k i s are constant parameters, we finally get the following topologically equivalent system: Since the condition (γ 2 + 3w 0 + 3) > 0 has been assumed at first, this system indicates an attractor nature for CP2, as desired. We had this anticipation that the evolution on the center manifold will reduce to s = 0 + O(s 5 ) (i.e. s is approximately constant on the center manifold) because the equation E = 0 at CP2 is satisfied by every arbitrary value of E . In passing, we found that the current status of the universe can be achieved by CP2. The special case (γ 2 + 3w 0 + 3) = 0 , for which the related center manifold is two dimensional, implies that the coordinates of CP2 be (x, y, E) = (1 + α, 0, m 3 ) which obviously has a strong relation with CP3. Indeed, under the aforementioned condition, one has 'CP2 = CP3', provided m 3 = 0 . Let us now consider this special status of CP2. For simplicity and for being in conformity with CP3, let us take m 3 = 0 . Transforming the system to the origin, i.e. (x, y, E)| CP2 −→ (X, Y, R) = (0, 0, 0) , and going to a new coordinate (u, v, s) using the system is cast in the standard form. Utilizing a two-dimensional mapping as where k i s are constant parameters, one may reach a topologically equivalent system having the following form in the polar coordinate, (s, v) = (r cos(θ) , r sin(θ)) , In short, only in one direction we have attraction and our system is totally unstable, stemming from the utter instability of its center manifold. Therefore, under the condition (γ 2 + 3w 0 + 3) = 0 , CP2 is a saddle point. Pursuant to the justification in footnote [4], this behavior may also be acceptable for the present status of the universe. CP3: According to the effective EoS parameter of CP3, (45), four options are of interest for investigating: 1. γ 2 + 3w 0 = 1 : This condition leads to the radiation-dominated era (w eff. = 1/3 ) and it is a special form of CP1 -when CP1 has one-dimensional center manifold, that is For having this CP in the physical region, it is sufficient to set α > −1 , but, as corroborated already, in this situation the related CP is completely unstable and repeller. Therefore, a saddle point expectation makes this option not an apt candidate for the radiation dominated era.
2. γ 2 + 3w 0 = −3 : Under this condition, the special position of CP2 is recovered, viz., Thus, in this case, CP3 would be a saddle point -just in one direction we have attraction and the twodimensional center manifold is completely unstable at the CP3-origin. To wrap the argument, according to w eff. = −1 and the aforementioned discussion, one may accept this point as a candidate for the present time.
3. γ 2 + 3w 0 = −3.09 : This condition is an interesting option as it yields 5 w eff = −1.03 with three negative eigenvalues indicating a complete attraction: However, the condition γ 2 +3w 0 = −3.09 makes CP3 a good candidate for the present time but this condition causes that CP1 (with w eff | CP1 = 1/3 ) be a repeller point that is not acceptable since it must be a saddle CP.
4. γ 2 + 3w 0 = 0 : This condition corresponds to the matter-dominated era (w eff. = 0 ). The related eigenvalues turn out to be Therefore, CP3 is a saddle critical point that is completely acceptable behavior for a CP of the matterdominated era. Now, let us take other CPs into account using the aforementioned condition: Utilizing the condition γ 2 + 3w 0 = 0 , a special form of CP2 (when it has a one-dimensional center manifold) is recovered. Inserting γ 2 + 3w 0 = 0 into (61), the related system can be cast in the following form which is topologically equivalent to the main system: where s = 0 + O(s 5 ) is the center manifold which is constant (Its meaning and reason have been explained before). Therefore, CP2, under the aforementioned condition, is an attractor and is suitable for the current position of the universe. However, under the condition γ 2 + 3w 0 = 0 , CP3 and CP2 are very good candidates for the matter-dominated era and the current time, respectively, but this condition leads to a peculiarity: Due to the positivity of eigenvalues of CP1, the CP will be of repelling nature, however, we were anticipating a saddle point because of the radiation dominance.
According to all discussions and detailed analysis performed above, now it is clear that there are many options for matching with observational data and plotting. But, among them, there are only some cases with no anomalies. It is not hard to show that we cannot tune all CPs in a way that all the last three main stages of the universe be recovered, 6 for it needs to hold the following conditions simultaneously which is impossible: In the best situation, one may be able to keep two of the three aforementioned stages. Among some presentable cases, we single out the constant parameters as which lead to sufficiently good cases: where m 3 is an arbitrary constant. Three Poincaré sections x − y , x − E , and y − E have been presented in fig. 2. The physical region in these phase portraits are the positive parts: where (x 1 , x 2 , x 3 ) = (x, y, E) . Thanks to the critical line E = m 3 , three critical points can be observed in the phase portrait of x − y plane. This critical line has clearly been illustrated in x − E Poincaré section together with CP3. As is clear from this phase portrait, this critical line is an attractor -a hole line. The present position of the universe is on this line (i.e. the point (κ 2 H 2 0 ) −1 ). The trajectory of the evolution of the universe tends to this point. In y − E phase portrait none of CPs are observable. Note that the vector fields in fig. 2(c) may give us the incorrect idea that there are two CPs, (0, 0) and (7/2, 0) , of the system under consideration. In fact, this fallacy follows from just a visual error. In fig. (3), it has been tried to illustrate a historical picture based on the studied model. We have presented it in a two-dimensional plane to present a better view. In what follows, we explain its form in a three-dimensional space. The studied model (i.e. "A simplified form of EoS") with our selections to the values of constant parameters, indicates some parts of the universe evolution: Let us assume that the universe arrived at the unstable CP3 (or at least, a neighborhood of this CP) and then got repelled from 7 . This repulsion is somewhat tenable, for no stage of the evolution of the universe includes a fluid with such EoS dominance. During this repulsion and before getting to the next CP (i.e. CP1), as is clear in fig. 3, the amount of y has increased while x has decreased. Because y is proportional to the density of the radiation component, hence it makes the universe ready to enter to the radiation-dominated era. After this process, the universe reaches the radiation-dominated era at CP1 with w eff. = 1/3 . At this CP, the density of the radiation component is greater than the density of the viscous fluid that confirms the domination of the radiation component. This CP is of a saddle nature, hence, beside attraction property, it has a repulsion feature that is 6 We cannot tune the achieved CPs so that they form the following picture: • A saddle point with w eff. = 1/3 ; • A saddle point with w eff. = 0 ; • An attractor or a saddle point with w eff. = −1 .
7 Note that, without loss of generality, we can set the initial conditions so that the system could start from CP3. Hence, its nature (i.e. completely unstable feature) does not cause an anomaly.  needed to get out of this area. Overtime, the effect of this repulsion becomes obvious and the universe not only exits this CP but also the density of the radiation components decreases, then our universe is attracted to CP2 with w eff. = −1 (i.e. the present time). Note that CP2 is on the critical line whose one of its points is presented in this figure. Indeed, the final CP that has been shown in this figure, is a projection of the point being the height Since all the eigenvalues related to this point are negative, hence it is an attractor and the universe will be stable at this CP. The density of the viscous fluid is greater than the density of the radiation component, indicating the domination of dark energy at this CP. Surely, if this diagram for the universe evolution be acceptable, then between CP1 and CP2 our universe must pass a CP with w eff. = 0 . Sadly, our system yields no such CP. We pose "Why don't we have a CP representing the matter-dominated era?" Because (76)-(77) do not hold simultaneously, so reaching a perfect picture of the last three stages of the evolution of the universe is impossible. Assuming the F (T ) -gravity and the interaction models are correct theories, then two propositions may be examined to see if they could resolve this problem or not: 1-Adding baryonic matter component that has ρ m = 0 and P m = 0 . 2-Enhancing the action through a scalar field part. In this paper, we are content with the first suggestion merely. The second one is deferred to another paper. In what follows, it is shown that neither with E nor without E , the first suggestion does not work. Both only lead to w eff. = −1 and w eff. = 0 and one further EoS which cannot be fixed to a radiation-dominated era: the following dynamical system is obtained: with the following EoS The fixed points with their associated eigenvalues and EoS parameters are as follows: CP1: CP2: CP3: where χ 1 = 3w 0 + 3αw 0 + w 1 κ 2 . As observed, CP1 and CP3 may be fixed to the matter-dominated era and current stage of the universe, respectively. If one tries to set CP2 as the radiation-dominated era, the condition must be imposed. This condition makes unacceptable nature (repeller) for CP2 while other CPs are appealing candidates for the corresponding eras: It is worthwhile to mention that the radiation-dominated era cannot be reached even approximately. In detail, let us exert the following condition which contains some fluctuations through : where is an arbitrary constant. It leads to It is clear that there is no real value for leading to true behavior (i.e. saddle) for the radiation-dominated era. Indeed, if we want to keep the matter-dominated era, then we lose the radiation-dominated era and the vice-versa. Both should be saddle simultaneously which is unattainable through this system.
• With E : For the generalized case, we define These yield the following dynamical system: This system leads to the following CPs with their associated EoS and eigenvalues: CP1: CP2: CP3: where χ 1 = 3w 0 + 3αw 0 + w 1 κ 2 and m 3 is an arbitrary constant. It is clear that both cases (i.e. 'with E ' and 'without E ') are equivalent. Hence, its explanations and justifications for the inefficiency of adding a pressureless matter part are the same as the previous case. Therefore, we conclude that adding further dimension, E , does not present anything special, but two properties: By comparison with the previous case, it is easily understood that the added dimension, E , behaves as a repeller direction in general -the eigenvalue of this direction is +1 in all CPs of this subsection. This special repulsion refers to the expansion of the universe. Furthermore, the positive eigenvalues indicate that the Hubble parameter decreases with time.
In most of this paper, we encounter E = (κ 2 H 2 ) −1 = 0 for CPs. Most of these CPs belong to the radiation, matter, and current stages of the universe. Hence, one may argue that the corresponding CPs are not good candidates for the associated eras. But, note that in most cases of interest, this interpretation is wrong. This boils down to the fact that this dependent variable is just an imposition. When our dynamical system and EoS parameter are of the form 8 where k i s are constant parameters and F i s are functions of x , y , and z , all CPs with E = 0 may be discussed without serious attention to E . Note that F 1 (x, y, z) is common among them. Furthermore, if we are to keep E = 0 (which yields F 1 = 0 and w eff. = −1 ), then for reaching CPs, we must have F 2 (x, y, z) = F 3 (x, y, z) = F 4 (x, y, z) = 0 which does not tend to any CP in most cases of interest. Hence, in most cases, we have F 1 (x, y, z) = 0 and it causes E = 0 . If this additional coordinate E were not defined, then CPs would not be in the position they are now. Therefore, we can ignore the value of coordinate E in interpreting provided existing above conditions. The benefit of adding such extra dimension has been explained above: The positivity of the eigenvalues of this extra dimension not only refers to the expanding nature of the evolution of the universe but also indicates the decreasing nature of the Hubble parameter with time (E ∝ H −2 ). Due to the fact that E does not explicitly appears in (107)-(109), eq. (110) can be solved analytically at a given CP. For example, for a CP i one has: or put differently where E 0 and H 1 are the constants of integration. The positive eigenvalue of E is equivalent to the positive value of C * . The system cannot reach the zero value of E because at a critical point the corresponding to E eigenvalue is always positive. So H goes over time from infinite value to zero or some constant value. Note the form of EoS parameter and compare it with the equation E , then you will find that E = const. can occur only for the de-Sitter expansion. However, in order to comply with other CPs, one may also take E = 0 for such CPs as well, but since this constant can be fixed to the current stage of the universe, hence we would like to keep it non-zero. We know that the effective value of EoS parameter is not exactly −1 , but it is close to −1 , so, for example, if, according to Plank data [41], we have w eff. = −1.03 at the CP of the present time, then again we obtain E = 0 .
Besides the above discussions, it can be predicted for other types of systems like that if the power of E in F 1 is positive, then all CPs will have E = 0 excluding those points which belong to the de sitter era (i.e. w eff. = −1 ). In the cases in which the power of E is negative (e.g. 1/E |n| ), then there will be only one CP that belongs to the de sitter era (i.e. w eff. = −1 ) with E = const.
In the following, sections we encounter such systems.

More complicated forms of the EoS parameter
In this sub-section, some complicated forms of the EoS parameter are studied. What we are going to do is that in these forms of EoS the powers of ρ v and H in f and G , are increased and decreased. Generally, they may be cast in the It has been demonstrated that this EoS yields finite-time singularities of a cosmological system [32]. Using f (ρ v ) + G(H) = Aρ α1 v + BH 2β we arrive at: Note that α = α 1 ; α is the coefficient of the torsion in the model while α 1 is the power of the density.
The corresponding EoS parameter would also be Due to the presence of many unknown constant parameters especially the powers α 1 and β , this case may not be considered analytically. If both α 1 and β are given, then analyzing this system is feasible enough. Hence, let us consider this system by choosing some values for α 1 and β . Furthermore, two cases seem to be notable in scrutinizing all the cases of this sub-section: B = 0 and B = 0 .

Case-1:
As the first case, we take α 1 = 1 and β = 2 . With these choices the cosmological system reads: Under the above selections, the effective equation of state would then be: As observed, unlike the previous case ("A simplified case"), the amount of EoS here is affected by all three coordinates x , y , and E . ♣ Case B = 0 : In general, the system presents a two-dimensional critical curve (i.e. xE = C where C = −B/(3A) ) for a given value of B/A: Thus, for a fixed time of the evolution of the universe, i.e. for a given H , the system has only one CP because E = (κ 2 H 2 ) −1 . The conditions (AB) < 0 and A = 0 are required to put this critical curve in the physical region. On this curve, the density of the radiation component is zero while the density of the viscus fluid depends upon the values of A , B , and H . Correspondingly, the eigenvalues for the system become: It is clear that only one of the eigenvalues is affected by the stage of the universe. In the limit m 3 → ∞ (or equivalently H → 0 ), one has λ 1 = −3A. So, if we take A > 0 and B < 0 , then it means that our universe will be attracted to this point in the future. Due to λ 3 = 0 , the related center manifold is at least one-dimensional. In the case that 3Am 3 + B = 0 , the corresponding center manifold would be two-dimensional. The value of the EoS for the critical curve is which corresponds to de sitter era and the present time.
Assuming λ 1 = 0 , the system has a one-dimensional center manifold. Regardless of whether λ 1 is negative or not, we study it and finally discuss the stability and these problems. Again, let us consider the problem in CP-origin (i.e. the system is transformed into a coordinate (X, Y, R) so that the position of the CPs be (0, 0, 0) in it). Then using 10 the system is presented in a new coordinate (u, v, s) . After a tedious computational process, we get the following system which is topologically equivalent to the system at the CP-origin: It is observed that the center manifold s is unstable, hence on each point of this critical curve, this system is totally unstable and the critical curve is a saddle curve. Based on the value of λ 1 , the number of the direction of attraction and repulsion will change; if λ 1 > 0 , then in two directions we have repulsion, while for λ 1 < 0 just in one direction (on the center manifold) we have repulsion. Tuning 11 m 3 = (κ 2 H 2 0 ) −1 and λ 1 < 0 , the CP may be a good candidate for the current status of the universe as our universe is unstable. Note that, however, both λ 1 < 0 and λ 1 > 0 are acceptable, but setting λ 1 < 0 speeds up the attraction to this point (indicating acceleration) and it slows down the exit from this era. As an example, three different Poincaré sections of the phase-portrait of the above system have been presented in fig. 4 which illustrates our discussions for a negative value of λ 1 .
♣ Case B = 0 : In this case, indeed, the viscous part is removed. The critical points and their corresponding eigenvalues and EoS parameters are as follows: where m 3 is an arbitrary value. CP1: CP1 is not a critical point, but it is a critical line and its EoS parameter is −1 . Hence, at first sight, it seems better to set m 3 = (κ 2 H 2 0 ) −1 , and try to make this CP as an attractor/a saddle. But, according to the zero densities of viscous fluid and radiation, it is not a good candidate for the present status of the universe. Anyway, the system is topologically equivalent to the following system at CP1-origin: It means that if A > 0 , then in two directions we have attraction. On the center manifold the related variable, s, is approximately constant and it may be for the point that E imposes a center manifold. Because its equation is satisfied automatically by x and y , not E , so it makes E = 0 for all values of E . In the case A < 0 , the corresponding CP would be a saddle point for all values of m 3 .

CP2:
In order to put CP2 in the physical region, one must set α > −1 and m 3 > 0 . According to the EoS of CP2, it has this potential to be in the matter/radiation/current stage of the universe. Because E = m 3 , thus m 3 should be fixed to the era in which we desire to insert this point.
• If it belongs to the matter-dominated era, then we must set A = 1 which leads to (λ 1 , λ 2 , λ 3 ) = (−4, 3, 3) . For sure, putting A = 1 makes this CP a saddle point which is completely acceptable for a matter-dominated era. Note that, here with the assumed conditions, namely α 1 = 1 , A = 1 , and B = 0 , we have Indeed, under the aforementioned conditions, the viscous fluid is converted to a pressureless matter. Furthermore, ρ v > ρ r at this CP makes it a good candidate for the matter-dominated era.
• If a radiation-dominated era is desired, then A = 4/3 must be adopted. It is worthy to note that at this CP, the density of the radiation component is zero, but under the conditions A = 4/3 , B = 0 , and α 1 = 1 , the viscous fluid is converted to the radiation, namely p v = −ρ v + f + G = ρ v /3 . As a matter of fact, at this CP, the universe is filled only with the radiation. This selection, A = 4/3 , yields (λ 1 , λ 2 , λ 3 ) = (0, 4, 4) . It has a one-dimensional center manifold. If the system be transformed as (x, y, E)| CP 2 −→ (X, Y, R) = (0, 0, 0) and then we use the following transformations, we arrive at the following topologically equivalent system where we exploited the following mappings: v =h 2 (s) = k 4 s 2 + k 5 s 3 + k 6 s 4 + O(s 5 ).
As is clear, the center manifold, s, is approximately constant and this CP is a repeller. Therefore, this CP is not a good point for the radiation era.
• For considering CP2 as a candidate for the current status of the universe, it is better to take A = −0.03 instead of A = 0 because of high accuracy 12 . The case A = −0.03 leads to w eff. | CP2 = −1.03 and (λ 1 , λ 2 , λ 3 ) = (−4.09, −0.09, −0.09) . Therefore, it is an attractor fixed point which is a good option for the present status of our universe, especially with regards to the point that ρ v > ρ r , representing the domination of dark energy for the current position of the universe. It is easily shown that A = 0 for CP2 is also a good candidate for the present time. The case A = 0 yields w eff. | CP2 = −1 and (λ 1 , λ 2 , λ 3 ) = (−4, 0, 0) . Since one of the eigenvalues is negative, hence, without the use of the center manifold theorem, one can state that this point is a saddle or an attractor point. Both are justifiable for the current status of the universe. Therefore, it is a good candidate for the current stage of the universe as well. However, the case A = −0.03 is better than A = 0 .
CP3: For this point, the conditions of being in the physical region are α > −1 and m 3 > 0 . According to the value of the EoS, 1/3 , this point belongs to a radiation-dominated era and thus it must be a saddle point, mandating A > (4/3) to yield (λ 1 , λ 2 , λ 3 ) = ('A negative number ', 4, 4) . In addition, m 3 can be fixed to the radiation era (i.e. m 3 = (κ 2 H 2 r ) −1 in which H r is the Hubble parameter corresponding to the radiation-dominated era). Under the aforementioned conditions, this point may be a good candidate for the radiation era, especially, adherent to the fact that at this CP we have ρ r > ρ v which indicates the domination of radiation component.
In fig. 5, three Poincaré sections of phase space of our system have been presented by setting A = 1.35 , and α = 1 .

Case-2:
As the second example, let us consider the system for the selections α 1 = 2 and β = 1 . The resulting dynamical system and the associated EoS parameter turn out to be ♣ The case B = 0 : In this case, instead of a critical point, we have a critical curve (i.e. E = Cx 2 where C = (−9A/(Bκ 2 )) > 0 ) for a given amount of A/B : where m 1 is an arbitrary constant. Interestingly, all CPs of this curve yield In order to put these CPs in the physical region, three conditions are required: m 1 > 0 , B = 0 , and AB < 0 . The corresponding eigenvalues take the forms: The analysis of the CPs requires the use of the center manifold theorem because λ 2 = 0 . Assuming m 1 = 2(1+α) , this system has a one-dimensional center manifold. Performing a transformation as (x, y, E)| CPs −→ (X, Y, R) = (0, 0, 0) and then applying the transformation our system will be expressed in a new coordinate (u, v, s) . By the use of two mappings we arrive at a system that is topologically equivalent to the original system at its CPs: The center manifold is repeller and hence it makes the CPs unstable and saddle. Setting engenders attraction property in two directions and hence under this condition and according to the unstable nature of our universe, one may conclude that the CPs in which we have set Bκ 2 H 2 0 + 9Am 2 1 = 0 (for putting them in the current time) are good candidates for the present status of the universe. It is worthwhile to mention that for non-zero values of m 1 , one has ρ v > ρ r which refers to the domination of dark energy.
As an example, three different Poincaré sections of the phase space of the above system have been presented in fig. 6 by setting κ 2 = 1 , A = 1 , B = −9 , m 1 = 1 , and α = 1 .
♣ The case B = 0 : Setting B = 0 , reduces the system to the following straight line of critical points: x = 0, y = 0, E = an arbitrary value ≡ m 3 , with the EoS parameter and these eigenvalues According to the value of EoS parameter, this case may be set to the current time when one takes m 3 = (κ 2 H 2 0 ) −1 , but this CP would not be a good case for the present status of the universe because ρ v = ρ r = 0 . Indeed, these CPs correspond to a vacuum universe. Anyway, since λ 1 = −4 < 0 , so this critical line has an attraction in one direction, implying that all CPs of this line would be attractor or saddle points. Let's pause and consider this problem deeply using the center manifold theorem. Clearly, the corresponding center manifolds are two-dimensional. Using two transformations, first (x, y, E)| CP s −→ (X, Y, R) = (0, 0, 0) and then one may get the topologically equivalent system: where we have used the mapping As seen, the plane of the center manifold, v − s , is completely unstable and all CPs are saddle.
As an example, three different Poincaré sections of the phase space of the above system have been presented in fig. 7 in which we have set κ 2 = 1 , A = 1 , m 3 = 1 , and α = 1 .

Case-3:
As the third case, let us consider the system with the choices α 1 = 2 and β = 2 , then the dynamical system takes the following form: The associated EoS parameter also reads: ♣ The case B = 0 : In this case, all the points of x − E plane can be the fixed point to our system: where m 3 is an arbitrary constant, and = ±1 . However, both = +1 and = −1 are acceptable from a mathematical point of view, but only = +1 is acceptable from a physical perspective. Besides this condition, one must have AB ≤ 0 , A = 0 , and m 3 > 0 to maintain these CPs in the physical region. The corresponding EoS parameter at the CP-plane turns out to be w eff. | CPs = −1. The eigenvalues of the Jacobian matrix at the CPs would be According to the value of the EoS parameter, the conditions λ 1 < 0 , m 1 > 0 , and m 3 = (κ 2 H 2 0 ) −1 ( H 0 is the current value of the Hubble parameter) should be adopted. In fact, λ 1 < 0 enhances the number of attractor directions, m 1 > 0 guarantees the domination of dark energy, and the condition m 3 = (κ 2 H 2 0 ) −1 fixes the CP at the current time. Under these conditions, the related CP would be a good case for the present time.
Let us consider this CP deeply using the center manifold theorem. First, we transform the system into the CP-origin (i.e. (x, y, E)| CP −→ (X, Y, R) = (0, 0, 0) ). Then utilizing the relations our system can be transformed into a new coordinate (u, v, s) . Using two mappings h 1 and h 2 as v =h 2 (s) = k 4 s 2 + k 5 s 3 + k 6 s 4 + O(s 5 ), one arrives at the following topologically equivalent system Clearly, for the CP under the aforementioned conditions (i.e. λ 1 < 0 , & m 1 > 0 , & m 3 = (κ 2 H 2 0 ) −1 ), this system indicates attractor behavior for the present stage of the universe which is justifiable. The center manifold, s, is approximately constant, hence it is neither attractor nor repeller. Indeed, a particle on the center manifold is in the rest status so that there is no force on it and the particle continues to move at the speed it has arrived at this point.
As an example, three different Poincaré sections of the phase space of the above system have been illustrated in fig. 8 in which we have chosen κ 2 = 1 , A = 1 , B = −9 , m 1 = 1 , m 3 = 1 , and α = 1 .

Case
The value of B

The number of CPs
The values of the effective EoS of the CPs

An important conclusion of "More complicated forms of the EoS" and a remedy for its main problem
First of all, check Table I. In case B = 0 , the critical points of 'case-1 ' for different values of the constant parameters could span a plane (i.e. x − E ), but in that case, all the points of that plane had w eff. = −1 . Therefore, if one tunes the constant parameters for the present time, the system gives us only one critical point and it may describe only one stage of the universe. In the case B = 0 we saw that it could act in a better way with more fixed points and stages. On the other hand, this occurred only for 'case-1 '. Perhaps, one may argue that such a behaviour stems from the viscous part of the system restricting the behavior of the system and not generating the last stages of the universe in the current model. But, the power of density, α 1 , also may have such behavior if it is greater than 1 . For sure, f (ρ) = Aρ α1=+1 is better than f (ρ) = Aρ α1=+2 in the current model of study. According to the above achievements, it is inspired that the absence of some CPs and stages of the evolution of the universe may be due to the fact that we have considered those cases that their α 1 s are positive. Nevertheless, α 1 with negative values could accentuate the significance of the viscous part. Note that the negative values of β cannot be related to the main problem, because the viscous part, H 2β , is highly likely to be the case at early times, so we deduce that β > 0 would be its representative.
Let us delve deeper into the subject by specifying some values for the constant parameters first. Let α 1 = −1 and β = +1 . Fortunately, these selections lead to an analytically solvable case without further quantifying the constant parameters. The CPs and corresponding eigenvalues and associated EoS parameters are as follows: CP1: CP2: CP3: (x, y, E) = (1 + α, 0, 0) , w eff. | CP3 = 1 3 Clearly, the role of the viscous part is significant. Especially, in CP3, this impression is very obvious as it can be used to put this CP in an arbitrary stage of the evolution of the universe. Furthermore, as is observed, α (the coefficient of torsion in the model) acts as assistance for tuning. CP1: CP1 can be fixed to the current stage of the universe by adopting the following conditions which guarantee almost all aspects of it (domination of dark energy, fitting to the current stage, attractor property, and being in the physical region): CP2: CP2 should be regarded as a candidate for the radiation-dominated era because of the value of the EoS parameter. The conditions reassuring the presence of this CP in the physical region are: Under the condition the domination of the radiation component is satisfied, and the condition makes this CP as a saddle point. The common domain of all these conditions tends to a good CP for the radiation era. CP3: By setting α > −1 , this point falls into a physical region in which dark energy dominates. The status of CP3 is affected by the conditions imposed on CP1 and CP2. After tuning the constant parameters, the situation of this CP is clarified. According to the EoS of this CP, it seems that it is flexible, but one must note that this point should be in a dark-energy-dominated era because of ρ v > ρ r . Furthermore, the viscous fluid could not be converted to other types of components (i.e. matter and radiation) due to the imposed conditions. According to the above explanations and examinations, the main problem has been solved by switching to negative values of α 1 . This change means that the Chaplygin is a better candidate to couple with the viscous part than Van der Waals gas. Last but not least, since the negative values of α 1 have strong correspondence with Chaplygin gas, the best value for α 1 may be in the range −1.8 ≤ α 1 ≤ −0.2 . In ref. [23], it was indicated that the constant power in the generalized Chaplygin gas model, namely ν in P = −σ 2 /ρ ν , might not be absolutely constant, but could be slowly varying with cosmic time. The reason for this theory was to reach the unification theory. More precisely, the power of the density is a time-dependent variable falling so slowly from an initial constant value within 0.2 ≤ ν ≤ 1.8 to zero. However, it is better to take the aforementioned initial domain of ν as |ν| 0.05 when we focus only on the dark sector [23]. These discussions guide us to develop the EoS of viscous fluid by making α 1 and β dependent upon the cosmic time: where σ 1 , σ 2 , and B are constant parameters. Considering such intricate generalized cases are beyond the scope of this paper, hence we terminate this section here.

Some interacting models
In this section, some non-trivial interacting models are studied. To begin with, using a power-law term a modified type of interaction is investigated in detail. Then, the oscillating type of interaction ('Cosine' type) which has been recently introduced by S.D. Odintsov et al. [30] is considered widely. Furthermore, some new types of oscillating interactions are suggested instead of the aforementioned type. It is illustrated that through some enhancements and modifications to the oscillating type, more proportionate results are derived in comparison to that of the basic model. Finally, the benefits of the oscillating types of the interacting models are indicated that are: "Unification between the early inflation and late-time-accelerated expansion" and "The crucial role of the oscillating types of interactions in getting out of some special eras that were under some dominations".
6.1. Interaction through modified term (using some power-law terms) As the first non-trivial case, the model which has been studied in sect. 4 is developed here by utilizing the multiplier H 2j where j is a constant [30]. Let's perform the replacement [30] in the equations of ρ v and ρ dm . Again, let us assume the following form for the bulk viscous coefficient ζ : The dynamical system, in this case, is expressed as follows: 3 Defining the dimensionless variables as we get in which The effective EoS in terms of unknown function U for this case takes the form Now, using f = ρ v (1 + w 0 ) and G = w 1 H 2 , the equations of our dynamical system become: Furthermore, the effective EoS parameter turns out to be: Investigating this case without any assumption is challenging. This behavior was predictable according to the simpler case which has been considered in sect. 4. Hence, we prefer to go about some special cases. Like ref. [30], let us consider three cases: j = 0 , j = −1 , and j = +1 . After studying these cases, some discussions about j ≤ 0 and j > 0 are performed and the behaviors of these general cases are clarified.
6.6.1. The case j = 0 The case j = 0 corresponds to the usual interaction. In fact, when one ignores E , the equations of sect. 4 are recovered, so those points with E = 0 , are common between this case and the one considered in sect. 4. For this case, all CPs and eigenvalues can be reached analytically but they are 'VLT' such that we cannot present them excluding one CP. The presentable CP is as follows: Pay attention to its correspondence with (27). This CP belongs to the radiation-dominated era because of w eff. | CP. = 1/3 . The eigenvalues of this CP are 'VLT' except one eigenvalue λ 1 = +4 . Other eigenvalues have many constant parameters that may be fixed such that they tend to a saddle nature for this CP. Anyway, according to λ 1 = +4 , one may conclude instability for this CP -it can be a saddle or a repeller point (note that saddle is acceptable). It is interesting to note that λ 1 = +4 is the eigenvalue of E -direction and it refers to the expansion nature of the universe and decreasing nature of Hubble parameter with time. It may be concluded that the expansion feature of the universe causes its unstable nature. Other CPs and their eigenvalues and EoS parameters are 'VLT' so that they propose nothing interesting excluding one common property: All other CPs are in x − y plane for which we have E = z = 0 . It means that all the remained CPs belong to the accelerated era and, moreover, they can be used for fixing some relations between ρ v and ρ dm . Indeed, all other CPs are on the surface which describes interactions between dark matter and dark energy. It is notable that on this surface, these CPs show zero value for the density of the radiation. The discussiblity of the CPs of the above system comes through fixing some relations among constant parameters or tuning constant parameters numerically merely. As an example, let us present a case with the fixed condition q = −w 1 κ 2 . Under this condition, four CPs are obtained as follows: where χ 2 = (κ 4 w 2 1 + 6κ 2 w 1 α + 6κ 2 w 1 + 9α 2 − 36αζ 0 + 18α − 36ζ 0 + 9) 1/2 . CP1: It is clear that this CP belongs to the radiation-dominated era because of its EoS. At this point, dark energy has no any effect. For maintaining this point in the physical region, we must set the following condition: Note that after the first part of the above condition, there is " > 0 " while after the second, we see " ≥ 0 ". The reason is that this CP is in the radiation-dominated era hence the associated ρ r = 0 is meaningless. Furthermore, this condition provides the domination of the radiation component.
We have a fixed relation between dark matter and radiation at this CP: According to λ 1 = +4 , this CP is unstable. To obtain a physical behavior for CP1, we must take one of the remained eigenvalues less than zero. This leads to saddle nature for CP1 which is desirable. CP2, CP3, and CP4: The conditions for having these CPs in physical regions are clear; it is sufficient to keep the non-zero coordinates positive. According to the coordinates of these CPs, CP3 can be fixed for the dark-energy-dominated era while CP2 and CP4 are neither suitable for the present stage nor for the radiation-dominated era. In fact, CP2 and CP4 belong to an era in which dark matter dominates. If it is desired to convert the dark matter to pressureless matter, then putting ζ 0 = 0 would suffice by which the pressure of dark matter vanishes, and the interaction between dark matter and dark energy becomes the interaction between a pressureless matter and dark energy. For example, by the use of the selections where m 4 is an arbitrary constant. In this example, dark matter is converted to a pressureless matter by choosing ζ 0 = 0 . For reaching suitable and physical points, it was preferred to remove the interaction by setting q = 0 . It is clear that, for α > −1 , CP1 and CP2 are good candidates for the present time and matter-dominated era, respectively. CP3, however, has EoS = 1/3 and the density of the radiation component is greater than others when α > −1 , but because it is a repeller point, hence it is not a good option for the radiation era. CP4 corresponds to a vacuum universe which is unacceptable. Thus, two salient stages of the universe were achieved.
In the above example, a fixed relation, q = −w 1 κ 2 , caused some presentable solutions. Now let us consider the case j = 0 by specifying some values for constant parameters.
As observed, CP1 with EoS= +1 is a repeller fixed point. At this point, all densities are zero. CP2 is in the radiation-dominated era with saddle nature that is acceptable. Note that at CP2, all densities are zero excluding radiation density ( ρ r = 0 ). This indicates the domination of radiation. Adopting the condition 0 < m 1 < 3 , CP3 can lie in the physical region. According to its EoS, this CP belongs to a matter-dominated era. The related eigenvalues are indicative of a saddle property which is acceptable. Some Poincaré sections of the phase space of this system have been presented in fig. 9. As another example with the above-mentioned data-set, let us generalize our equations by adding the baryonic matter, ρ m , to the collection of densities. Defining we get the following equations The effective EoS parameter for this system then reads w eff. = −1 With the previous selections for parameters, the following CPs with their corresponding eigenvalues and EoS parameters are found: CP1: (x, y, z, g) = (0, 0, 0, 0), w eff. | CP1 = +1, (λ 1 , λ 2 , λ 3 , λ 4 ) = (2, 3, 3, 3) (224) CP2: (x, y, z, g) = (0, 0, 3/2, 0), w eff. | CP2 = +1/3, CP1 and CP2 are similar to the previous example and therefore their interpretations are the same. The effective EoS parameter of CP3 depends upon the values of x ∝ ρ v and g ∝ ρ m . Note that the position of y is affected by x and g . In fact, the density of dark matter is affected by the amounts of densities of both dark energy and baryonic matter. It is clear that the value of EoS is proportional to the density of dark matter. For being in a physical region, the condition 0 ≤ (m 1 + m 2 ) ≤ 3 must be adopted besides m 1 ≥ 0 and m 2 ≥ 0 . Therefore, the dark matter causes the existence of an upper bound on the summation of densities of dark energy and baryonic matter.
According to the elements of coordinates of CP3 and its EoS parameter and eigenvalues, this point cannot be fixed to the radiation-dominated era, but it is a good option for the present time.
In the generalized case, namely by defining E = (κ 2 H 2 ) −1 and adding to our equations, the same results are repeated with E| all CPs = 0 and the following positive eigenvalues, which indicate the expansion nature of the universe evolution and decreasing property of the Hubble parameter: Overall, our data analysis indicates that we cannot set all the constant parameters so that all these CPs belong to different crucial stages of the universe -in matching data, some of CPs do not belong to well-known eras, or some of CPs are repeated, for example, two CPs for the radiation-dominated era, or in most of the cases, some of points ruled out because of strange physical meanings, etc. In a nutshell, after examining many cases under the aforementioned condition (i.e. q = −w 1 κ 2 ), it is concluded that we cannot expect to find a data set so that these CPs illustrate a complete schema of the crucial stages of the evolution of the universe. For the case j = −1 there is a similar situation. However, all of CPs can be found analytically, but they are 'VLT' excluding one CP: The EoS parameter of this CP is 1/3 which indicates that this CP is in the radiation-dominated era, hence it must be a saddle point. One of the associated eigenvalues of this CP is +4 showing its unstable nature. Other eigenvalues are 'VLT', but they could be tuned so that at least one of them will be negative, making it a saddle point then.
Other CPs have a common property: They are on y − z plane because we have x = E = 0 for all of them. As mentioned earlier, the plane of y − z is the plane of the interaction between dark matter and dark energy. Clearly, at these CPs one can fix some relations between the densities of dark matter and dark energy. Moreover, these CPs are in the accelerated era.
For reaching a presentable set of CPs, again there are two ways like the previous case. Note that the number of different CPs is affected by our fixed relations among parameters and our selections for the values of constant parameters.
Let us here we take w 0 = −1 . This condition leads to three presentable CPs with their properties as follows: where χ 3 = (κ 2 w 2 1 − 6κ 2 w 1 α − 6κ 2 w 1 + 9α 2 + 36αζ 0 + 18α − 36ζ 0 + 9) 1/2 . CP1: CP1 is an excellent candidate for the radiation-dominated era. Besides maintaining the coordinates of this point greater than zero, one must note that the value of z must be greater than others for providing the domination of the radiation component. Also, we must care about the amounts of eigenvalues; at least one of λ 1 or λ 2 must be negative because this CP should be a saddle point.
CP2 and CP3: According to the coordinates of CP2 and CP3, these points can be fixed for the accelerated era -from the onset of acceleration of the universe up to now. For this end, besides tuning the values of the EoS parameters, one must take the points (excluding the present point) as a saddle point, for saddle points can lie on the evolutionary trajectory of the universe. The associated CP of a point scribed to the present time will be of attractor or saddle nature. Furthermore, at a CP of this era, the value of x must be greater than the value of y due to the acceleration feature of expansion and the domination of dark energy.

The case j = +1
For the case j = +1 , no solution was found. Even, under special conditions like w 0 = w 1 = 0 , there is no analytical solution for the model. In j = +1 case, E appears in the denominator of the equations of x , y , and z . On the other hand, the evolution equations cannot be satisfied when E = 0 . We know that a solution comes through E = 0 and since it leads to infinity and singularity for our system, hence this system is not able to provide any CP.

Discussion about different values of j
In brief, it is argued that only for j ≤ 0 analytical solutions may be found. It is interesting to mention that among the analytical solutions for the case j < 0 , all solutions are 'VLT' excluding one common solution, (227), which belongs to the radiation-dominated era. The 'VLT' solutions are in the interaction plane of dark matter and dark energy namely x − y plane, for we have E = z = 0 . Hence, the critical points can be used to fix some relations between the densities of dark matter and dark energy at these CPs. Furthermore, these CPs are in the accelerated era. In the case j = 0 , three coordinates x ∝ ρ v , y ∝ ρ dm , and z ∝ ρ r are affected by the interaction part (q ∝ Q ), while in the case j < 0 , we do not observe such affects. It is due to the fact that in the case j = 0 , the contribution of the interaction term in the equations appears without other variables as coefficients and it causes no effect by the amounts of other variables. It serves in the equations as a shift value, while in the case j < 0 , the interaction term appears as qE |j| , and for the satisfaction of the equations we must have E = 0 (imposed by the equation of E ), hence this term has been removed by the coefficient E . Therefore, the interaction does not play any role in the fixed points of the systems with j < 0 .
It is interesting to indicate that the upper bound on j (for the existence of analytical solutions) is increasable from j ≤ 0 to j ≤ 3/2 through changes in the definitions of the dimensionless variables: Combine q with the multiplier H 2j to get hence the bound j ≤ 0 is increased to j ≤ 3/2 . Both bounds demonstrate that the power-law interactions in the related dynamical systems should be of increasing nature with time to get more critical points. The benefit of defining q instead of the aforementioned combination in our case is that the behaviors of other types of interaction (e.g. oscillating types that are studied in the next sub-section) are predictable. For more merits check the next sub-section.

Oscillating dark energy-dark matter interaction
In this sub-section, the oscillating model of interaction is studied. This type of interaction ('Cosine' function in which its argument is proportional to H 2 or equivalently E −1 ) has recently been introduced in ref. [30]. Furthermore, some enhancements about this type of interaction are suggested and examined. The dynamical equations of the cosine type of interaction may be easily obtained by the replacement [30] −QH 2j → −Q cos h 0 H 2 in the equations of the previous power-law interaction. The equations would then be: where cos h 0 H 2 = cos h 0 /κ 2 E . It has been assumed in ref. [30] that the argument of cosine function changes monotonically from π/2 to 0 . The starting point, π/2 , is equivalent to the early inflationary era while the endpoint corresponds to the present time. This kind of behavior can be understood by taking sufficiently small values of the parameter h 0 [30]. Indeed, the cosine model under the above assumptions indicates that at early times (i.e. inflationary era; H 2 ∼ M Pl ; t → 0 ) there is no interaction between dark matter and dark energy, and as the universe ages this interaction starts to grow until the maximum value at H = 0 (equivalently t → ∞) [30]. That is why besides the general case, two asymptotic behaviors are of interest: cos h 0 H 2 = 0 (The inflationary epoch), and cos h 0 H 2 = 1 (The present time and late-time acceleration).
Without assuming any condition, the system of evolution may not be solved analytically. This is because the interaction term, is like interaction H 2j . As earlier mentioned, for positive values of j , there is no analytical solution, while for j ≤ 0 the system is analytically solvable. Hence, by switching the problem could be solved. In this case, cos(h 0 κ 2 E) , all CPs can be achieved analytically. Excluding one CP, all other solutions are 'VLT' and they are in x − y plane (i.e. for all CPs we have E = z = 0 and hence they cannot belong to the radiation-dominated era). The presentable CP with EoS = 1/3 is as follows: The corresponding eigenvalues are 'VLT' except +4 case, implying instability. Due to the existing some constant parameters, this CP can be fixed as a saddle point which is completely acceptable.
It is interesting to note another oscillating case namely In this case, the situation is similar to the previous case. But in this case, the coordinates of CPs are shifted with respect to the above solutions because of the phase difference between 'Cosine' and 'Sine' functions. For instance, the CP of the radiation era would be: As observed, the interaction term does not affect this fixed point while in the previous case it has an active role. It is due to the phase difference between the two types of interactions. Note that since all CPs in both cases have E = 0 , hence the interaction term only in 'Cosine'-type has an effect, not in 'Sine'-type, thanks to sin(0) = 0 (which removes the interaction term from the coordinate of the CP).
Since most of the CPs of the general case were of 'VLT' type, hence we have to consider the problem by one of two aforementioned approaches. As an example, let's set w 0 = −1 for the modified 'Cosine' type of interaction (232). Under this condition, all CPs are presentable: CP2 and CP3: (1 + α) Besides z CP1 > x CP1 ≥ 0 and z CP1 > y CP1 ≥ 0 , at least one of λ 1 or λ 2 must be less than zero to maintain CP1 as a physical CP (a saddle CP) to the radiation-dominated era. According to the coordinates, EoS parameters, and eigenvalues, CP2 and CP3 have this potential to be fixed as candidates for the present status of the universe or even earlier stages up to the inflection point where the acceleration of the expansion of the universe started. Thus, there are many options for these points. But the best seems to be the current situation of the universe (i.e. w eff. = −1.03 ± 0.03 ). To guarantee the domination of dark energy for the current stage of the universe, one must set x| CP > y| CP > 0 . This case corresponds to the inflationary epoch of the universe. The results of this case are exactly the results of sin(h 0 κ 2 E) (i.e. eq. (235) etc). The reason is that all CPs have E = 0 . Nonetheless, since the inflationary era is our objective in this sub-section, hence let us single out some parameters and fix some relations to reach this end. If one takes then the following CPs are found: where = ±1 and η is an arbitrary constant. The effective EoS parameter of these points would then be Clearly, at the limit like (i.e. very close to zero) our goal is fulfilled. Depending upon the values of the parameters, both or one of the CPs can be regarded as a candidate for the inflationary epoch. The corresponding eigenvalues are as follows: Respecting the imposed limitation (i.e. (239)), λ 1 and λ 2 will be negative which undertake the attractor nature of the points -at least, we have attraction in two directions. So, the values of the constant parameters specify the stability/instability of the fixed points according to the values of λ 3 and λ 4 . This asymptotic case corresponds to the case j = 0 in the power-law interaction which has been studied above carefully. Nevertheless, let us consider the possibility of the description of the current status of the universe with non-zero E like ref. [30]. However, the equation governing the evolution of E is dependent and it is not the original and crucial equation. But let us take it into account to see under which conditions our objective is feasible. To this end, we rewrite the equations in the following form where U = (3(1 + w 0 )x + 3y + w 1 κ 2 + 4z)/3 . If U → 0 , then one has Therefore, the CP of this system would be where m 4 taken to be m 4 = (κ 2 H 2 0 ) −1 , H 0 being the current value of the Hubble parameter. These coordinate points stipulates ζ 0 → 0 because of U → 0 . In fact, the current stage of the universe can be achieved when the pressure of the bulk viscous approaches to zero.
As mentioned above, the lack of a general solution makes cos(h 0 H 2 ) = cos(h 0 /(κ 2 E)) of no favor. So, improvements seem necessary. In what follows, some remedies are suggested and finally, we examine it in the system of ref. [30] and show some missed CPs which cover other stages of the evolution of the universe. 1. Using where C is an arbitrary constant, instead of These modified models have at least two benefits: i. This C allows those CPs with E = 0 to be included in the category of CPs. Thus, the number of our CPs can be increased through this modification.
If it is preferred to restrict the evolution range of the argument of cosine and sine functions in the above models such as ref. [30], the range of this limitation is not the same as the range of ref. [30]. The evolutionary equivalent ranges have been presented in Table II. As is clear from Table II, the asymptotic behaviors of all five models are the same. Therefore, the asymptotic solutions that were considered above, hold true for all five types of interactions. As a result, all five types of interaction are successful in describing both early inflation and late-time-accelerated expansion (i.e. the unification of early and late accelerations is achieved through these models).
If we study this generally, then only one CP is achieved in a closed form: The effective EoS parameter of this CP is −1 indicating the present stage of the universe.
• Now, if, for example, we change the model of interaction as then we get to the following CPs: CP2: w eff. | CP2 = 1 3 , x = κ 2 w 1 + q 1 − 3w 0 , y = −(3ζ 0 + q), z = 0; CP3: w eff. | CP3 = w 0 + κ 2 w 1 + q 3x 03 , x = 3ζ 0 + 3w 0 − κ 2 w 1 + χ 5 6w 0 ≡ x 03 y = κ 2 w 1 − 3ζ 0 3w 0 − κ 2 w 1 + q 3w 0 x 03 , z = 0; CP4: w eff. | CP4 = w 0 + κ 2 w 1 + q 3x 04 , where χ 5 = (κ 4 w 2 1 − 6κ 2 ζ 0 w 1 + 6κ 2 w 0 w 1 + 9ζ 2 0 + 18ζ 0 w 0 + 12qw 0 + 9w 2 0 ) 1/2 . Obviously, the number of CPs has been increased and some missed-out stages have been covered in our approach. First of all, note the correspondence between (253) and (255). Unlike x and y , z is different between both, but these two CPs are almost equivalent due to the existence of h 0 . Thus, the CP of the previous interaction can be recovered here. CP2 can be related to the radiation-dominated era because of its EoS parameter. CP3 and CP4 can be fixed to arbitrary stages. It is interesting to note that the effective EoS parameter and the density of dark matter, in both CP3 and CP4, are affected by the value of the dark energy density. Furthermore, the difference between the effective EoS parameter of CP4 and CP3 is proportional to the difference between their dark-matter densities, viz., w eff. | CP4 − w eff. | CP3 = −w 0 (y| CP4 − y| CP3 ) ; as well. Interestingly, the CPs of this type of interaction are similar to the previous type. Hence, the previous discussions are valid here, even the relation (259) between CP3 and CP4. As a final point, note that thanks to the existence of h 0 and C , CP1, i.e. (261), is equivalent to both (253) and (255). Indeed, one solution is common among different types of interacting models.
According to all examinations performed in this paper, we conclude that not only our proposed models are generative of CPs equivalent to that of ref. [30] but also give rise to further CPs. Furthermore, as showed in 'General case', our models yield CPs that are absent in ref. [30]. Therefore, the proposed extensions to interactions seem to be useful. 6.6.5. A nice feature and benefit of the modified oscillating interactions Remember that this CP belongs to the radiation-dominated era. Now, note the amplitude of this oscillating term in x , y , and z . When the density of radiation becomes less than other components, then we exit the radiationdominated era, hence, the interaction term, q cos h 0 H −2 , must be positive for reaching this goal. It means that the boundaries should be corrected as 0 < q cos h 0 H −2 < 1 3 .
These bounds clarify one of the benefits of restricting the arguments of interacting models to positive eras. Clearly, as the universe ages, the effect of this oscillating term on the reduction of the density of the radiation component is higher than dark matter because of the coefficients ( −9 and −1 ). This guarantees exiting the radiation era.
On the other hand, the density of dark energy raises with time under the assumed conditions, and it causes the late-time-accelerated expansion and the domination of the dark energy. The same process may be performed for other types of modified interaction introduced in this paper.

Conclusion
This paper was devoted to investigating a multi-fluid universe using the dynamical system approach in precise detail. The model of this study is a special and well-known form of F (T ) -gravity in FRW background geometry. First, a universe filled with radiation and a perfect fluid with a non-trivial EoS was studied based on the standard approach of dynamical systems. Then, in sect. 4, by adding dark matter to the collection of density-set and assuming usual and well-known interaction between dark matter and dark energy, the system was considered. It is notable that in this study, the bulk viscous pressure of the dark matter fluid was taken into account according to the best fit of observational data. Despite that the interaction was assumed between dark matter and dark energy, the results demonstrate its effect on another component, radiation, as well (It should be noted that this direct effect of interaction on all components, was observed in all interactions studied in this paper). It was stated that excluding one CP which belongs to the radiation-dominated era, other CPs are of VLT-type and all these CPs are on the interaction plane (i.e. they belong to the accelerated era and can be used to fix some relations between the densities). Next, in sect. 5, by employing one extra dimension, E = κ −2 H −2 , and considering the effective equation of state parameter, generalized forms of cosmological fluid were studied in two main parts: 'A simplified form of EoS' and 'More complicated forms of EoS'. In the former one, indeed, at a special limit, the case of sect. 4 was recovered. It was found that we cannot set the constant parameters so that all crucial stages of the universe evolution be given. We also generalized this case by adding baryonic matter to the density-collection and found that the problem still prevails. Also, it was noticed that the eigenvalue of the extra dimension E is positive in all CPs. It refers to the universe expansion and reduction in the Hubble parameter overtime. Furthermore, it was found that the value of E is zero in all CPs excluding those with EoS = −1 . In fact, if one, according to the observational data, takes the current value of the effective EoS parameter equal to −1.03 or any value which is close to it (excluding −1 ), then 0 will be the value we find for E . Indeed, rather than the value of the variable itself, it is the eigenvalues that represent physical information.
In the latter one, a more complicated form of EoS of dark energy candidate, p v = −ρ v + Aρ α1 v + BH 2β , was considered. It was indicated that for different positive values of α 1 and β , all CPs belong to the present time with EoS = −1 . Furthermore, these models were also studied for B = 0 . The solutions of B = 0 were better than B = 0 since they could cover more stages of the universe evolution. The big question in studying the complicated forms was the lack of other stages in case B = 0 . It was illustrated that by switching to negative values of α 1 in B = 0 , the problem may be solved. In fact, this change of α 1 value is equivalent to the shift from Van der Waals to Chaplygin gas. Therefore, more stages of the universe evolution coverage by Chaplygin gas make it a better candidate than Van der Waals to couple with the viscous part. Furthermore, it was argued that the case in which α 1 and β are time-dependent variables can be taken into account. Then, in sect. 6, some interesting models of interaction between dark matter and dark energy were studied. This consideration was done for two types. First, using the multiplier H 2j the usual interaction of sect. 4 was enhanced. This type of interaction has been suggested in ref. [30]. Generally, without any assumption on j , studying this case is challenging. Hence, it was split into three main parts; j = −1 , j = 0 , and j = +1 . After studying these cases analytically and numerically, it was concluded that only for j ≤ 0 there are analytical solutions, not for j > 0 , meaning that the power-law interactions in the related dynamical systems should be of increasing nature with time to get more critical points. The reasons for this event, were discussed in the text. Finally, the cosine