{'fraction_non_alphanumeric': 0.039219061201406474, 'fraction_numerical': 0.014876195628119699, 'mean_word_length': 4.224149004495826, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scaling solutions from interacting fluids

We examine the dynamical implications of an interaction between some of the fluid components of the universe. We consider the combination of three matter components, one of which is a perfect fluid and the other two are interacting. The interaction term generalizes the cases found in scalar field cosmologies with an exponential potential. We find that attracting scaling solutions are obtained in several regions of parameter space, that oscillating behaviour is possible, and that new curvature scaling solutions exist. We also discuss the inflationary behaviour of the solutions and present some of the constraints on the strength of the coupling, namely those arising from nucleosynthesis.


I. INTRODUCTION
There has been a growing appreciation of the importance of the asymptotic behaviour of cosmological models [1]. Indeed, unless there is a case for a cosmic coincidence [2], the features of the dynamics should be associated with some stationary regime which should be obtained without fine-tuning of the initial conditions.
A particular case which has attracted a great deal of interest concerns the possibility of obtaining cosmological scaling solutions with self-interacting scalar fields. It was shown that the exponential potentials yield the remarkable feature that the dynamics of the scalar field self-adjusts to that of matter so that the corresponding energy densities become proportional [3][4][5][6][7]. In these solutions one envisages the interplay between the scalar field and matter, instead of focusing on the dynamics of models exclusively dominated by a selfinteracting scalar field as was the case in most of the models of inflation. It was shown that these solutions attract all the other phase space trajectories in the case of flat space models [7] and hence provide the late-time asymptotic behaviour for the scalar field cosmologies under consideration. This gives, on the one hand, a possible answer for why a non-vanishing scalar field does not introduce radical changes with respect to the usual Einstein-de Sitter rate of expansion of the universe. On other hand, it may contribute an explanation to the difference between the actual density of matter and the critical energy density of the spatially flat isotropic models. Furthermore, the scalar field component would also fulfill the convenient role of delaying the time of matter-radiation equality which would help fitting the power spectrum of large scale structure [8][9][10].
So far the emphasis has been placed on the role of given families of potentials. In the literature we find mainly two sorts of potentials underlying the scaling behaviour: the exponential potentials and a class of power law potentials [2,3,11] respectively. However it is worth noticing that the solutions corresponding to the latter set of potentials were only shown to hold in the regime where the perfect fluid component fully dominates the expansion and the energy density of the scalar field is negligible. Regarding the solutions with an exponential potential the self-adjustment of the behaviour of the energy density of one or more scalar fields with that of matter has been investigated in Friedmann-Robertson-Walker (FRW) models [3][4][5][6][7]12,13] both with and without curvature, in spatially homogeneous, but anisotropic models [14], and in FRW models in scalar-tensor gravity theories (also referred to as non-minimal coupling) [15][16][17][18]. It was also shown by two of us [19] that every positive and monotonous potential which is asymptotically exponential yields a scaling solution as a global attractor.
Recently Billyard and Coley [20] have included an interaction term which transfers energy from the scalar field to the matter fields. In this work we analyse a related question, namely the role of the interaction between two of the components of the universe in promoting the scaling behaviour. We show that in the case of the self-interacting scalar field the self-tuning is a direct result of the energy transfer between the purely kinetic part and the vacuum-like part of the field accomplished by the gradient of the potential. We generalize this case by allowing for an interaction between two otherwise barotropic fluids and we examine its implications for the cosmological behaviour of homogeneous and isotropic universes. We consider a phenomenological interaction term between two of the matter components with the form ∝ H λ ρ α 1 ρ β 2 , where α, β and λ are constants which on dimensional grounds have to satisfy λ + 2(α + β − 1) = 0, and which, as we will show below, extends the case found in the scalar field with an exponential potential. Such a kind of energy exchange between two components guarantees the global conservation of energy-momentum, that is, satisfies the contracted Bianchi identities and is akin to the type of dissipative terms considered in single component models where either a bulk viscosity term or matter creation are considered [21][22][23][24][25]. We show that the consideration of such an interaction term yields a variety of situations where scaling solutions emerge. We produce a full classification of those cases where non-trivial scaling solutions (NTSS) are attractors. They include as particular case the scalar field cosmologies already referred to. Moreover we show that oscillating behaviour is also possible. In this latter case although the mean expansion rate is of a power-law type the various matter components oscillate. This paper is organized as follows. In the next section we review the scaling behaviour associated with scalar field cosmologies. We show how we may view this effect as the result of an interaction between the two limit situations where the energy of the scalar field lies in its kinetic part and the alternative case where it lies in its vacuum-like part associated with its potential. We show how the transfer of energy between these is promoted by the gradient of the potential. We also consider scalar field cosmological models within the extended framework of non-minimally coupled gravity theories. The Brans-Dicke (BD) scalar-tensor theory is remarkable in that in the conformally transformed Einstein frame the original coupling between the BD scalar field and the curvature scalar of the space-time is traded into a re-defined scalar with an exponential potential which is now coupled to matter. In scalartensor theories there is then an intrinsic interaction between the conformally transformed scalar field energy density and the matter fields energy densities which lies at the root of the scaling behaviour of some solutions. In fact we consider the matter fields to be a combination of a perfect fluid and a radiation fluid and find a new type of scaling behaviour both in the Einstein and Jordan frames.
In section 3 we generalize the scalar field models to a combination of three matter components where one of them is a perfect fluid and the other two are interacting. We perform a qualitative analysis which reveals that non-trivial scaling solutions arise in a number of situations and we classify them.
Finally in section 4 we summarize and discuss in the context of our general model a number of issues such as inflationary behaviour, curvature scaling solutions, bounds from nucleosynthesis and a phenomenological approach to the decay of massive particles out of equilibrium.

A. Minimal-coupling theories
We consider a homogeneous and isotropic flat Friedmann-Robertson-Walker (FRW) universes filled with a perfect fluid, characterized by p = (γ − 1)ρ, where 0 ≤ γ ≤ 2 is a dimensionless constant, and a self-interacting scalar field within the framework of Einstein's theory of general relativity. The field equations then read where the overdots stand for the derivatives with respect to the time t, and H =ȧ/a (we use units in which c = 1 = 8π/m 2 p , where m p = ( √ G) −1 is the Planck mass and G the gravitational constant). Another equation which is useful, albeit not independent of the former, isḢ Defining ρ K =φ 2 /2 and ρ V = V (ϕ), the scalar field equation (2) may be cast as a system of two equationsρ which shows how the gradient of the potential promotes an interaction between two limit "perfect fluids", namely one associated with the scalar field's kinetic energy, ρ K , that may be characterized by γ K = 2 and hence is akin to a stiff fluid, and the other associated with the potential energy, ρ V , which has a vacuum like character and γ V = 0. In the absence of this interaction, when V is constant (the case of an effective cosmological constant), there is a smooth evolution from the early time domination of the perfect fluid with a higher value of γ to the late time domination of the perfect fluid with the lower γ [26,28]. In this case, there are no scaling solutions as the asymptotic solutions occur with the vanishing of one of the perfect fluid components. Let us further define ρ c and ρ d as which are, respectively, the energy density and the pressure of the self-interacting scalar field [26] in the co-moving observer frame. It then follows from Eqs. (5,6) The second equation (10) is the evolution equation for the scalar field pressure and we see that it is this one which explicitly involves the interaction between ρ 1 and ρ 2 as discussed above.
When we take V (ϕ) to be an exponential potential, i.e., V (ϕ) = V 0 exp (νϕ), where ν is a constant, and introduce τ = ln a 3 as the new time variable, equations (9) and (10) become We further introduce the new variables x and y defined as which correspond to the density parameters associated with ρ c and with ρ d , respectively. Notice that in previous works on scaling solutions with exponential potentials [7,12,13] it became popular to use expansion normalized variables which are square roots of density parameters. The only reason for that choice is to avoid the √ ρ K which appears in the Eqs. (5,6). Since, in the next sections, we shall be interested in models where the interaction term will involve other powers of the individual energy densities, apart from the power 1/2, it is preferable to use the density parameters themselves given their immediate connection to observations. With the definitions (13), (14) we obtain where δ = −ν/ √ 3 and we have made use of the fact that the Friedmann constraint equation for the flat models, Eq. (1), yields Ω = 1 − x − y.
It is straightforward to see that we just need to consider the triangle in the x, y plane bounded by the invariant lines x = y, x = −y and x = 1. The vertices of this triangle are trivial critical points. The origin (0, 0) represents the vanishing of the scalar field energy density and hence of its pressure as well (ρ c = ρ d = ρ K = ρ V = 0, the universe scales as a ∝ t 2/(3γ) and ρ ∝ a −3γ ). The point (1, 1) represents the case where the kinetic energy of the scalar field dominates and thus corresponds to the massless field case (ρ = ρ V = 0, the universe scales as a ∝ t 1/3 and ρ K ∝ a −6 ). Finally, the point (1, −1) represents a spurious solution introduced by the change of variables we performed, for which ρ = ρ K = 0 and ρ V = 3H 2 . Clearly, this is not a solution of the original equations (1,2,3). Notice though that in cases where the potential exhibits an underlying non-vanishing vacuum energy this critical point corresponds to the late-time domination of the vacuum energy and we have the well-known de Sitter exponential behaviour [26,27].

B. Non-minimally coupled scalar fields
Scaling solutions in non-minimally coupled theories have also been found in the literature [15][16][17][18]. These theories, which can be formulated as general scalar-tensor gravity theories, are based on the Lagrangian [32][33][34][35] where R is the usual Ricci curvature scalar of a spacetime endowed with the metric g ab , Φ is a scalar field, ω(Φ) is a dimensionless coupling function, U(Φ) is a function of Φ, and L m is the Lagrangian for the matter fields. They provide a most natural generalization of Einstein's general relativity (GR), and their investigation enables a generic, model-independent approach to the main features and cosmological implications of the unification schemes which involve extra-dimensions. A distinctive feature of these theories is the coupling of the dynamical scalar field Φ to the scalar curvature R, which implies that the gravitational constant is now a function of Φ, in fact G = Φ −1 . The archetypal of these theories is Brans-Dicke theory (BD), for which ω(Φ) is constant [36]. The Lagrangian (17) corresponds to the so-called Jordan frame, in which the matter fields satisfy the equivalence principle (hence their energy-momentum tensor satisfies ∇ b T ab = 0). This means that the L m terms do not explicitly involve the scalar field Φ. By means of an appropriate conformal transformation of the space-time metric g ab to the so-called Einstein frame, where Φ * is a constant allowing the normalization of Newton's constant in the latter frame, we recover a minimally coupled theory where the coupling of Φ to the curvature is traded into a coupling of the redefined scalar field with the matter fields. In fact Φ → ϕ [37,38] such that where α = ( 2ω(ϕ) + 3) −1 . The field equations for the flat FRW universes with a perfect fluid reduce then to the simple form [39] where the overdots stand for the derivatives with respect to the conformally transformed timet,Ṽ (ϕ) = Φ 2 * U(Φ(ϕ))/(8πΦ 2 (ϕ)),M (ϕ) = µ (Φ(ϕ)/Φ * ) −(2−3γ/2) and µ is the constant defined by µ ≡ ρa 3γ which fixes the initial conditions for the energy-density of the perfect fluid (note that for the sake of making clear the effect of the conformal transformation on the coupling we have written the latter equations with the 8π/Φ * in spite of of our choice of units; in what follows we shall again set it equal to 1). From Eq. (21) and the definition of M(ϕ) it is apparent that if the matter sources are radiation fields (the case in whichM = µ is constant) or, alternatively, vacuum (the case in which µ = 0 and henceM = 0) the coupling vanishes, which translates the fact that the latter cases are conformally invariant.
It is a simple matter to see from Eq. (19) that in the BD case, where ω (and hence α) is constant, Φ ∝ exp √ 2 αϕ so thatM is exponential and so isṼ if the original potential U(Φ) is a power-law in Φ. This means that the coupling encapsulated in the M(ϕ) function amounts in this case to a modification of the constant coefficients of the dynamical system associated with the general relativistic case yielding scaling solutions for suitable values of the parameters [15][16][17][18].
Here we shall restrict ourselves to the original BD theory where, besides having a constant ω (or equivalently α), we also have a vanishing U(Φ) in Eq. (17). We will continue to focus on the flat Friedmann-Robertson-Walker model, but instead of just having a perfect fluid we will also consider the simultaneous presence of radiation (hereafter we shall drop the tildes referring to the conformal frame quantities). Denoting the perfect fluid (respectively, radiation) energy density and pressure by ρ m and p m (respectively, ρ r and p r ), and the energy density and the pressure of the redefined scalar field ϕ by ρ ϕ and p ϕ , the field equations are and the evolution equations for the energy densities of the three fluids are given by rendering explicit the interaction that now leads to a transfer of energy between the matter perfect fluid and the stiff fluid (φ 2 /2) associated with the massless scalar field ϕ. Since p m = (γ − 1)ρ m and given the definitions of ρ m and of ρ ϕ , we see that the interaction term has the same form as in the case of the minimally coupled scalar field with an exponential potential. Indeed, we may write for each component an equation of state of the form p m = (γ − 1)ρ m − p int and p ϕ = ρ ϕ + p int , where the interaction term can be cast as There are two major differences with regard to the minimally coupled case. On one hand, the dimensionless constant δ that measures the strength of the interaction depends on γ according with On the other hand, the γ-index of the non-interacting fluid is now fixed to take the value 4/3 of radiation and γ 1 remains equal to 2, and γ 2 now takes the value of the free parameter γ. Therefore taking into account these differences the dynamical system (15,16) adopts the form It is a simple matter to see that as before we obtain two NTSS. One at the point (x = 1,ȳ = δ 2 /2 −1) which corresponds to a late-time domination of the coupled components, the scalar field and the γ-fluid. Radiation is depleted and the universe scales as power-law x where the three components are simultaneously present and the energy densities of the two that are linked by the interaction have adjusted to the radiation behaviour. Thus the scale factor evolves as a ∝ t 1/2 as in a universe dominated by radiation, but with the major difference that both the BD-scalar fluid and the γ-fluid remain in relevant proportions (this scaling solution corresponds in fact to the C RM fixed point in the paper of Amendola [40]).
If we reverse the conformal transformation we are able to recover the exact solutions in the Jordan frame (which is commonly taken to be the physical frame). The case of the fixed point on the line x = 1 corresponds to the power-law solutions found by Nariai [41] and includes the dust solution (γ = 1) derived by Brans and Dicke [36]. The other non-trivial scaling solution corresponds again to a power-law in the Jordan frame, but differs from the single-fluid exact solutions found in the literature [42,43]. For this we have where a and t are the scale factor and the time coordinate in the Jordan frame. This exact power-law solution is new and it differs from the corresponding behaviour of the radiation BD models for which a ∝ t 1/2 .

III. COSMOLOGIES WITH INTERACTING FLUIDS
In this section we shall construct a general model which includes the two previous classes of models as particular cases. We assume the matter content of the universe to be the combination of three components. One that evolves without interaction with the other two and which satisfies the usual barotropic p = (γ − 1)ρ law (in the particular case of the non-minimal coupling models of the last sub-section it would correspond to the radiation component). The two remaining components are mutually coupled by an interaction term of the type where i, j = 1, 2, and α, β, λ are constants, and η is a dimensionless constant. The ± sign of the interaction term means that if it takes one of the signs for one of the components it necessarily takes the opposite sign for the other, thus ensuring the overall conservation of the energy-momentum tensor ∇ b (T ab (1) +T ab (2) ) = 0 as required by the Bianchi contracted identities. Moreover, using dimensional considerations, it is a simple matter to verify that α, β and λ have to satisfy λ + 2(α + β − 1) = 0 for the interaction to be a pressure. The dependence of the interaction term on the products of powers of the densities of the particle species reflects, on one hand, the fact that one expects that it should be proportional to the number of collisions between the particles of each species and thus on the product of their number densities, and, on the other hand, the fact that in the collisions the energies of the particles will be shifted. Furthermore the factor depending on H accounts for the characteristic time of the individual interactions. Since a detailed relativistic kinetic model of the interaction between two fluids is missing, we keep our model as general as possible by using the two free parameters of the set α, β and λ to allow for unknown aspects of the interaction. The model under consideration aims at a phenomenological description of transient phases of the universe which arise when some of the material components are not in thermal equilibrium. One must be wary that this is in fact a likely situation as the expansion of the universe is permanently trying to pull the matter fields out of equilibrium [21,44]. In this sense the type of interaction generalizes those found in the literature when viscous pressure and matter creation terms are considered in the literature on dissipative isotropic models [22][23][24]21,25]. Since our main purpose is to produce a classification of the cases which lead to scaling solutions, where two or more components self-adjust their behaviour so that their energydensities scale with the same rate, we shall be essentially concerned with the dynamical aspects of this model, rather than with the specifics of particular models. Nevertheless we emphasize that the post-inflationary reheating period, the decay of massive particle species into lighter ones and the general situations of decoupling of particle species provide examples of the situations where the present model may be applied. We will address some of these examples in the section IV.
We shall denote the energy density of the first component ρ X and assume that it has a barotropic equation of state characterized by a constant adiabatic index γ X which we shall leave as a free parameter.
The field equations read Here ρ j , p j , with j = 1, 2 are, respectively, the energy density and the pressure of the j th matter component measured by a comoving observer. We define so that the field equations become It is also convenient to introduce the density parameters as dynamical variables and the new time variable τ = ln (a/a o ) 3 . Then, from Eq. (43), ρ X /3H 2 = 1 − x and the field equations reduce to the following dynamical system where the prime denotes differentiation with respect to the τ variable. It is remarkable that the dimensional relation that must be satisfied by the parameters, λ + 2(α + β − 1) = 0, is precisely the one that renders the above system autonomous. From (43), (48) and (49), we have the restrictions 0 ≤ x ≤ 1, |y| ≤ x which define the phase space of the system, a triangle in the x, y plane bounded by the invariant lines x = y, x = −y and x = 1. Using the symmetry (y, γ d , α, β) → (−y, −γ d , β, α) we may restrict ourselves to the when γ d > 0, i.e., γ 1 > γ 2 . Rescaling time through the factor γ d and defining the new (50) and (51) may be simplified and rewritten in a more compact form aṡ The parameter δ measures the strength of the interaction, and |ǫ| < 1 indicates that γ X lies within the interval (γ 1 , γ 2 ). The cases of the previous section correspond to α = 1/2, β = 1. In the minimally coupled scalar field cosmologies, ǫ = γ − 1, whereas in the BD case ǫ = 2−3γ 3(γ−2) . The points (0, 0), (1, 1) and (1, −1) are always equilibrium points for this system. They correspond to the trivial scaling solutions which exist when η = 0. We shall say that an equilibrium point of (52) is a non trivial scaling solution (NTSS) if it is none of the above and furthermore it is stable. Our study of system (52) will be directed towards the search for non trivial scaling solutions. The additional symmetry (y, t, ǫ, α, β) → (−y, −t, −ǫ, β, α) is useful to reduce the number of different possible cases. The techniques we employ are the standard methods of qualitative theory of planar systems, assisted, in the degenerate cases, by numerical integration.
A. Dynamics on x = 1 When ρ X = 0, the interplay between the two interacting fluids is described by the single equationẏ It is straightforward to check that NTSS on x = 1 show up in the three cases depicted in Figure 1. Denoting generically by (x,ȳ) the coordinates of any of these fixed points, we see from Eqs. (43), (46), (48) and (49) that the scaling solution is characterized by a ∝ t 2/(3Γ) where Γ = γ c + γ dȳ .

B. Global dynamics for |ǫ| ≥ 1
When |ǫ| ≥ 1 there are no equilibria of (52) outside the boundary of the phase space. For ǫ < −1, i.e., γ X < γ 2 < γ 1 , the origin is a global attractor. For ǫ > 1, the line x = 1 becomes the global attractor, so that the stable equilibria described in the preceding subsection become NTSS for the full system, which behaves as shown in Figure 2.
When the non-interacting matter component behaves as one of the two interacting fluids, we have the degenerate cases ǫ = 1 (γ X = γ 1 ) and ǫ = −1 (γ X = γ 2 ). These cases also provide NTSS on the lines y = x and y = −x, respectively, as shown in Figure 3 for the cases when other equilibria on x = 1 coexist. In the cases of Figures 3.a) and 3.b), that is, when γ X = γ 2 and α > 1, ρ X would have been negligible in the past, while in the case of Figure 3.c), that is when γ X = γ 1 and β < 1, ρ X would have dominated in the past.
In the case of Figure 4.b), we have the new feature of a neutrally stable NTSS. In the case of Figure 4.b), the basin of attraction of the NTSS shares the phase space with the basin of attraction of (1, −1). D. The case |ǫ| < 1, α, β < 1 and α + β ≥ 1 As above, an equilibrium in the interior of the phase space exists only when (1+ǫ The dynamics in this case is as shown in Figure 5.a). With respect to Figure 2.c), the change in ǫ has turned the NTSS into the only attractor. When In the first case, there are no NTSS on the boundary of the phase space, and we have at (x, ǫx),x given by (54) a neutrally stable (resp. asymptotically stable) NTSS when ǫ = ǫ M (resp. ǫ < ǫ M ), see Figures 5.b) and 5.c). Notice that, in both situations, (0, 0) is an attractor for a set of positive measure of initial conditions. In the second case, the dynamics is as shown in Figures 5.d) and 5.e), according to whether ǫ < ǫ M or ǫ > ǫ M . In both situations, the NTSS is a global attractor.
E. The case |ǫ| < 1, α, β < 1 and α + β < 1 In this case, we shall have an equilibrium in the interior of the phase space at (x, ǫx with x given by (54) whenever (1 + ǫ) 1−α (1 − ǫ) 1−β > δ. When this happens, there are always equilibria on the line x = 1, which falls into the case of Figure 1.c), and the additional equilibrium is a saddle. The dynamics is shown in Figure 6.a). With respect to Figure 2.c), the NTSS survives as a partial attractor.
When (1 + ǫ) 1−α (1 − ǫ) 1−β ≤ δ there are no equilibria in the interior of the phase space and we may still have either Only the first case may provide NTSS, when, moreover, ǫ > 0, and the dynamics is then as shown in Figure 6.b). With respect to Figure 2.c), the dynamics is only slightly changed, and the NTSS.
F. The case |ǫ| < 1 and α ≤ 1 < β or α < 1 ≤ β An equilibrium in the interior of the phase space exists if and only if (1+ǫ) α−1 (1−ǫ) β−1 > 1/δ. In this case, it is easy to check that the equilibrium on x = 1 is a saddle, and that the additional equilibrium is a sink and hence a NTSS. In the complementary case, there can only be non trivial scaling solutions when the dynamics on x = 1 is as in Figure 1.b). Then, the NTSS on x = 1 is a sink for the full system, and a global attractor. The dynamics is represented on Figure 7. Finally, it is easy to check using the symmetries that the case β ≤ 1 < α provides no NTSS.
G. The case |ǫ| < 1, α = β = 1 and δ > 1 In this case, we have also a neutrally stable NTSS in the interior of the phase space. The dynamics is shown in Figure 8. With respect to Figure 5.b), the behaviour is similar, except that (0, 0) is now a saddle so that, in the present case, there are no attractors.

IV. DISCUSSION AND CONCLUSIONS
In this work we have produced a qualitative study of the dynamics of the flat FRW model with three matter components one of which is barotropic and the other two are coupled, exchanging energy. The consideration of the flat model is motivated by observational data that indicate that the universe is flat [45][46][47]. Since these results also suggest that the density parameter of matter (including baryonic matter and cold dark matter) only accounts for Ω m ≃ 0.35, some additional form of dark energy should also intervene. Irrespective of the particular proportions taken by each of the major material components the mere fact that any other components might be present apart from the usual barotropic radiation and matter fluids of the standard cosmological model [48,49] raises the worry that the usual expansion rates during the radiation and Einstein-de Sitter epochs of the standard model would be significantly disrupted. Attracting scaling solutions provide a possible way of reconciling the presence of several matter components with the standard model as the various energy densities at stake all scale with the same rate as that of the barotropic fluid if the latter is non-vanishing. For the general type of interaction we considered in this paper, we have presented all the classes of parameter values for which system (52) exhibits non trivial scaling solutions. Those which correspond to equilibria outside the line x = 1, that is, for which ρ X = 0, are especially interesting from the physical point of view. These are the cases of Figures 3, 4.b) and c), 5.b), c), d) and e), 7.a) and 8. Among these, the cases of Figures 4.b), 5.b) and 8 are also remarkable. They correspond to situations where for a set of positive (or even full) measure of initial conditions in phase space the solutions are periodic, and so the relative abundance of each of the three fluids oscillates. We would also like to remark that this class of models opens the possibility of having simultaneously scaling and periods of inflationary behaviour.

A. Inflation
In fact, it is easy to see that the condition for inflationary behaviour (ä > 0) is and thus defines a half-plane in the x, y plane. This condition does not depend on the values of the parameters λ, α, β and δ, it just depends on γ X , γ 1 and γ 2 , and adjusting the latter the triangle of the physical phase space may have partial or full overlap with the inflationary half-plane. When the boundary line y * (x) = ǫx + (2 − 3γ X )/(3γ d ) lies in the y > 0 halfplane and intersects the x = 1 vertical line above the point (1, 1) all the trajectories are inflationary. This situation arises when γ X < 2/3 and γ 1 < 2/3 (recall that we are assuming γ 2 < γ 1 ). The opposite situation arises when the boundary line y * (x) lies in the y < 0 half-plane and intersects the x = 1 vertical line below the point (1, −1) in which case no trajectory is inflationary. This happens for γ X > 2/3 and γ 2 > 2/3. In the intermediate situations, namely when γ 2 < 2/3 < γ 1 , the inflationary region corresponds to the portion of the phase space triangle below the y * (x) line. The intersections of this line with the triangle are at the points Let us consider what happens in the scalar field models previously discussed in Section II. In the minimal coupling case the inflationary solutions lie below the y * (x) = (γ − 1)x + (2 − 3γ)/3 line. This line intersects the x = 1 frontier of the triangle at y = −1/3 regardless of its slope defined by the value of γ. For γ = 1, that is, for dust, the line is horizontal and the inflationary region is a triangle. As we consider smaller values of γ, the slope of the line y * (x) becomes increasingly negative and in the limit case of γ = 0 it becomes −1. Notice that this limit value of γ does not yield NTSS solutions in the interior of the phase-space domain. It corresponds to a cosmological constant and it is equivalent to having a nonvanishing vacuum energy in the exponential potential. The converse happens as we consider values of γ > 1. The slope increases up to a maximum value when the limit γ = 2 is chosen and the boundary line is parallel to y = x. So there is always an inflationary region in the phase diagram of the models and this means that portions of the trajectories approaching the attractors will exhibit inflationary transients. Let us consider now the question of whether the scaling solutions fall within those regions. In the ν 2 < 3γ < 6 case, the only fixed point lies on the x = 1 vertical line. When ν 2 < 2 its y coordinate satisfies y < −1/3 and hence falls within the inflationary region. In the ν 2 > 3γ case, the stable scaling solution x = 3γ/ν 2 ,ȳ = 3γ(γ − 1)/ν 2 belongs to the inflationary region only if γ < 2/3. This means that we may have NTSS which exhibit inflationary behaviour. However, as the most interesting models from the viewpoint of the late time behaviour of the universe are those for which the perfect fluid has γ ≥ 1, namely γ = 1, 4/3, the remarkable issue is that a non-negligible set of solutions naturally undergo a finite period of inflation before reaching the attractor.
The models leading to oscillatory behaviour are also interesting in what concerns inflation. Consider the model associated with the figure 8, that is, the model characterized by |ǫ| < 1, α = β = 1 and δ > 1 . The NTSS is located atx = 1/δ,ȳ = ǫ/δ and it corresponds to an expansion that tracks the non-interacting perfect fluid so that a ∝ t 2/(3γ X ) . Therefore it is immediate to see that this NTSS falls within the inflationary region of the phase-diagram if γ X < 2/3 and off it otherwise. Taking into consideration what was expounded at the beginning of this subsection, in order for the inflationary half-plane to overlap the phase-plane triangle one requires that γ 2 be smaller than 2/3. Assuming γ X > 2/3 and γ 2 < 2/3 we have then that the oscillatory trajectories beyond certain radius from the fixed point will undergo cyclic periods of inflation since they periodically cross the inflationary region of the phase-plane. Notice that since the model under consideration corresponds to a non-linear oscillator the periods associated with the trajectories increase from the immediate (and extremely small) neighbourhood of the fixed point, where the linear approximation holds true and T ≃ 2π/ǫ(1 − 1/δ) (1 − ǫ 2 )/(δ − 1), to the regions close to the limits of the triangle where the period becomes infinite.

B. Curvature scaling solutions
An interesting case which also emerges from our analysis regards the effect of a nonvanishing spatial curvature in models with two interacting fluids. Indeed taking the particular case where γ X = 2/3 for the non-interacting perfect fluid from eq. (37) we see that this is equivalent to having a term ρ ∝ a −2 in the Friedmann equation. Conversely, it is easy to verify that the usual curvature term taken as ρ k = −3k/a 2 satisfies Eqs. (35), (36), (37) with γ = 2/3. Thus the case of two coupled fluids in k = −1 models falls within the scope of our study.
From the definition of ǫ we have then ǫ = (2 − 3γ c )/(3γ d ) and so, according to our results of Section III, curvature scaling solutions, as defined in [13], exist when γ 2 < 2/3 < γ 1 . In this case the scale factor evolves as a ∝ t and both ρ 1 and ρ 2 self-adjust to the ρ 1 , ρ 2 ∝ a −2 behaviour of the curvature term. If both γ 1 and γ 2 are greater than 2/3 the origin is a global attractor and we recover the usual asymptotic behaviour found when we have two non-interacting perfect fluids. The curvature term eventually dominates and we have the vanishing of ρ 1 and of ρ 2 . This solution corresponds to well-known Milne universe.
In the limit case where ǫ = −1, and hence γ 2 = 2/3, the curvature term dominates in the future, but only a limited set of solutions corresponds to the depletion of both the ρ 1 and the ρ 2 components. Indeed almost all solutions end up in the x = −y border line and they correspond to solutions with a ∝ t which occur with the depletion of the ρ 1 component. In this case there are no inflationary solutions in the sense that the scale factor does not evolve with a power greater than 1 (we are in the limit coasting model). In the alternative limit case, when ǫ = 1 and hence γ 1 = 2/3, the curvature was the dominating component in the past and the cosmological models evolve either towards scaling solutions characterized by the depletion of the ρ 2 component or towards scaling solutions where the curvature vanishes. In this case all the solutions are inflationary.

C. Nucleosynthesis
It was pointed out in previous works on scaling solutions [4][5][6][7] that the most stringent bounds on the admissible densities of the components that are present in addition to the usual perfect fluid are set by the primordial nucleosynthesis of light elements [49]. Since the perfect fluid is radiation, the attractor solution is characterized by the usual of expansion a ∝ t 1/2 and so any deviations from the standard model yields of the light elements are a consequence of the number of degrees of freedom N(t nuc ) which are due to the extra matter components [6]. The limits that beset this number can be translated into a permitted range of energy density associated to the additional matter contributions which is Ω extra < ∼ 0.13 − 0.2 [6].

D. Decay of massive particles
We now briefly consider the possibility of applying the present model to a transient regime during the early universe when two particle species interact. For instance the decay of some massive particle species into a lighter one occurring out of equilibrium. This question has been analysed in the literature (for a review see Chapter 5 of [49] and references therein) and usually it is assumed that the massive particles decay into relativistic particles that rapidly thermalize so that, on one hand one may consider them as being a part of the radiation component, and on the other hand one does not have to consider reverse processes.
The model envisaged in this work enables one to relax the assumption of thermal equilibrium of the lighter species. Following [49] if we denote by ψ the decaying massive particle species, the relevant equations which are usually adopted arė where ρ ψ is the energy density of the ψ species, ρ r is the energy density of the radiation fluid which includes the thermalized daughter products of the decays, and the time at which the decays take place is given by τ ∼ t ∼ H −1 .
If we apply our model it becomes possible to consider the transient stage during which the lighter particle species energy density increases due to the decays of the ψ's and is not yet in thermal equilibrium with the radiation. Since there is no well-established thermokinetic prescription for this situation the purpose of this analysis is mainly illustrative and a more detailed investigation of the specifics of the process is left to a future work. From the viewpoint of our model the usual treatment found in the literature may be associated with taking ρ X as being the radiation perfect fluid, ρ ψ to be characterized by γ 2 = 1 and the daughter products to be described by ρ 1 with γ 1 = 4/3. Then we take τ ∼ ηH λ+1 ρ α 1 ρ β−1 2 in Eq. (59) where the parameters are left free. Thus we have a model characterized by ǫ = 1 and we immediately know that the dynamics corresponds to the phase-plane of the Fig. 3 c) where only the trajectories above the separatrix that connects the (0, 0) singular point to the saddle scaling solution on the x = 1 border of the phase diagram evolve towards equilibrium with the radiation component. These solutions correspond to a depletion of the ρ 1 component and, hence, of the decaying massive particles.
If we relax the assumption that the daughter products are in equilibrium with radiation and we keep the γ 1 parameter free, we may still have a non-trivial scaling solution in the interior of the phase-plane domain, for which both of the components interacting through the decays swiftly adjust themselves to track the perfect fluid behaviour of radiation. We would have then their thermalization and there is no depletion of the massive particle species in this case. For this to happen the lighter species must be characterized by a γ 1 > 4/3 in order to satisfy the condition |ǫ| < 1. According with Eq. (54) the combined density parameter settles at a value determined by γ 1 , since in the present case it completely defines ǫ (we have ǫ = (5 − 3γ 1 )/(3(γ 1 − 1)), and also by the values taken for α, β and δ (the latter is given by . However as presented in our study additional conditions must also be met in order to have an attracting NTSS. The relevant cases are those represented in Figs. 4 c), 5 c) and d) and 7 a) and a common condition δ > (1+ǫ (60) If one selects particular values of the remaining parameters this means yet another constraint on δ.
The main point to emphasize is then that there exists the possibility within the framework of a flat model of having a dynamical thermalization of the three components and, hence, the energy densities of the particle species adding to the radiation energy density. In the conventional models the relics of the ψ's do not contribute to the latter and are subject to the so-called Lee-Weinberg bound [49].