Dynamics of interacting monomial scalar field potentials and perfect fluids

Motivated by cosmological models of the early universe we analyse the dynamics of the Einstein equations with a minimally coupled scalar field with monomial potentials $V(\phi)=\frac{(\lambda\phi)^{2n}}{2n}$, $\lambda>0$, $n\in\mathbb{N}$, interacting with a perfect fluid with linear equation of state $p_\mathrm{pf}=(\gamma_\mathrm{pf}-1)\rho_\mathrm{pf}$, $\gamma_\mathrm{pf}\in(0,2)$, in flat Robertson-Walker spacetimes. The interaction is a friction-like term of the form $\Gamma(\phi)=\mu \phi^{2p}$, $\mu>0$, $p\in\mathbb{N}\cup\{0\}$. The analysis relies on the introduction of a new regular 3-dimensional dynamical systems' formulation of the Einstein equations on a compact state space, and the use of dynamical systems' tools such as quasi-homogeneous blow-ups and averaging methods involving a time-dependent perturbation parameter. We find a bifurcation at $p=n/2$ due to the influence of the interaction term. In general, this term has more impact on the future (past) asymptotics for $p<n/2$ ($p>n/2$). For $p<n/2$ we find a complexity of possible future attractors, which depends on whether $p=(n-1)/2$ or $p<(n-1)/2$. In the first case the future dynamics is governed by Li\'enard systems. On the other hand when $p=(n-2)/2$ the generic future attractor consists of new solutions previously unknown in the literature which can drive future acceleration whereas the case $p<(n-2)/2$ has a generic future attractor de-Sitter solution. For $p=n/2$ the future asymptotics can be either fluid dominated or have an oscillatory behaviour where neither the fluid nor the scalar field dominates. For $p>n/2$ the future asymptotics is similar to the case with no interaction. Finally, we show that irrespective of the parameters, an inflationary quasi-de-Sitter solution always exists towards the past, and therefore the cases with $p\leq(n-2)/2$ may provide new cosmological models of quintessential inflation.


Introduction
The most successful theory of gravity is General Relativity which is based on the Einstein field equations. These equations form a non-linear system of partial differential equations for the dynamics of the metric of a four-dimensional Lorentzian manifold called spacetime. In the case of the standard cosmological models of very large scales, the spacetime metric is commonly assumed to be spatially homogeneous which reduces the Einstein equations to a system of ordinary differential equations. It is remarkable that many recent rigorous results about the dynamics of cosmological models result from the application of methods of dynamical systems, see e.g. [1,2] as well as [3] and references therein. Part of the challenge in applying those methods is the necessity to use conformal rescalings and appropriate normalizations that lead to dimensionless and regular autonomous finite-dimensional dynamical system on compact state spaces. For a review on the theory of dynamical systems in cosmology see [4,5,6].
According to the current cosmological paradigm, the universe went through a period of accelerated expansion soon after the initial (big bang) singularity, called the inflationary epoch [7,8]. Cosmological Inflation, has the ability to solve some of the puzzles of the big bang model, such as the horizon and flatness problems, and also provides probably the best-known mechanism to predict the observable nearly scaleinvariant primordial scalar and tensor power-spectrum [9]. The most popular physical mechanism driving inflation in the early universe consists of a scalar field (called inflaton) with monomial type potentials [10]. A drawback of the original (cold) inflation models is the separation that exists between the inflationary and the reheating periods. The process of reheating is crucial in inflationary cosmology and a natural answer to the problem is the so-called warm inflation first introduced by Berera [11]. In this scenario, the production of radiation occurs simultaneously with the inflationary expansion, through interactions between the inflaton and other fields on a thermal bath. The presence of radiation allows having a smooth transition between the inflationary phase and a radiation-dominated phase without a reheating separation phase. After the inflationary epoch, the standard early universe scenario then involves a period of radiation dominance until a time of decoupling between radiation and matter, after which the universe is dominated by matter, and then galaxies start to form as well as other cosmic structures. The radiation and matter-dominated eras are usually modeled by a perfect fluid with a linear equation of state. It is our purpose here to investigate this dynamical interaction between scalar fields and perfect fluids.
We will consider the Einstein equations for a spatially homogeneous and isotropic metric (the so-called Robertson-Walker metric) having a scalar field with monomial potentials interacting with perfect fluids with linear equation of state. Our main goal is to obtain a global dynamical picture of the resulting system of non-linear ordinary differential equations (ODEs) and in particular of its past and future asymptotics. Our analysis relies on the introduction of a new set of dimensionless bounded variables which results in a regular dynamical system on a compact state-space consisting of a 3-dimensional cylinder. This allow us to describe the global evolution of these cosmological models identifying all possible past and future attractor sets which, as we will see, in many situations, can be non-hyperbolic fixed points, partially hyperbolic lines of fixed points, or even bands of periodic orbits. So, our analysis will require in one hand blow-up techniques and center manifold theory around the non-hyperbolic fixed points and, on the other hand, averaging methods involving a time dependent perturbation parameter.
Our paper is mostly self-contained and is organized as follows: In Section 2 we explain how the non-linear system of ODEs is obtained from physical principles. In Section 3 we find the appropriate dimensionless variables that transform the ODE system into a 3-dimensional regular dynamical system on a compact state space. This construction naturally shows that a significant bifurcation occurs for p = n/2. We therefore split the analysis into three cases according to the different exponents of the scalar field potential and the interaction term: Section 4 treats the case p < n/2, where the future asymptotic analysis is further split in two distinct sub-cases corresponding to p = (n − 1)/2 and p < (n − 1)/2. The case p = n/2 is treated in Section 5 and the case p > n/2 in Section 6. Interestingly, in some situations we encounter non-hyperbolic fixed points whose exceptional divisor of the blow-up space consists of generalised Liénard systems. We provide proofs as well as conjectures about the global dynamics complemented by numerical pictures of representative cases. We present our conclusions in Section 7 where we also mention some physically interesting consequences of our results.

Non-linear ODE system
A spacetime (M, g) is a 4-dimensional Lorentzian manifold with metric g whose evolution is given by the Einstein equations where Ric denotes the Ricci curvature tensor, R the scalar curvature of the spacetime, while T is the energy-momentum tensor which encodes the spacetime physical content. When the spacetime is spatially homogeneous, the Einstein equations became ODEs and can be analysed using dynamical systems' methods. The set of solutions depends crucially on the right-hand-side of the equations, i.e. on the energy-momentum tensor. In turn, this depends on the physical scenario under consideration. Motivated by the warm inflation scenario of the early universe, here we assume a minimally coupled scalar field φ with self-interaction potential V (φ), interacting with a perfect fluid. The evolution equations can be derived from an action principle and the most general action for this case is given by where as standard, we use greek indices µ, ν, ... = 0, 1, 2, 3 for each coordinate in spacetime. Here (∂φ) 2 := (∂ µ φ)(∂ µ φ) with ∂ µ φ = g µν ∂ ν φ whereas L pf is the Lagrangian density of the perfect fluid and L int describes the interaction between the scalar field and the thermal bath. By varying the Lagrangian density with respect to the metric we obtain the Einstein equations (1) with stress-energy tensor components where and u µ denotes the unit 4-velocity vector field of the perfect fluid, with ρ pf > 0 and p pf being the fluid energy density and pressure, respectively. The stress-energy tensor for the scalar field can be written in a perfect fluid form with the identifications u µ (φ) = (∂ µ φ)/ −(∂φ) 2 , The total stress-energy tensor T µν obeys the conservation law ∇ ν T ν µ = 0. However each component of the total stress-energy tensor, T (φ) µν and T (pf) µν , is not conserved, in contrast to the case when the scalar field does not interact with the thermal bath. In the presence of interactions where Q (φ) µ and Q (pf) µ describe the energy exchange between the scalar field and the perfect fluid. It follows from the energy-momentum conservation equation that In this work we consider a phenomenological friction-like interaction term for which where we assume that Γ = Γ(φ) is a function of the scalar field φ only. In more general warm inflationary models, the function Γ can also depend on the thermal bath temperature [14,15,17,18,19], although recent studies suggest that temperature dependence is redundant [22]. Equation (4a) then gives the modified energy "conservation" equation and the Euler equation (ρ pf + p pf )u µ ∇ µ u ν + u ν u µ ∇ µ p pf + ∇ ν p pf = 0. (9b) The above system is closed once an equation of state relating the pressure and the energy density is given.
Here we assume that the fluid obeys a linear equation of state where for example, γ pf = 1 corresponds to a dust fluid, and γ pf = 4/3 to a radiation fluid. The value γ pf = 0 corresponds to the case of a positive cosmological constant and γ pf = 2 to a stiff fluid, both yielding significant dynamical bifurcations. Equation (4b) yields the wave equation where g is the usual D'Alembertian operator associated with the metric g. Motivated by the current cosmological models, as mentioned in the Introduction, we will use a flat spatially homogeneous and isotropic metric g, called Robertson-Walker (RW) metric, that in the Cartesian coordinates (t, x, y, z) ∈ (t − , t + ) × R 3 takes the form g = −dt 2 + a(t) 2 (dx 2 + dy 2 + dz 2 ), where a : (t − , t + ) → R + is a C 2 positive function of time t called scale-factor whose evolution and maximal existence interval will be determined by the Einstein equations. A solution is said to be global to the past (future) if t − = −∞ (t + = +∞). The Einstein equations coupled to the nonlinear scalar field equation (11) and the energy conservation equation for the fluid component (9a), form the following non-linear ODE system for the unknowns {a, H, φ, ρ pf }:ȧ = aH, together with the Gauss (Hamiltonian) constraint where is called Hubble function and a dot denotes differentiation with respect to time t. For expanding cosmologies H > 0. Note also that the equation for the scale factor a(t) decouples leaving a reduced system of equations for the unknowns {H, φ, ρ pf }. The scale-factor can then be obtained by quadrature a = a 0 e Hdt . The Γ(φ) term appearing in (13c) and (13d) acts as a friction term which describes the decay of the scalar field φ due to the interactions encoded in the Lagrangian L int . Here we assume monomial scalar field potentials which are popular examples of inflaton models V (φ) = (λφ) 2n 2n , λ > 0 , n = 1, 2, 3, ...
The exponent 2p reflects the parity invariance of the interaction term, and the condition µ > 0 ensures that the second law of thermodynamics is satisfied (see e.g. [12,20,21]). To summarise, we will analyse the system (13) for the unknowns {H, φ, ρ pf } having the free parameters {n, p, λ, µ, γ pf } besides the initial conditions {ρ 0 , φ 0 ,φ 0 , H 0 }. Note that n, p and γ pf are dimensionless parameters while λ and µ have dimensions. We shall see ahead that it is the dimensionless ratio (see (22) ahead) of these two quantities that plays an important role on the qualitative behaviour of solutions.

Dynamical systems' formulation
In order to obtain a regular dynamical system on a compact state-space, we start by introducing dimensionless variables normalized by the Hubble function H (which is positive for ever expanding models) where c = 6 n−1 n 1 2n λ is a positive constant. We also introduce a new time variableÑ defined by where When δ = 0, i.e. p ≤ n 2 , thenÑ = N = ln(a/a 0 ) is the number of e-folds N from some reference epoch at which a = a 0 , i.e. N = 0. When written in terms of the new variables, the system (13) reduces to a regular local 3-dimensional dynamical system where the constraint equation is used to globally solve for Ω pf . Since Ω pf > 0, the above equation implies that Ω pf is bounded as Ω pf ∈ (0, 1), while Σ φ ∈ (−1, 1) and X ∈ (−1, 1). The positive dimensionless constant ν > 0, is given explicitly by and q is the usual deceleration parameter defined viaḢ = −(1 + q)H 2 , i.e.
where we introduced the Hubble normalized scalar-field energy density and the scalar-field effective equation of state γ φ , defined by Moreover, since γ pf ∈ (0, 2), it follows from (21) and (23) that with limits q = −1 when Σ φ = 0 and Ω pf = 0; q = 2 when X = 0 and Ω pf = 0; and q = 1 2 (3γ pf − 2) when X = 0 and Σ φ = 0. These special constant values of q correspond to well-known solutions: (quasi-)de-Sitter (dS) solution when q = −1, kinaton or massless scalar field self-similar solution when q = 2 and whose scale factor is given by a(t) = t 1/3 , and the flat Friedmann-Lemaître (FL) self-similar solution when q = 1 2 (3γ pf −2) with scale factor given by a(t) = t 2 3γ pf . Although the constraint (21) is used to solve for Ω pf , it is nevertheless useful to consider the auxiliary equation for Ω pf (equivalently Ω φ = 1 − Ω pf ) given by While the variables (X, Σ φ ) are bounded, the variableT becomes unbounded (T → +∞) when H → 0. In order to obtain a regular and global 3-dimensional dynamical system, we introduce so that T ∈ (0, 1) with T → 0 asT → 0, and T → 1 asT → +∞, as well as a new independent time variable τ defined by where This leads to a regular and global 3-dimensional dynamical system for the state-vector (X, Σ φ , T ) where q is given by (23). The auxiliary equation (27) written in terms of the new time variable τ becomes The state space S is thus a 3-dimensional space consisting of a (deformed when n > 1) open and bounded solid cylinder without its axis The state space S can be regularly extended to include the axis of the cylinder with Ω pf = 1 (Ω φ = X 2n + Σ 2 φ = 0) which is an invariant boundary subset as follows from (31), and the outer shell of the cylinder which consists of the pure scalar field boundary subset, Ω pf = 0 (Ω φ = X 2n + Σ 2 φ = 1). Due to the interaction term when ν = 0, the Ω pf = 0 boundary subset is not invariant for the flow. Furthermore at Ω pf = 0 it follows that Since ν > 0, this shows that the surface Ω pf = 0 not being invariant, it is future-invariant, which motivates the following definition: The complement of such orbits in S are said to be of class A.
Class B orbits enter the state-space S by crossing the outer cylindrical shell with Ω φ = Σ 2 φ + X 2n = 1. Moreover S can be regularly extended to include the invariant boundaries T = 0 and T = 1, which is essential, since all possible past attractor sets for class A orbits are located at T = 0 and all possible future attractors for both class A and B orbits are located at T = 1 as follows from the following lemma: Lemma 1. The α-limit set of class A interior orbits in S is located at {T = 0}, while the ω-limit set of all interior orbits in S is located at {T = 1}.
Proof. Since 1 + q ≤ 0, then T is strictly monotonically increasing in the interval (0, 1), except when q = −1 in which case Therefore the points in S with q = −1 are just inflection points in the graph of T (τ ). By the monotonicity principle, it follows that there are no fixed points, recurrent or periodic orbits in the interior of the state space S, and the α-limit sets of class A orbits are contained at {T = 0} and ω-limit sets of all orbits in S are located at {T = 1}.
Thus the global behavior of both classes of orbits can be inferred by a complete detailed description of the invariant subsets T = 0 and T = 1, which are associated with the past and future limits H → +∞ and H → 0, respectively. Due to their distinct properties, we split our analysis into three cases: p < n/2, p = n/2 and p > n/2. 4 Dynamical systems' analysis when p < n 2 When p < n 2 the global dynamical system (30) takes the form where we recall q = −1 + 3Σ 2 φ + 3 2 γ pf Ω pf and the auxiliary equation (31) becomes

Invariant boundary T = 0
The induced flow on the T = 0 invariant boundary is given by Thus and the intersection of the T = 0 invariant boundary with the pure perfect fluid subset Ω pf = 1, i.e. the axis of the cylinder, consists of the fixed point where q = 1 2 (3γ pf − 2), corresponding to the flat FL self-similar solution. The linearisation around FL 0 yields the eigenvalues 3 2n γ pf , − 3 2 (2 − γ pf ) and 3 2n γ pf with eigenvectors the canonical basis of R 3 . Since 0 < γ pf < 2, FL 0 has two positive real eigenvalues and a negative real eigenvalue, being a hyperbolic saddle, and the α-limit point of a 1-parameter set of class A orbits in S.
On the intersection of the invariant boundary T = 0 with the subset Ω pf = 0 there are two equivalent fixed points with q = 2 corresponding to the self-similar massless scalar field or kinaton solution. The linearisation of the full system around these fixed points yields the eigenvalues 3 n , 3(2−γ pf ), and 3 n , with generalised eigenvectors (1, 0, 0), (0, 1, 0) and (∓1, 0, 1). Since γ pf ∈ (0, 2), then K ± are hyperbolic sources and the α-limit points of a 2-parameter set of class A orbits in S.
Finally there are other two equivalent fixed points given by that correspond to quasi-de-Sitter states with q = −1. The linearisation around dS ± 0 yields the eigenvalues −3γ pf , −3 and 0 with eigenvectors (1, 0, 0), (0, 1, 0) and (0, ∓ n 3 , 1) respectively. The fixed points dS ± 0 have two negative real eigenvalues (since γ pf > 0) and a zero eigenvalue corresponding to a center manifold. Due to the monotonicity of T it is clear that a single class A orbit originates from each dS ± 0 into S corresponding to a 1-dimensional center manifold. This center manifold corresponds to what is usually called in the physics literature the inflationary attractor solution, see e.g. [23,24] and references therein. In order to simplify the analysis of the center manifold we use instead system (20) for the unbounded variableT , and introduce the adapted variablesX This leads to the transformed adapted system where the fixed points dS ± 0 are now located at the origin of coordinates (X,Σ φ ,T ) = (0, 0, 0), and F , G and N are functions of higher order terms. The 1-dimensional center manifold W c at dS ± can be locally represented as the graph h : E c → E s , i.e. (X,Σ φ ) = (h 1 (T ), h 2 (T )), satisfying the fixed point h(0) = 0 and the tangency dh(0) dT = 0 conditions (see e.g. [28]). Using this in the above equation and usingT as an independent variable, we get 2 . The problem of finding the inflationary attractor solution amounts to solving the previous system of non-linear ordinary differential equations. Although in general the existence of an explicit solution for the above system is not expected, it is possible to approximate the solution by a formal truncated power series expansion inT : with a i , b i ∈ R. Plugging (45) into (44a)-(44b), using the expansions (X ± 1) 2n = 1 ± 2nX + 2n 2n−2 where δ a b is the Kronecker delta symbol, and solving the resulting linear system of equations for the coefficients of the expansions yields, asT → 0, Therefore, it follows that to leading order on the center manifold which shows explicitly that dS ± 0 are center saddles with a unique class A center manifold orbit originating from each fixed point into the interior of S.
We now show that on T = 0 the above fixed points are the only possible α-limit sets for class A orbits in S, and that the orbit structure on T = 0 is very simple: Then the {T = 0} invariant boundary consists of heteroclinic orbits connecting the fixed points as depicted in Figure 1.
Proof. It is straightforward to check that {Σ φ = 0} and {X = 0} are invariant 1-dimensional subsets consisting of heteroclinic orbits FL 0 → dS ± 0 and K ± → FL 0 respectively. Therefore the two axis divide the (deformed) circle with boundary X 2n + Σ 2 φ = 1, consisting of the heteroclinic orbits K + → dS ± 0 and K − → dS ± 0 , into 4-invariant quadrants. On each quadrant there are no interior fixed points and hence, by the index theorem, no periodic orbits. Since closed saddle connections do not exist, it follows by the Poincaré-Bendixson theorem that each quadrant consists of heteroclinic orbits connecting the fixed points. Moreover, in this case, the {T = 0} invariant boundary admits the following conserved quantity which determines the solution trajectories on T = 0.  Theorem 1. Let p < n 2 . Then the α-limit set for class A orbits in S, consists of fixed points on {T = 0}. In particular as τ → −∞ (N → −∞), a 2-parameter set of orbits converge to each K ± , with asymptotics with C X , C Σ > 0, and C T > 0 constants. A 1-parameter set of orbits converges to FL 0 with asymptotics with C X and C T > 0 constants, and a unique center manifold orbit converges to each dS ± 0 with asymptotics When p = 1 2 (n − 1) we also get from (46d) the asymptotics Ω pf = 2νn 2 27γ pf 1 − 2n 3 N −2 . For p < 1 2 (n − 1) one needs to go higher order on the center manifold of Ω pf .
Proof. The proof follows from Lemmas 1 and 2, and the above local analysis of the fixed points. Remark 1. The solutions of class A which approach K ± behave asymptotically as the self-similar massless scalar field or kinaton solution, and the ones approaching FL 0 as the self-similar Friedmann-Lemaître solution whose asymptotics towards the past exhibit well-known (big-bang) singularities. In the context of early cosmological inflation the physical interesting solution is the center manifold originating from each dS ± 0 whose asymptotics for the variables (H, φ, ρ pf ) of the original system (13), are given by as t → −∞.
with p < n 2 .

Invariant boundary T = 1
On the T = 1 invariant boundary, the system (34) reduces to and the auxiliary equation (35) for Ω pf , satisfies The analysis can be divided into two particular sub-cases: p < parameterised by X 0 ∈ [−1, 1]. In addition to L 1 there exists another line of fixed points when p > 0, parameterised by Σ φ0 ∈ [−1, 1]. We shall refer to the non-isolated fixed point at the origin of the T = 1 invariant set as FL 1 , and the end points of L 1 with X = ±1 as dS ± 1 . The description of the induced flow on T = 1 when p < 1 2 (n − 1) is given by the following simple lemma: Lemma 3. When p < 1 2 (n − 1), the set {T = 1} \ L 1 for p = 0, and the set {T = 1} \ {L 1 ∪ L 2 } for p > 0 are foliated by invariant subsets X = const. consisting of regular orbits which enter the region Ω pf > 0 by crossing the set Ω pf = 0 and converging to the line of fixed points L 1 , as depicted in Figure 2.
Proof. When p < 1 2 (n − 1), the system (53) admits the following conserved quantity   Theorem 2. Let p < 1 2 (n − 1). Then the ω-limit set of all orbits in S is contained on L 1 . In particular: i) If p < 1 2 (n − 2), then as τ → +∞, a 2-parameter set of orbits converges to each of the two fixed points dS ± 1 on the line L 1 with X 0 = ±1, and a 1-parameter set of orbits converges to FL 1 with X 0 = 0.
ii) If p = 1 2 (n − 2), then as τ → +∞, a 2-parameter set of orbits converges to each of the two fixed points S ± on the line L 1 with , and a 1-parameter set of orbits converges to FL 1 with X 0 = 0.
Proof. The first statement follows from lemmas 1 and 3 if p = 0, while for p > 0, it is shown in Lemma 5, Subsection 4.3.1 (by doing a cylindrical blow-up of L 2 on top of the blow-up of FL 1 ), that no solution trajectories converge to the set L 2 \ FL 1 .
The linearised system around L 1 has eigenvalues 0, −νX 2p 0 and 0, with associated eigenvectors (1, 0, 0), (0, 1, 0) and (0, − 2 ν (p − 1 2 (n − 1))X 2(n−p)−1 0 δ n−2p 2 , 1). On the {T = 1} invariant boundary the line of fixed points L 1 is normally hyperbolic, i.e. the linearisation yields one negative eigenvalue for all X 0 ∈ [−1, 1], except at X 0 = 0 when p > 0 (where the two lines intersect), and one zero eigenvalue with eigenvector tangent to the line itself, see e.g. [29]. On the closure of S, the line L 1 is said to be partially hyperbolic [31]. Each fixed point on the line, including the point at the center when p = 0, has a 1-dimensional stable manifold and a 2-dimensional center manifold, while the point with X 0 = 0 is non-hyperbolic for p > 0. In this case the blow-up of FL 1 is done in Section 4.3.1.
To analyse the 2-dimensional center manifold of each partially hyperbolic fixed point on the line, we start by making the change of coordinates given bȳ which takes a point in the line L 1 to the origin (X,Σ φ ,T ) = (0, 0, 0) withT ≥ 0. The resulting system of equations takes the form where F , G and N are functions of higher order. The center manifold reduction theorem yields that the above system is locally topological equivalent to a decoupled system on the 2-dimensional center manifold, which can be locally represented as the graph h : E c → E s withΣ φ = h(X,T ) which solves the nonlinear partial differential equation subject to the fixed point and tangency conditions h(0, 0) = 0, and ∇h(0, 0) = 0, respectively. A quick look at the nonlinear terms suggests that we approximate the center manifold at (X,T ) = (0, 0), by making a formal multi-power series expansion for h of the form Solving for the coefficients it is easy to verify thatã i0 are identically zero, so that h can be written as a series expansion inT with coefficients depending onX, i.e.
, with the origin (0, 0) being a nilpotent singularity. Since the coefficientc 01 (X) = 0 for all X 0 , then the formal normal form is zero with and Φ an analytic function. The phase-space is the flow-box multiplied by the functionT * , with the direction of the flow given by the sign of b 01 , see Figure 3(a). For X 0 = 0 (when p = 0), then b 11 = 3γ pf 2n > 0 and c 01 = − 3γ pf 2n < 0 which after changing the time variable to d/dτ =T −1 d/dτ yields a hyperbolic saddle, see and after changing the time variable to d/dτ =T −1 d/dτ the origin is a hyperbolic sink, see Figure 3(c). For p < 1 2 (n−2), the coefficient b 01 vanishes at X 0 = 0 and X 0 = ±1, being negative for X 0 ∈ (−1, 0) and positive for X 0 ∈ (0, 1). For b 01 = 0 the phase-space is again as depicted in Figure 3(a) with the direction of the flow given by the sign of b 01 , i.e. of X 0 .
When X 0 = 0 (and restricting to p = 0), b 01 = 0 and b 11 = 3γ pf 2n > 0 which after changing the time variable to d/dτ =T −1 d/dτ yields that FL 1 is a hyperbolic saddle, see Figure 3(b). For p > 0 the blow-up of FL 1 can be found in Section 4.3.1, where it is shown that a 1-parameter set of interior orbits in S end at FL 1 , see also Remark 5. For and the origin is a semi-hyperbolic fixed point with eigenvalues −3γ pf , 0 and associated eigenvectors (1, 0) and (− n 3γ pf ν , 1). To analyse the 1-dimensional center manifold W c at (X,T ) = (0, 0) we introduce the adapted variableX =X + n 3γ pf νT . The 1-dimensional center manifold can be locally represented as the graph h : E c → E s , i.e.X = h(T ), satisfying the fixed point h(0) = 0 and tangency dh(0) dT = 0 conditions, usingT as an independent variable. Approximating the solution by a formal truncated power series expansion h(T ) = N i=2 a iT i , a i ∈ R, and solving for the coefficients yields to leading order on the center manifold Therefore for X 0 = ±1, the origin is the ω-limit set of a 1-parameter set of orbits on the 2-dimensional center manifold, see Figure 3(d).
Remark 2. Let p < (n − 1)/2. The asymptotics for solutions of (13) converging to dS ± 1 are given by while those converging to S ± are given by It is possible to deduce the asymptotics for the 1-parameter set of solutions converging towards the fixed point FL 1 when p = 0, while for p > 0, the asymptotics can be obtained by the analysis of the blow-up of FL 1 done in Section 4.3.1, see Remark 5.
The global state-space picture for solutions of the dynamical system (34) when p < 1 2 (n − 1) is shown in Figure 4. The solid numerical curves correspond to the center manifold of dS ± 0 , and the dashed numerical curves to solution curves originating from the source K − . All these solutions end at the generic future attractors dS ± 1 when p < 1 2 (n − 2) or S ± when p = 1 2 (n − 2).

Case
In this case there is a single fixed point lying in the intersection of the T = 1 invariant boundary with the pure perfect fluid invariant subset Ω pf = 1 given by The linearisation around this fixed point on the T = 1 invariant boundary yields the characteristic polynomial λ 2 + νδ n 1 λ + nδ n 1 = 0. Since ν > 0, then for n = 1 (p = 0), FL 1 has two eigenvalues with negative real part being a hyperbolic sink on T = 1 (stable node if ν ≥ 2, and a stable focus if 0 < ν < 2), while on the full state space FL 1 has a 1-dimensional center manifold with center tangent space E c = (0, 0, 1) , i.e. consisting of the Ω pf = 1 invariant subset. Therefore FL 1 is the ω-limit point of a 2-parameter set of orbits, converging to FL 1 tangentially to the center manifold when ν ≥ 2, or spiraling around the center manifold when 0 < ν < 2. When p = 0 and n > 1 all eigenvalues of the fixed point FL 1 are zero. The blow-up of the fixed point FL 1 when p > 0 is done in Subsection 4.3.2.
Lemma 4. Let p = 1 2 (n − 1). Then the {T = 1} invariant boundary consists of orbits which enter the region Ω pf > 0 by crossing the set Ω pf = 0 and converging to the fixed point FL 1 at the origin as τ → +∞, see Figure 5.
Proof. It suffices to note that the bounded function Ω pf is strictly monotonically increasing along the orbits, except at the axis of coordinates Σ φ = 0 or X = 0 when p > 0. However since dΣ φ /dτ = 0 on Σ φ = 0 and dX/dτ = 0 on X = 0, except at origin where the axis intersect, it follows by LaSalle's invariance principle that (Σ φ , X) → (0, 0) and Ω pf → 1. In fact, when p = 1 2 (n − 1), the system (53) admits the following conserved quantity on T = 1: which determine the solution trajectories on the T = 1 invariant boundary.
The global state-space picture for solutions of the dynamical system (34) when p = 1 2 (n − 1) is shown in Figure 6. The solid numerical curves correspond to the center manifold of dS ± 0 and the dashed numerical curves to solution curves originating from the source K − . All these solutions end up at the generic future attractor FL 1 .

Blow-up of FL 1 when p > 0
To analyse the non-hyperbolic fixed point FL 1 of the dynamical system (34) when p > 0, we start by relocating FL 1 at the origin, i.e. we introduceT = 1 − T , to obtain the dynamical system where we recall In order to understand the dynamics near the origin (X, Σ φ ,T ) = (0, 0, 0), which is a non-hyperbolic fixed point for p > 0, we employ the spherical blow-up method, see e.g. [32,27,33]. I.e. we transform the fixed point at the origin to the unit 2-sphere S 2 = {(x, y, z) ∈ R 3 : x 2 + y 2 + z 2 = 1}, and define the blow-up space manifold as B := S 2 × [0, u 0 ] for some fixed u 0 > 0. We further define the quasi-homogeneous blow-up map which after cancelling a common factor u 2p(n−2p) (i.e. by changing the time variable toτ defined by d/dτ = u −2p(n−2p) d/dτ , where recall p < n 2 , with p > 0) leads to a desingularisation of the non-hyperbolic fixed point on the blow-up locus {u = 0}. Since Ψ is a diffeomorphism outside of the sphere S 2 × {u = 0}, which corresponds to the fixed point (0, 0, 0), the dynamics on the blow-up space It usually simplifies the computations if instead of standard spherical coordinates on B, one uses different local charts κ i : B → R 3 such that ψ i : Ψ • κ −1 i and the resulting state vector(-fields) are simpler to analyze. We choose six charts κ i such that where ψ 1± , ψ 2± , and ψ 3± are called the directional blow-ups in the positive/negative x, y, and z-directions respectively. It is easy to check that the different charts are given explicitly by The transition maps κ ij = κ j • κ −1 i then allow us to identify fixed points and special invariant manifolds on different charts, and to deduce all the dynamics on the blow-up space. For example, in this case, we will need some of the following transition charts Since the physical state-space hasT ≥ 0, we are only interested in the region {z ≥ 0}, i.e. the union of the upper hemisphere of the unit sphere S 2 with the equator of the sphere {z = 0} which constitutes an invariant boundary. This motivates that we start the analysis by using chart κ 3+ , i.e. the directional blowup map in the positive z-direction, on which the northern hemisphere is mapped into the invariant plane of coordinates (x 3 , y 3 ) = (x, y, 1). After cancelling a common factor u 2p(n−2p) 3 (i.e. by changing the time we obtain the regular dynamical system where In these coordinates the equator of the sphere is at infinity and it is better analysed using charts κ 1± and κ 2± . Moreover, the above system is symmetric under the transformation (x 3 , y 3 ) → −(x 3 , y 3 ) and, therefore, it suffices to consider the charts in the positive directions. To study the points at infinity, we notice that both the directional blow-ups in the positive x and y directions already tell how such local chart must be given. To study the region where x 3 becomes infinite, we use the transition chart κ 3+1+ : where To study the region where y 3 becomes infinite, we use the chart κ 3+2+ : to obtain the regular dynamical system where The general structure of the blow-up space B for the two different cases with p < 1 2 (n − 1) and p = 1 2 (n − 1) is shown in Figure 7. In the first case, we shall see ahead that the points R ± are still non-hyperbolic and, therefore, we need a further blow-up, while in the second case the number of fixed points on the equator depends on whether ν < 2n, ν = 2n or ν > 2n. The lines of fixed points L ± 1 and L ± 2 correspond to the lines of fixed points L 1 and L 2 on the cylinder state space S with X 0 < 0, or X 0 > 0, and Σ φ0 < 0, or Σ φ0 < 0, respectively, i.e. To obtain a global phase-space description, instead of projecting the upper-half of the unit 2-sphere on the z = 1 plane, we shall project it into the open unit disk x 2 + y 2 < 1 which can be joined to the equator (the unit circle on {z = 0}). In this way we obtain a global understanding of the flow on the Poincaré-Lyapunov unit cylinder D 2 × [0,ū 0 ] for some fixedū 0 > 0. Usually, a generalised angular variable is used on the invariant subset {u 3 = 0}, see e.g. [34,35]. Here we use a different type of transformation, based on [24,25], which makes the analysis somewhat simpler: where is bounded and analytic for θ ∈ [0, 2π), satisfying F (θ) ≥ 1 (with F ≡ 1 when p = 0) and F (0) = √ 2p + 1. The above transformation leads to Making a further change of time variable fromτ 3 to ξ defined by we obtain the dynamical system where On the {ū = 0} invariant boundary, the right-hand-side of the above dynamical system is regular and can be regularly extended up to {r = 0} and {r = 1}, while forū > 0 it can be extended to the invariant boundaries {r = 0} and {r = 1} at least in a C 1 manner. The general structure of the Poincaré-Lyapunov cylinders in both cases when p < 1 2 (n − 1) and p = 1 2 (n − 1) is shown in Figure 8. All fixed points are thus hyperbolic or semi-hyperbolic and, in particular, when p < 1 2 (n − 1) the line L 2 no longer exists, see Figure 8.

Positive z-direction
In this case all fixed points of (85) are located at the invariant subset {u 3 = 0}, where the induced flow is given by and where we introduced the notation Since γ pf ∈ (0, 2), it follows that K ∈ −∞, 1 1+2p , i.e. K < 1 for all p > 0.

Remark 4. By changing coordinates to
where which arises from the second-order Liénard-type differential equation Introducing the functions the energy of the system is E = 1 There is vast amount of literature on Liénard-type systems, see e.g., [27,33,36,37] and references therein. The most difficult problem in this context concerns the existence, number, relative position and bifurcations of limit cycles arising in the Liénard equations.
The fixed points of (96) are the real solutions to In this case there are at most three fixed points. The fixed point at the origin of coordinates whose linearised system has eigenvalues , there is a bifurcation leading to a center manifold associated with a zero eigenvalue.
To analyse the center manifold we introduce the adapted variableȳ 3 The center manifold reduction theorem yields that the resulting system is locally topological equivalent to the 1-dimensional decoupled equation on the center manifold, which can be locally represented as graph , solving the nonlinear differential equation subject to the fixed point h(0) = 0 and tangency dh dȳ3 (0) = 0 conditions, respectively. Approximating the solution by a formal truncated power series expansion, and solving for the coefficients, yields to leading order and therefore the one dimensional center manifold is stable. In addition to M there are two more fixed points when K > 0, whereas no additional fixed points exist when K ≤ 0. When K > 0, the fixed points are The linearisation around these fixed points yields the eigenvalues with associated eigenvectors It follows that S ± are hyperbolic saddles.

Fixed points at infinity
On the positive x-direction, for p < 1 2 (n−1), the system (87) when restricted to the invariant subset {u 1 = 0} reduces to Furthermore, we are interested in the invariant subset {z 1 = 0} on {u 1 = 0}, resulting in which has only one fixed point at the origin The linearisation yields the eigenvalues λ 1 = −ν, λ 2 = 0, and λ 3 = 0 with associated eigenvectors v 1 = (1, 0, 0), v 2 = (0, 1, 0), and v 3 = (0, 0, 1). The zero eigenvalue in the u 1 -direction is associated with a line of fixed points parameterized by constant values of u 1 = u 0 > 0, and which corresponds to the half of the line of fixed points L 1 with X 0 > 0, and denoted by L + 1 , see Figure 7(a). Thus on the {u 1 = 0} invariant set, the fixed point P + is semi-hyperbolic, with the center manifold being the invariant subset {y 1 = 0}, where the flow is given by and so, on {u 1 = 0}, P + is the ω-limit point of a 1-parameter set of orbits. By the symmetry of the system on the (x 3 , y 3 ) plane, an equivalent fixed point P − exists in the negative x-direction.
On the positive y-direction, when p < 1 2 (n − 1), the flow induced on the invariant subset {u 2 = 0} is given by We now focus on the invariant subset {z 2 = 0} on {u 2 = 0}, which results in and has only one fixed point at the origin whose linearised system has all eigenvalues zero. The zero eigenvalue in the u 2 direction is due to the line of fixed points parameterized by constant values of u 2 = u 0 > 0, and which corresponds to the half of the line of fixed points L 2 with Σ φ0 > 0, and denoted by L + 2 , see Figure 7(a). Nevertheless it follows from (110), and (114), that the equator of the Poincaré sphere consists of heteroclinic orbits R + → Q ± , and R − → Q ± .
To blow-up R + and, more generally, the complete line L + 2 ∪ R + , we perform a cylindrical blow-up, i.e. we transform each point on the line to a circle for some fixed u 20 > 0 and s 0 > 0. We further define the quasi-homogeneous blow-up mapΨ and choose four charts such thatψ In fact we only consider the semi-circle with w ≥ 0 since z 2 ≥ 0, which means that we only need to consider the blow-up in the positive w-direction, i.e. the directional blow-up defined byψ 2+ .
We start with the v-direction {v = ±1} which after cancelling a common factor s 2p(n−2p−1) 1± (i.e by changing the time variable d/dτ 2 = s 2p(n−2p−1) 1± d/dτ 1± ) leads to the regular dynamical system . (118) The above system has the fixed points whose linearised system has eigenvalues − ν n−2p−1 , − ν n , and (n−2p) n(n−2p−1) ν with eigenvectors the canonical basis of R 3 , and the fixed points where only Q + − exists in the region w 1± ≥ 0. The eigenvalues of the linearised system around Q + − are 2pν n(2p+1) , ν, and − ν n with associated eigenvectors the canonical basis of R 3 . In the positive w-direction {w = 1}, and after cancelling a common factor s 2p(n−2p−1) 2+ (i.e. by changing the time variable d/dτ 2 = s 2p(n−2p−1) 2+ d/dτ 2+ ), we get the regular dynamical system This system only has the fixed point Q + − which is located at v 2+ = − 2p+1 with associated eigenvectors the canonical basis of R 3 . The blow-up of L + 2 ∪ R + is shown in Figure 9. The blow-up of the equivalent non-hyperbolic set L − 2 ∪ R − follows by symmetry considerations from which we deduce the existence of the equivalent fixed points T − ± , and Q − + . Moreover, since all fixed points

Global phase-space on the Poincaré-Lyapunov disk
The previous results can be collected in a global phase-space by employing the Poincaré-Lyapunov compactification. This compactification has the advantage that all fixed points are hyperbolic or semi-hyperbolic and, in particular, on the Poincaré-Lyapunov cylinder the line L 2 is absent. When p < 1 2 (n − 1), we obtain from (95) that the induced flow on the {ū = 0} invariant subset is given by sin 2θF (θ) cos 2p θ − cos 2(2p+1) θ r .
At {r = 0} lies the fixed point M which is the origin of the (x 3 , y 3 ) plane. The fixed point M is a hyperbolic saddle for K < 0, a center-saddle when K = 0 and a hyperbolic source when K > 0. When K > 0 we have two additional saddle fixed points S ± that are located at The points at infinity in the (x 3 , y 3 ) plane are now located in the {r = 1} invariant set. The hyperbolic sinks P ± and the hyperbolic sources Q + − and Q − + are given by Theorem 4. Let p < 1 2 (n − 1) with p > 0. Then for all ν > 0, the Poincaré-Lyapunov disk consists of heteroclinic orbits connecting the fixed points M, P ± , Q + − , Q − + , and S ± when they exist, with the separatrix skeleton as depicted in Figure 10.
Proof. First notice that {y 3 = 0} is an invariant subset consisting of heteroclinic orbits M → P ± which splits the phase-space into two invariant sets {y 3 > 0} and {y 3 < 0}. On each of these invariant sets there are no fixed points when K ≤ 0, and if K > 0 there is a single fixed point which is a saddle. Therefore by the index theorem there are no periodic orbits in each of these regions. Since closed saddle connections cannot exist either, then by the Poincaré-Bendixson theorem the ω and α-limit sets of all orbits in the Poincaré-Lyapunov disk are the fixed points M, P ± , Q + − , Q − + , and S ± when they exist and the phase-space consists of heteroclinic orbits connecting these fixed points. In particular, from the local stability properties of the fixed points, it is straightforward to show that when K ≤ 0 there are two separatrices Q + − → M and Q − + → M which further split the regions y 3 > 0 and y 3 < 0 into two invariant subsets. The flow on these subsets is trivial. When 0 < K < 1 1+2p there are four separatrices Q + − → S − , M → S − and Q − + → S + , M → S + which further split the invariant regions y 3 > 0 and y 3 < 0 into four invariant subsets where the flow is also trivial.

Remark 5.
It is interesting to obtain the asymptotics for the orbits on the cylinder S towards FL 1 when p < (n − 1)/2. For example when 0 < γ pf < 2n n+1 , there exists a 1-parameter family of orbits in S with the following asymptotics towards FL 1 as τ → +∞ with C ∈ R a constant, which is obtained via the linearised solution at M restricted to the 2-dimensional stable manifold when K < 0. Similarly, one can obtain the asymptotics towards FL 1 when K = 0, by making use of the flow on the center manifold of M, and when 0 < K < 1 1+2p via the linearised solution at S ± restricted to the 2-dimensional stable manifold.

Positive z-direction
Setting p = 1 2 (n − 1) in (85) leads to where Since n is odd, the system is symmetric under the transformation (x 3 , y 3 , u 3 ) → (−x 3 , −y 3 , u 3 ), and all fixed points lie on the invariant subset {u 3 = 0} where the induced flow is given by and K defined in (97) reduces to with K ∈ −∞, 1 n . Since n > 1, it follows that K < 1 for all n. Remark 6. The system (129) can be transformed into the Liénard type system (98), where now the functions f (x) and g(x) are given by The fixed points of (129) are the real solutions to The first equation admits at most five real solutions. The origin of coordinates is always a fixed point M : When K = 0 (γ pf = 2(n−1) n ) there is a bifurcation leading to a center manifold associated with a zero eigenvalue. To analyse the center manifold we introduce the adapted variablesȳ 3 = y 3 − n(1−K) 3(1−nK) x 3 . The center manifold can then be locally represented as the graph h : E c → E u i.e. x 3 = h(ȳ 3 ), which solves the differential equation subject to the fixed point h(0) = 0 and tangency dh(0) dȳ3 = 0 conditions respectively. Approximating the solution by a formal truncated power series expansion and solving for the coefficients yields and therefore M is a center-saddle in this case. The remaining four fixed points that may or not exist, depending on the parameters range, are with where the subscripts ± stand for A ± and the superscripts ± stand for the sign of x 3 at the fixed point.
The eigenvalues λ 1− and λ 2− of the linearisation around S ± − are real for all ν > 2n √ K. Moreover, λ 1− is negative and λ 2− is positive, so that S ± − are hyperbolic saddles. When ν = 2n √ K, with 0 < K < 1/n, it follows that λ 1 and λ 2 are always real and positive, so that the fixed points S ± 0 are unstable nodes.

Fixed points at infinity
We start with the positive x-direction. Setting p = 1 2 (n − 1), the flow induced on the {u 1 = 0} invariant subset is given by We are interested in the invariant subset {z 1 = 0} on {u 1 = 0}, where the system reduces to and yields two fixed points when ν > 2n: where we introduced the notation The linearised system at P + ± has eigenvalues λ 1 = −(n − 1)B ± , λ 2 = ∓2n B±+λ2 . Since P + + has λ 2 , λ 3 < 0 and λ 1 > 0, and P + − has λ 3 < 0 and λ 1 , λ 2 > 0, they are both hyperbolic saddles. But on {u 1 = 0}, P + − is a source and P + + is a saddle. In particular, a one-parameter set of orbits originates from P + − and a single orbit from P + + into the region {z 1 > 0}. When ν = 2n, the fixed points P + ± merge into a single fixed point P + 0 : where the eigenvalues reduce to λ 1 = n − 1, λ 2 = 0 and λ 3 = −1. Hence the z 1 = 0 axis is the center manifold of P + 0 . As before, symmetry considerations allow us to deduce the existence of equivalent fixed points P − ± and P − 0 when blowing-up in the negative x-direction. Finally, when 0 < ν < 2n there are no fixed points on the {z 1 = 0} axis.
Using the directional blow-up in the positive y-direction, the flow induced on {u 2 = 0} is given by Further restricting to the invariant subset {z 2 = 0} results in Since n is odd, there are two fixed points when ν > 2n, corresponding to P − ± studied above which are located at (x 2 , z 2 , u 2 ) = B 1/n ± , 0, 0 and, for ν = 2n, merge into a single fixed point P − 0 . When 0 < ν < 2n there are no fixed points on {z 2 = 0} ∩ {u 2 = 0}. A similar analysis in the negative y-direction shows that the only fixed points are P + ± and P + 0 . Therefore the equator of the Poincaré sphere has fixed points when ν ≥ 2n on the second and fourth quadrants, which are α-limit points for orbits in the northern hemisphere. When 0 < ν < 2n the equator consists of a periodic orbit. In order to study the stability of the periodic orbit at infinity in the (x 3 , y 3 ) plane when 0 < ν < 2n, we shall employ a Poincaré-Lyapunov compactification on the unit disk.

Global phase-space on the Poincaré-Lyapunov disk
When p = 1 2 (n − 1), we obtain from (95) that the induced flow on the {ū = 0} invariant subset is given by At {r = 0} lies the fixed point M which is the origin of the (x 3 , y 3 ) plane. The previous analysis showed that M is a hyperbolic saddle if K < 0, a center-saddle if K = 0, and a hyperbolic source if K > 0. The fixed points S ± ± are located at while the infinity of the (x 3 , y 3 ) plane is now located at {r = 1}. When ν > 2n there are four fixed points corresponding to P ± ± and given by which merge into two fixed points P ± 0 when ν = 2n, with which we have seen to be α-limit points for interior orbits on the disk. Finally, when ν ∈ (0, 2n) there exists a periodic orbit: Lemma 6. For all 0 < ν < 2n, with n > 1, n odd, Γ ∞ is an unstable limit cycle.
Proof. To analyse the periodic orbit we introduce the variable s = (1 − r)/r, so that the periodic orbit is now located at s = 0, and a new time variable ξ by d/dτ 3 = (1−r) r d/dξ which leads to a regular ODE that, close to s = 0, is given by where Denoting by s(θ, s 0 ) the solution of the above differential equation such that s(0, s 0 ) = s 0 , then close to where β 1 and β 2 solve the initial value problem The solutions are with The Poincaré return map near s = 0 is P (s 0 ) = s(2π, s 0 ). Since P (0) = 0, and for n odd, P (0) = e α(2π) < 1 and θ is strictly monotonically decreasing, then the periodic orbit Γ ∞ is an unstable limit cycle for all 0 < ν < 2n and n > 1.
Proof. The proof follows by the local analysis of the fixed points P ± ± , P ± 0 and Lemma 6. Figure 11, shows the three different types of orbit structure at the invariant boundary {r = 1} of the unit disk.
, then for all ν > 0 the ω-limit set of all interior orbits on the Poincaré-Lyapunov disk is contained on the set M ∪ S ± + . In particular, as ξ → +∞, exactly 2 orbits converge to the fixed point M and a 1-parameter family of orbits converges to each fixed point S ± + , the separatrix skeleton being as depicted in Figure 11.
Proof. From Proposition 1 every regular orbits on the (x 3 , y 3 ) plane remains bounded for all future times. The divergence of the vector field (129) yields (a) 0 < ν < 2n. which for K ≤ −1/n and ν > 0 does not change sign and vanishes at a set of measure zero. Therefore by the Bendixson-Dulac criteria there are no interior periodic orbits. Moreover the origin M is a saddle and S ± + are sinks for all ν > 0 (focus if ν < f − (K, n), or nodes if ν ≥ f − (K, n)). Since closed saddle connections cannot exist, it follows by the Poincaré-Bendixson theorem that the only possible ω-limit sets in this case are the fixed points S ± + and M. The last statement follows from the local stability properties of the fixed points.
Proof. Here we make use of the equivalent system (102) and apply Liénard's Theorem, see e.g. [27]. The functions g(x) and F (x) = is monotonically increasing to infinity asx → +∞.
Therefore the system has a unique stable interior limit cycle Γ in . Since in this case ν < 2n √ K < 2n, then by Lemma 6 the infinity consists of an unstable limit cycle. Moreover M is a hyperbolic source and therefore the only possible ω-limit set is the unique interior stable limit cycle Γ in .
Remark 8. We note that Liénard's theorem also gives the relative location of the interior stable limit cycle Γ in . Figure 13 shows the Poincaré-Lyapunov disk for K ∈ 0, 1 n and 0 < ν < 2n √ K, where the orbits accumulate at the interior stable limit cycle Γ in .  Remark 9. We now briefly discuss the missing case K ∈ (0, 1 n ) with ν ≥ 2n √ K. Numerical results suggest that, in this case, an interior stable limit cycle exists if ν ≤ (1 + nK) √ n, with S ± + being sources when ν < (1 + nK) √ n and centers with an unstable outer periodic orbit if ν = (1 + nK) √ n, while no interior periodic orbits exist when ν > (1 + nK) √ n. Recall that when K > 0, the fixed point M is a source and S ± − are saddles when they exist, i.e. when ν > 2n √ K. If ν = 2n √ K, then S ± + and S ± − merge into the fixed points S ± 0 , which are unstable strong focus. See Figure 14, for some representative cases. Most of the results about the existence of limit cycles for the Liénard system rely on the strong assumption thatxg(x) > 0 forx = 0, i.e. that the fixed point at the origin is the only fixed point of the system, see e.g. [38] and references therein. For recent works where such assumption is relaxed see e.g. [39] and references therein. Nevertheless, they do not seem enough to prove the numerical results discussed in Remark 9. This is the case, for instance, of Theorem 3.4 of [39] Remark 10. It is interesting to obtain the asymptotics for the orbits on the cylinder S towards FL 1 . For example when 0 < γ pf < 2n n+1 , there exists a one parameter family of orbits in S with the asymptotics towards FL 1 as in Remark 5 with n − 2p = 1.

Invariant boundary T = 0
The induced flow on the T = 0 invariant boundary is given by and from the auxiliary equation for Ω pf , we have that Therefore on the T = 0 invariant boundary, the system (158) admits five fixed points. The fixed point at the origin, with Ω pf = 1, is given by where q = 1 2 (3γ pf − 2) corresponds to the flat Friedmann-Lemaître solution, and whose linearisation yields the eigenvalues 3 2n , − 3 2 (2 − γ pf ) and 3 2n with associated eigenvectors the canonical basis of R 3 . This fixed point has one negative and two positive real eigenvalues, being a hyperbolic saddle from which originates a 1-parameter family of Class A orbits in S (see Definition 1).
On T = 0, the subset Ω pf = 0 (X 2n + Σ 2 φ = 1) is not invariant (except at Σ φ = 0 or X = 0), but it is future-invariant as follows from (162). On this subset there are four fixed points. The first two equivalent fixed points are given by and correspond to massless scalar field states (or kinaton states) with q = 2. The linearisation around these fixed points yields the eigenvalues 3 n , 3 n , and 3(2 − γ pf ), with associated eigenvectors the canonical basis of R 3 . It follows that K ± are hyperbolic sources, so that from each K ± originates a 2-parameter family of Class A orbits in S. The other two equivalent fixed points are and correspond to a quasi-de-Sitter state with q = −1. The linearisation yields the eigenvalues −3γ pf , −(3 + ν) and 0 with eigenvectors (1, 0, 0), (0, 1, 0) and (0, ∓ n 3+ν , 1). These fixed points have two negative real eigenvalues and a zero eigenvalue, possessing a 2-dimensional stable manifold contained in the boundary T = 0, and a 1-dimensional center manifold (the inflationary attractor solution). Just as in section 4.1, the monotonicity of T implies that dS ± 0 are center-saddles with a unique orbit, the center manifold orbit, entering the state-space S. However, due to its relevant physical meaning, it is important to obtain approximations for the center manifold solution.
In order to analyse the center manifold of dS ± 0 we use instead system (20) for the unbounded variableT for p = n 2 , and introduce the adapted variables which put the fixed points dS ± 0 at the origin of coordinates (X,Σ φ ,T ) = (0, 0, 0). This leads to the system where F , G and N are functions of higher-order terms. The 1-dimensional center manifold W c at dS ± 0 can be locally represented as the graph h : E c → E s , i.e. (X,Σ φ ) = (h 1 (T ), h 2 (T )), satisfying the fixed point, h 1 (0) = 0 = h 2 (0), and the tangency conditions dh1(0) dT . Hence usingT as an independent variable, we get Finding the attractor solution amounts to solve the above system of non-linear ordinary differential equations. We can however approximate the solution by performing a formal power series expansion where a i , b i ∈ R. Inserting (169) into (168) subject to the fixed point and tangency conditions, and solving the resulting linear system of equations for the coefficients, results in , Therefore, it follows that to leading order on the center manifold which shows explicitly that dS ± 0 are center saddles with a unique class A center manifold orbit originating from each fixed point into the interior of S.
We now show that on T = 0 the fixed points FL 0 , K ± and dS ± 0 are the only possible α-limit sets for Class A orbits in S, and that the orbit structure on T = 0 is very simple: Then the {T = 0} invariant boundary consists of heteroclinic orbits connecting the fixed points and semi-orbits crossing the set Ω pf = 0 and converging to dS ± 0 , as depicted in Figure 15. Proof. The proof is identical to the proof of Lemma 2, although in this case there are no conserved quantities, and the set Ω pf = 0 is not invariant but future invariant, except at K ± and dS ± 0 .
Theorem 7. Let p = n 2 . Then the α-limit set for class A orbits in S consists of fixed points on {T = 0}. In particular as τ → −∞, a 2-parameter set of orbits converges to each K ± , a 1-parameter set orbits converges to FL 0 , and a unique center manifold orbit converges to each dS ± 0 . Proof. The proof follows from Lemma 1, Lemma 7 and the local stability analysis of the fixed points. Remark 11. When p = n/2, the asymptotics for the inflationary attractor solution originating from dS ± 0 are given by as t → −∞.

Invariant boundary T = 1
When p = n 2 , the induced flow on the T = 1 invariant boundary is given by This system has only one fixed point at Ω pf = 1 given by Lemma 8. Let p = n 2 . Then the {T = 1} invariant boundary is foliated by periodic orbits P Ω φ characterized by constant values of Ω φ , centered at FL 1 , see Figure 16.
Proof. On {T = 1} the auxiliary equation for Ω φ reads from which the result follows.
Proof. The proof relies on Lemma 1 together with generalised averaging techniques based on the methods introduced in [24,26], see also [40] and [41]. Standard averaging techniques and theorems can be found in [42] for the periodic case and [43] for the general case. In these theorems a key role is played by a perturbation parameter . Roughly, a differential equation of the form dx dτ = f (x, τ, ) with > 0 is approximated by the truncated averaged equation at = 0 , i.e. dȳ dτ = f (x, ., 0) . In the present situation the role of -parameter is instead played by the function (τ ) = T (τ ) − 1.
Therefore, we have to prove an averaging theorem for the case where is not a constant, but a variable that slowly goes to zero. Each periodic orbit on T = 1 ( = 0) has an associated time period P (Ω φ ), so that for a given real function f , its average over a time period associated with Ω φ is given by In what follows we will need to compute several averaged quantities such as Σ 2 φ and X n Σ 2 φ . We therefore use a different formulation which is better adapted to the problem at hand. So we introduce new polar variables (r, θ) that solve the constraint equation Ω φ = Σ 2 φ + X 2n , where Ω φ can be seen as the square of the radial coordinate, i.e. r = Ω φ , and (X, Σ φ ) = (r 1 n cos θ, rG(θ) sin θ), where G(θ) satisfies G(θ) ≥ 1 (with G ≡ 1 when n = 1) and G(0) = √ n, see [24,25]. The resulting system of equations takes the form and At = 0, it follows that and therefore θ is strictly monotonically increasing in r ∈ (0, 1]. In this case, the average of a real function f over a time period, associated with Ω φ = const., is given by where Γ[x] is the usual Γ-function. From (178) we obtain where Ω φ = r 2 . Note that the above implies Σ 2 φ = n X 2n , which is in accordance with the result in [24] obtained by averaging the dynamical system. In particular, it follows that the scalar field equation of state (181a) has an average while for the interaction term (181b), we get both independent of Ω φ . The general idea of the averaging method is to start with the near identity transformation and then prove that the evolution of the variable y is approximated, at first order, by the solutionȳ of the truncated averaged equation. The evolution equation for y is obtained using equations (179a) and (179c), leading to Setting where, for large times, the right-hand-side is almost periodic and has an average that is zero so that the variable g is bounded. Then we can use the fact that (1 + ∂g ∂y ) −1 ≈ 1 − ∂g ∂y + O( 2 ) to get dy dτ = f (y) + 2 h(y, g, τ, ) where and h(y, g, , τ, Dropping the higher-order terms in in (189), we study the truncated averaged equation coupled to an evolution equation for : This system has a line of fixed points at = 0 which can be removed by introducing a new time variable where now the = 0 invariant subset has the isolated fixed point whose linearisation yields the eigenvalues 3(γ pf − γ φ ) and −3γ pf /2, with associated eigenvectors (1, 0) and (0, 1) respectively. When γ pf > γ φ , and since σ φ ∈ (0, 1) and ν > 0, there is a second fixed point The linearised system at F 2 has eigenvalues with associated eigenvectors v 1 = (1, 0) and v 2 = (0, 1). Hence F 2 is a hyperbolic sink while F 1 is a hyperbolic saddle. Notice that in the absence of interactions, i.e. when ν → 0, the fixed point F 2 has the limitȳ = 1 as in [23,24]. When γ pf = γ φ , F 2 merges into F 1 , leading to a center manifold with flow given by dȳ/dτ = −3/2ν σ φ ȳ 2 , so that the solutions converge to F 1 tangentially to the = 0 axis. When γ pf < γ φ , F 1 is the only fixed point being a hyperbolic sink.
Next, we prove that the solutions y of the full averaged system (189) have the same limit as the solutions y of the truncated averaged equation when τ → +∞ and, hence, so does r and subsequently Ω φ . In order to do that we define sequences {τ n } and { n }, with n ∈ N, as follows with lim τ n = +∞ and lim n = 0, since (τ ) → 0 as τ → +∞. Notice that for small enough, y is monotone and bounded and, therefore, has a limit as τ → +∞. Then, we estimate |η(τ )| = |y(τ ) −ȳ(τ )|, where y and y are solution trajectories with the same initial conditions as for τ > τ n+1 , and where C and M are positive constants. By Gronwall's inequality Hence for τ − τ n ∈ [0, 1/ n ], i.e. τ ∈ [τ n , τ n+1 ], it follows that with K a positive constant. Letting n → +∞ implies that η → 0 as τ → +∞. Therefore y andȳ have the same limit as τ → +∞, i.e. the fixed point F 1 or F 2 . Finally, from equation (187), using the triangle inequality and the fact that → 0 as τ → +∞, it follows that r (and hence also Ω φ ) has the same limit as y.
The global state-space picture when p = n/2 for the three different future asymptotic regimes is shown in Figure 17. The solid numerical curves correspond to the center manifold solutions of dS ± 0 , and the dashed numerical curves to orbits originating from the source K − . 6 Dynamical systems' analysis when p > n 2 When p > n 2 the global dynamical system (30) reduces to where we recall q = −1 + 3Σ 2 φ + 3 2 γ pf Ω pf and the auxiliary equation for Ω φ becomes where now the effective interaction term σ φ is given by (202)

Invariant boundary T = 0
When p > n 2 the induced flow on the T = 0 invariant boundary reduces to and Thus, the subset Ω pf = 0 is not invariant but future invariant except at Σ φ = 0 or X = 0, which are the points of intersection of the subset Ω pf = 0 with the lines of fixed points with X 0 ∈ [−1, 1] and L 2 : with Σ φ0 ∈ [−1, 1]. We shall refer to the non-isolated fixed point at the origin of the T = 0 invariant set as FL 0 = L 1 ∩ L 2 , the end points of L 1 with X = ±1 as dS ± 0 , and the end points of L 2 with Σ φ0 = ±1 as K ± . The description of the induced flow on T = 0 is given by the following lemma: Lemma 9. When p > n 2 , the set {T = 0} \ {L 1 ∪ L 2 } is foliated by invariant subsets X = const. consisting of regular orbits which enter the region Ω pf > 0 by crossing the set Ω pf = 0 and converging to the line of fixed points L 1 as τ → −∞, see Figure 18.
Proof. When p > n/2 the system (203) admits the following conserved quantity  Theorem 9. Let p > n 2 . Then, the α-limit set of all orbits in S is contained in the set dS ± 0 ∪ FL 0 ∪ K ± . In particular, as τ → −∞, a 2-parameter set of orbits converges to each fixed point K ± , a 1-parameter set of orbits converges to FL 0 and a single orbit converges to each of the fixed points dS ± 0 . Proof. By Lemma 1, the α-limit set of all orbits in S is located at T = 0 and the description of this boundary is given in Lemma 9. We start by analysing the line of fixed points L 1 . The linearised system around L 1 has eigenvalues 0, −νX 2p 0 and 0, with associated eigenvectors (1, 0, 0), (0, 1, 0) and ( 3(2p−n) , 0, 1). On the {T = 0} invariant boundary, the line of fixed points L 1 is normally hyperbolic, i.e. the linearisation yields one negative eigenvalue for all X 0 ∈ [−1, 1], except at FL 0 where the two lines intersect with X 0 = 0, and one zero eigenvalue with eigenvector tangent to the line itself, see e.g. [29]. On the complement of S, the set L 1 \ FL 0 is said to be partially hyperbolic, see e.g. [31]. Each fixed point on this set has a 1dimensional stable manifold and a 2-dimensional center manifold, while the fixed point FL 0 with X 0 = 0 is non-hyperbolic. To analyse the 2-dimensional center manifold of each partially hyperbolic fixed point on the line, we make the change of coordinates given bȳ which takes a point in the line L 1 to the origin (X,Σ φ ,T ) = (0, 0, 0) withT ≥ 0. The resulting system of equations takes the form where F , G and N are functions of higher-order. The center manifold reduction theorem yields that the above system is locally topological equivalent to a decoupled system on the 2-dimensional center manifold, which can be locally represented as the graph h : E c → E s , i.e.Σ φ = h(X,T ), which solves the nonlinear partial differential equation Solving for the coefficients of the expansion one sees that all coefficients of typeã i0 are identically zero, so that h can be written as a series expansion inT with coefficients depending onX, i.e.
where for example After a change of time d/dτ =T 2p−n d/dτ , the flow on the 2-dimensional center manifold is given by For X = 0, the coefficient b 01 only vanishes at X 0 = ±1, being negative for X 0 ∈ (−1, 0) and positive for X 0 ∈ (0, 1). When b 01 = 0 the origin (0, 0) is a nilpotent singularity. Since the coefficient c 01 (X) = 0 for all X 0 , then the normal formal form is zero with and Φ an analytic function. The phase-space is the flow-box multiplied by the functionalT * , with direction given by the sign of b 01 , see Figure 19(a). When X 0 = ±1 we have that b 11 = −3γ pf < 0, b 02 = 0, c 01 = 0, c 02 = 0 and c 12 < −3γ pf . After changing the time variable to d/dτ = T −1 d/dτ , we obtain so that the origin is a semi-hyperbolic fixed point with eigenvalues −3γ pf , 0, and associated eigenvectors (1, 0), and (− n 3γ pf ν , 1). To analyse the center manifold we introduce the adapted variableX =X + n 3γ pf νT . The 1-dimensional center manifold W c at (0, 0) can then be locally represented as the graph h : E c → E s , i.e.X = h(T ), satisfying the fixed point h(0) = 0 and tangency dh(0) dT = 0 conditions, usingT as an independent variable. Approximating the solution by a formal truncated power series expansion and solving for the coefficients yields to leading order on the center manifold So for X 0 = ±1, the origin is the α-limit set of a single orbit, the inflationary attractor solution, see Figure 19(b). Therefore on the set L 1 \ FL 0 only the fixed points dS ± 0 with X 0 = ±1 are α-limit sets for Class A interior orbits in S, being unique center manifolds.
The conclusions about the non-hyperbolic fixed point FL 0 can be found in Section 6.2, where the blow-up of FL 0 is done as well as the blow-up of the line L 2 . In particular, Lemma 10 in Section 6.2.2, states that no interior orbit in S converges to the set L 2 \ {K ± ∪ FL 0 } and that the fixed points K ± are sources.  Remark 12. When p > n/2, the asymptotics for the inflationary attractor solution originating from dS ± 0 are given by

Blow-up of FL 0
To analyse the non-hyperbolic fixed point FL 0 we use the unbounded dynamical system (20), which for p > n 2 gives where we recall In order to understand the dynamics near the origin (X, Σ φ ,T ) = (0, 0, 0), which is a non-hyperbolic fixed point we employ the spherical blow-up method. I.e. we transform the fixed point at the origin to the unit 2-sphere S 2 = {(x, y, z) ∈ R 3 : x 2 + y 2 + z 2 = 1} and define the blow-up space manifold as B := S 2 × [0, u 0 ] for some fixed u 0 > 0. We further define the quasi-homogeneous blow-up map where ψ 1± , ψ 2± and ψ 3± are called the directional blows ups in the positive/negative x, y, and z-directions respectively. It is easy to check that the different charts are given explicitly by The transition maps k ij = κ j • κ −1 i will allow to identify special invariant manifolds and fixed points on different local charts, and to deduce all dynamics on the blow-up space. In this case we will need some of the following transition charts: Similarly to the blow-up of FL 1 , we are only interested in the region {z ≥ 0}, i.e. the union of the upper hemisphere of the unit sphere S 2 and the equator of the sphere {z = 0} which constitutes an invariant boundary subset. This motivates that we start the analysis by using chart κ 3+ , i.e. the directional blow-up map in the positive z-direction, on which the the northern hemisphere is mapped into the plane z = 1 and the equator of the sphere is at infinity, which is better analysed using the charts κ 1+ and κ 2+ . Figure 20 shows the blow-up space of FL 0 when p > n 2 . Later, instead of projecting the upper-half of the unit 2-sphere

Positive z-direction
We start with the positive z-direction {z = 1} which after cancelling a common factor u 2p(2p−n) 3 (i.e. by changing the time variable d/dÑ = u 2p(2p−n) 3 d/dτ 3 ) leads to the regular dynamical system where All fixed points are located at the invariant subset {u 3 = 0}, where the induced flow is given by and where we introduced the notationK together with the change of time variable d/τ 2 = z 2 d/dτ 3 , i.e. d/dÑ = u 2p(2p−n) 2 d/τ 2 , which leads to the regular dynamical system The induced flow on the invariant subset {u 2 = 0} is given by which has only one fixed point R + : The linearised system at R + has all eigenvalues equal to zero. One of these zero eigenvalues is due to the line of fixed points L + 2 in the u 2 direction, which corresponds to the half of the line L 2 with Σ φ0 > 0, see Figures 20 and 21.
To blow-up the non-isolated set of fixed points L + 2 ∪ R + , we perform a cylindrical blow-up, i.e. we transform each point on this set to a circle We choose four charts such thatψ and recall that only the half circle with w ≥ 0 is of interest since z 2 ≥ 0 and, therefore, only the blow-up map in the positive w-directionψ 2+ needs to be considered. We start with the v 1 -direction {v 1 = ±1} which, after cancelling the common factor s 2p(2p−n) 1± (i.e. by changing the time variable d/dτ 2 = s 2p(2p−n) 1± d/dτ 1± ), leads to the regular dynamical system where In the physical region w ≥ 0, the above two systems have only one fixed point each, given by and whose linearisation gives the eigenvalues ν n , 0 and − ν n , with associated eigenvectors the canonical basis of R 3 . The fixed points T + ± have the center manifold s 1± = 0, i.e. the w 1± -axis. On the center manifold, the flow is just so that T + ± are center-saddles. In the positive w-direction and after cancelling the common factor s 2p(2p−n) 2+ (i.e. by changing the time variable d/dτ 2 = s 2p(2p−n) 2+ d/dτ 2+ ), leads to the regular dynamical system where The above system has only the fixed point whose linearisation gives the eigenvalues , with associated eigenvectors the canonical basis of R 3 . So, on the {u 2 = 0} subset, Q + is a hyperbolic source.
The eigenvalues of the linearised system around K + are 3 2p , 3 2np and 3(2 − γ pf ), with associated eigenvectors the canonical basis of R 3 . Since 0 < γ pf < 2, all eigenvalues are real and positive so that K + is a hyperbolic source. The blow-up of L + 2 ∪ R + is shown in Figure 21. Due to the symmetry of the system in the (x 3 , y 3 ) plane, the blow-up of the equivalent non-hyperbolic set L − 2 ∪ R − follows by symmetry considerations yielding identical results, and therefore equivalent fixed points Q − ± and N − . Hence we have the following result: Lemma 10. No interior class A orbit in S converges, for τ → −∞, to the points on the set L 2 \ {FL 0 ∪ K ± }, while a two-parameter set of class A orbits in S converges to each K ± . where F (θ) = 1 − cos 2(2p+1) θ 1 − cos 2 θ .
Proof. First notice that {y 3 = 0} and {x 3 = 0} are invariant subsets consisting of heteroclinic orbits M → P ± and Q ± → M respectively. These separatrices split the phase-space into four invariant subsets (corresponding to the four quadrants). Since on these quadrants there are no fixed points, then there are also no periodic orbits and the result follows by the Bendixson-Poincaré theorem.

Invariant boundary T = 1
The induced flow on the T = 1 invariant boundary is given by and this system has only one fixed point at Ω pf = 1 given by Lemma 11. Let p > n 2 . Then the {T = 1} invariant boundary is foliated by periodic orbits P Ω φ characterized by constant values of Ω φ , and centered at FL 1 , as depicted in Figure 24.
The global state space S for p > n/2 can be found in Figure 25, where the solid numerical curves correspond to the center manifold of each dS ± 0 and the dashed numerical curves to solutions originating from one of the equivalent hyperbolic sources K − . All solutions converge to the global future attractor which is either a limit cycle for γ pf ≥ 2n n+1 or the Friedmann-Lemaître fixed point FL 1 when γ pf < 2n n+1 .

Concluding remarks
Cosmological models are commonly modeled by flat Robertson-Walker metrics which contain either a perfect fluid with linear equation of state p pf = (γ pf − 1)ρ pf , γ pf ∈ (0, 2), or a scalar field φ with a particular form of the potential V (φ). Here, motivated by warm inflation models of the early universe, we considered that both matter models are present, interact with each other and furthermore allow for a large family of potentials V (φ) = (λφ) 2n 2n , n ∈ N, λ > 0, with a friction-like interaction term of the form Γ(φ) = µφ 2p , p ∈ N ∪ {0}, µ > 0. This brings new mathematical challenges as new non-linearities arise in the resulting ODE system comparing to the non-interacting case with µ = 0, discussed in detail in [24]. Our analysis relied on reformulating the Einstein equations as a new regular 3-dimensional dynamical system on a compact state-space. This construction showed that a significant dynamical bifurcation occurs for p = n/2, and that the positive dimensionless parameter ν = n6 2p−(n−1) µ λ n , with ν = 0 when µ = 0, plays a fundamental role on the qualitative description of the models. The analysis was then split into the three distinct cases p < n/2, p = n/2, and p > n/2. For all these models we were able to prove rigorously their main global dynamical properties which we complemented with numerical pictures. Our results also have interesting physical consequences: First, as is common in interacting cosmologies where there is an energy exchange between scalar fields and fluids, see e.g. [15,20,16], there are solutions for which the perfect fluid energy density is generated via energy transfer from the scalar field, thereby starting its evolution at some finite time t * , where ρ pf (t * ) = 0 with H(t * ) > 0 finite, and ρ pf (t) > 0 for all time of existence t > t * . We called these solutions Class B, while Class A solutions have past asymptotics associated with the limit H → +∞.
For Class A solutions, and towards the past, we showed that irrespective of the model parameters an inflationary attractor solution always exist and that it is associated with a 1-dimensional unstable center manifold of a quasi-de-Sitter fixed point, while generically the asymptotics are governed by the massless scalar-field or kinaton solution. However, the influence of the interaction term on the inflationary leading order dynamics differs for the three different cases p < n/2, p = n 2 or p > n 2 , see Remarks 1, 11, and 12. For Class A and B solutions, and towards the future, the asymptotics are more intricate: When p < 1 2 (n − 2) the generic future attractor is the de-Sitter solution, a results that seems unknown in the literature and that might offer a new physical model of quintessential inflation. When p = 1 2 (n − 2) the generic future attractor corresponds to new scaling solutions where the deceleration parameter satisfies −1 < q S ± < −1 + 3γ pf 2 , which interestingly, for some parameter values can also lead to late time accelerated expansion (q S ± < 0). When p = 1 2 (n − 1) the scale factor behaves as the Friedman-Lemaître (FL) solution as t → +∞, while the matter fields have simple asymptotics when (p, n) = (0, 1) depending on if ν < 2, ν = 2 or ν > 2, while for p > 0 there are various kinds of asymptotic behaviour which can be obtained by the analysis of generalised Liénard type systems. When p = n 2 the future asymptotics are either fluid dominated if γ pf ≤ 2n n+1 in which case the scale-factor behaves as the FL solution, or has an oscillatory behaviour, if γ pf > 2n n+1 , where neither the fluid nor the scalar field dominates. Finally when p > n 2 the future asymptotics are similar to the non-interacting case with ν = 0.
Overall, the interaction term has a larger impact on the future asymptotics when p < n/2, and on the past asymptotics when p > n 2 .