Dynamics of the cosmological models with perfect fluid in Einstein-Gauss-Bonnet gravity: low-dimensional case

In this paper we performed investigation of the spatially-flat cosmological models whose spatial section is product of three- ("our Universe") and extra-dimensional parts. The matter source chosen to be the perfect fluid which exists in the entire space. We described all physically sensible cases for the entire range of possible initial conditions and parameters as well as brought the connections with vacuum and $\Lambda$-term regimes described earlier. In the present paper we limit ourselves with $D=1, 2$ (number of extra dimensions). The results suggest that in $D=1$ there are no realistic compactification regimes while in $D=2$ there is if $\alpha>0$ (the Gauss-Bonnet coupling) and the equation of state $\omega<1/3$, the measure of the initial conditions leading to this regime is increasing with growth of $\omega$ and reaches its maximum at $\omega \to 1/3 - 0$. We also describe some pecularities of the model, distinct to the vacuum and $\Lambda$-term cases -- existence of the isotropic power-law regime, different role of the constant-volume solution and the presence of the maximal density for $D = 2$, $\alpha<0$ subcase and associated features.

in the Lagrangian which are quadratic with respect to the curvature; particular examples include R 2 and R µν R µν terms [6] in the Lagrangian of the Virasoro-Shapiro model [7,8], R µνλρ R µνλρ term [9] for the low-energy limit of the E 8 × E 8 heterotic superstring theory [10] and others. An important milestone was a discovery that there exist only one combination of the above terms which leads to ghost-free nontrivial gravitation interaction -the Gauss-Bonnet (GB) term [11] L GB = L 2 = R µνλρ R µνλρ − 4R µν R µν + R 2 .
This term was already known at that time, as it was first discovered by Lanczos [12,13] (and therefore sometimes it is referred to as the Lanczos term). It is an Euler topological invariant in (3+1) and lower dimensions while in (4+1) and higher it gives nontrivial contribution to the equations of motion. Soon after the idea of [11] was supported and extended in [14] on the case of higherthen-quadratic terms. So that the unified theory in the final form should include contributions from all possible curvature powers. In this regard the Lovelock gravity [15], which is the sum of all possible Euler topological invariants, which give nontrivial contribution to the equations of motion in particular number of dimensions, is of special interest. Of course, the Einstein-Gauss-Bonnet (EGB) gravity, which is under consideration in the current paper, is a particular case of more general Lovelock gravity.
Our everyday experience, as well as experiments in particle and space physics, clearly indicate that there are only three spatial dimensions (for example, in Newtonian gravity in more then three spatial dimensions there are no stable orbits while we are rotating around the Sun for ages).
So that to bring together the extra-dimensional theories and the experiment, we need to explain where are these extra dimensions. The commonly accepted answer is that the extra dimensions are compactified on a very small scale -similar to the "curling" of extra dimension into a circle in the original Kaluza-Klein paper. But this answer, in turn, gives rise to another one -how come that they became compact? The answer to this question is not so simple. The first attempts to answer this question involves solution known as "spontaneous compactification" [16,17]. Similar solutions but more relevant to cosmology were proposed in [17,18] (see also [19]). More natural way to achieve compactified extra dimensions is "dynamical compactification". The works on this direction involves different approaches [20,21] and different setups [22,23]. Also, apart from the cosmology, the studies of extra dimensions involve investigation of properties of black holes in Gauss-Bonnet [24][25][26][27][28][29][30][31] and Lovelock [32][33][34][35][36] gravities, features of gravitational collapse in these theories [37][38][39], general features of spherical-symmetric solutions [40], and many others.
As we shall see, the equations of motion for EGB and for more general Lovelock gravity are quite complicated and it is nontrivial to find exact solutions within them, so that one usually apply ansatz of some sort to obtain them. For cosmology the usual ansätzen are power-law and exponential -the former resemble Friedmann stage while the latter -accelerated expansion nowadays or inflationary stage in Early Universe. Power-law solution were studied in [17,41] and more recently in [19,[42][43][44][45], so that we can say that we have some understanding of their dynamics. One of the first considerations of the exponential solutions in the considered theories was done in [46], the recent works include [47], the separate description of the cases with variable [48] and constant [49] volume; see also [50] for the discussion of the link between existence of power-law and exponential solutions as well as for the discussion about the physical branches of the solutions. We have also described the general scheme for building all possible exponential solutions in arbitrary dimensions and with arbitrary Lovelock contributions taken into account [51]. Deeper investigation revealed that not all of the solutions found in [51] could be called "stable" [52]; see also [53] for more general approach to the stability of exponential solutions in EGB gravity.
The above approach gave us asymptotic power-law and exponential solutions, but not the evolution. Without full evolution we cannot decide if the asymptotic solution which we found realistic or not, could it be reached from vast enough initial conditions or is it just some artifact. To answer this question we need to go beyond exponential or power-law ansätzen and consider generic form of the scale factor. In this case, though, one would often resort to numerics. In [54] we considered the cosmological model with the spatial part being a product of three-and extra-dimensional spatially curved manifolds and numerically demonstrated the existence of the phenomenologically sensible regime when the curvature of the extra dimensions is negative and the EGB theory does not admit a maximally symmetric solution. In this case both the three-dimensional Hubble parameter and the extra-dimensional scale factor asymptotically tend to the constant values (so that threedimensional part expanding exponentially while the extra dimensions tends to a constant "size").
In [55] the study of this model was continued and it was demonstrated that the described above regime is the only realistic scenario in the considered model. Recent analysis of the same model [56] revealed that, with an additional constraint on couplings, Friedmann-type late-time behavior could be restored.
This investigation was performed numerically, but to find all possible regime for a particular model, the numerical methods are not the best practice. We have noted that if we consider the model with spatial section as the product of three-and extra-dimensional spatially flat subspaces, the equations of motion simplify and become first-order instead of the second one 1 . In this case the dynamics could be analytically described with use of phase portraits and so we performed this analysis. For vacuum EGB case it was done in [57] and reanalyzed in [58]. The results suggest that in the vacuum model has two physically viable regimes -first of them is the smooth transition from high-energy GB Kasner to low-energy GR Kasner. This regime exists for α > 0 (Gauss-Bonnet coupling) at D = 1, 2 (the number of extra dimensions) and for α < 0 at D 2 (so that at D = 2 it appears for both signs of α). Another viable regime is the smooth transition from high-energy GB Kasner to anisotropic exponential regime with expanding three-dimensional section ("our Universe") and contracting extra dimensions; this regime occurs only for α > 0 and at D 2. Apart from the EGB, similar analysis but for cubic Lovelock contribution taken into account was performed in [59,60].
The similar analysis for EGB model with Λ-term was performed in [61,62] and reanalyzed in [58]. The results suggest that the only realistic scenario is the transition from high-energy GB Kasner to anisotropic exponential solution, it requires D 2, see [58,61,62] for exact limits on (α, Λ). The low-energy GR Kasner is forbidden in the presence of the Λ-term so the corresponding transition do not occur.
In the studies described above we have made two important assumptions -we considered both subspaces being isotropic and spatially flat. But what will happens in we lift these conditions?
Indeed, the spatial section being a product of two isotropic spatially-flat subspaces could hardly be called "natural", so that we considered the effects of anisotropy and spatial curvature in [63].
The effect of the initial anisotropy could be demonstrated on the example of vacuum (4 + 1)dimensional EGB cosmology with Bianchi-I-type metric (all the directions are independent), where the only future asymptote is nonstandard singularity [43]. Our analysis [63] suggest that the transition from GB Kasner regime to anisotropic exponential expansion (with expanding three and contracting extra dimensions) is stable with respect to breaking the symmetry within both threeand extra-dimensional subspaces. However, the details of the dynamics in D = 2 and D 3 are different -in the latter there exist anisotropic exponential solutions with "wrong" spatial splitting and all of them are accessible from generic initial conditions. For instance, in (6 + 1)-dimensional space-time there are anisotropic exponential solutions with [3 + 3] and [4 + 2] spatial splittings, and some of the initial conditions in the vicinity of E 3+3 actually end up in E 4+2 -the exponential solution with four and two isotropic subspaces. In other words, generic initial conditions could easily end up with "wrong" compactification, giving "wrong" number of expanding spatial dimensions (see [63] for details).
The effect of the spatial curvature on the cosmological dynamics could be dramatic -say, positive curvature changes the inflationary asymptotic [64,65]. In EGB gravity the influence of the spatial curvature (also described in [63]) reveal itself only if the curvature of the extra dimensions is negative and D 3 -in that case there exist "geometric frustration" regime, described in [54,55] and further investigated in [56].
The current manuscript is, from one hand, a direct continuation of our ongoing study of all possible cosmologically sensible cases in EGB gravity: vacuum models (see [57]) and models with Λ-term (see [61,62]; both cases with revised presentation of the results are reviewed in [58]). Both of them are important milestones on our way to understand the cosmological dynamics of EGB (and more general Lovelock) gravity, but if we want to describe our Universe (and this is the goal of any realistic physical theory regardless of the field), then we must make the theory as close as possible to the reality. And in the reality we observe not vacuum and not only Λ-term, but ordinary matter as well. So that in current paper we started to consider another important matter source for EGB gravity -matter in form of perfect fluid -making the study of the entire EGB case more complete. It worth mentioning that we already performed consideration of some particular GB and EGB models with matter source in form of perfect fluid: some initial considerations in [47], some deeper studies of (4+1)-dimensional Bianchi-I case in [43] and deeper investigation of power-law regimes in pure GB gravity in [45]. In the last cited paper we obtained the results on which we will rely in the current study and which we are going to use to explain some of the features. The most important of these results is the fact that ω = 1/3 is the critical value for the equation of state which interchange the dominance of GB and matter terms: for ω < 1/3 past asymptote is GB-dominated and future asymptote is matter-dominated, for ω > 1/3 the situation changespast asymptote is matter-dominated while future asymptote is GB-dominated. But there is one important issue which disallow us from directly implementing these results to our current study -now we have EGB gravity, so that the low-energy regime now is described by Einstein-Hilbert Lagrangian and for it the limiting equation of state is ω = 1 (Jacobs solution [66]). And since we do not consider ω > 1, all of our low-energy regimes are the same so the difference lies only in high-energy ones. So that the analogue with GB case is not direct and there could be other interesting features and observations.
The manuscript is structured as follows: in the following Section II we write down equations of motion, then consider the dynamics for D = 1 in Section III and for D = 2 in Section IV, discuss the obtained results and the interesting observations in Section V and finally draw conclusions in Section VI.

II. EQUATIONS OF MOTION
The structure of the general Lovelock gravity [15] is the following: its Lagrangian consists of where δ i 1 i 2 ...i 2n j 1 j 2 ...j 2n is the generalized Kronecker delta of the order 2n. One can note that L n is the Euler invariant in D < 2n spatial dimensions and so it does not give nontrivial contribution into the equations of motion. So the Lagrangian density for any given spatial dimensions D is sum of all Lovelock invariants (1) up to n = D 2 which give nontrivial contributions into equations of motion: where g is the determinant of metric tensor and c n are coupling constants of the order of Planck length in 2n dimensions and summation over all n in consideration is assumed. We consider the metric in the form As we mentioned earlier, we are interested in the dynamics in Einstein-Gauss-Bonnet gravity, which is second order Lovelock gravity and so we consider n up to two (n = 0 is the boundary term while n = 1 is Einstein-Hilbert and n = 2 is Gauss-Bonnet contributions). In this paper we studying the cosmological dynamics in the presence of the matter source in form of perfect fluid with an equation of state p = ωρ, so the standard notation for its stress-energy tensor is assumed: Substituting metric (3) into the Lagrangian and following the usual procedure with use of (4) gives us the equations of motion: the dynamical equation that corresponds to h, and the constraint equation, Looking at (7)-(9) one can notice that the structure of the equations depends on the number of extra dimensions D (terms with (D − 2) multiplier nullifies in D = 2 and so on), so that we shall consider all these distinct cases separately. Formally, the system should be complimented by the continuity equatioṅ but the resulting system (7)-(10) is overdetermined, so for practical purposes we skip (10) and use (7)-(9) for analysis.
In the previous papers, dedicated to either vacuum or Λ-term cases in EGB or cubic Lovelock gravity [57][58][59][60][61][62], we solved the constraint equation (9) with respect to either H or h, resulting in one to three different branches of the solution with several parameters. Then the obtained expressions are substituted into (7)-(8) and the system was solved with respect toḢ andḣ. As a result, we have bothḢ andḣ as a single-variable function with several parameters and investigated its phase portrait to find all the asymptotic regimes, stable points etc. In the current paper, though, this method will not work, because with the same number of independent equations we have one more variable -the density ρ. So for practical purposes we use the following procedure: we substitute density ρ from (9) into (7) and (8) with use of the equation of state p = ωρ; the resulting system is first-order with two equations and two variables and so we can plot directional vector phase portrait to visualize the dynamics. The asymptotes as well as exact solutions are derived analytically from In the course of the paper we use the results we previously obtained in [45]. In particular, we demonstrated that: for ω < 1/3 past GB asymptote is preserved while future asymptote is destroyed by matter and is replaced by the corresponding matter-dominated regime. For ω > 1/3 the situation is opposite -the past asymptote is destroyed and replaced by matter-dominating regime while the future asymptote remains GB-dominated.
Before proceeding to the particular cases, let us introduce the notations we are going to use through the paper. We denote Kasner regime as K i where i is the total expansion rate in terms of the Kasner exponents p i = (2n − 1) where n is the corresponding order of the Lovelock contribution (see, e.g., [19]). So that for Einstein-Hilbert contribution n = 1 and p i = 1 (see [67]) and the corresponding regime is K 1 , which is usual low-energy regime in vacuum EGB case (see [57,58]).
For Gauss-Bonnet n = 2 and so p i = 3 and the regime K 3 is typical high-energy regime for EGB case (again, see [57,58,61,62]). For cubic Lovelock n = 3 and so p i = 5 and the regime K 5 is one of the high-energy regimes in that case [59,60].
Another power-law regime is what is called "generalized Taub" (see [68] for the original solution).
We mistakenly taken it for K 3 in [57], but then in [58] corrected ourselves and explained the details.
It is a situation when for one of the subspaces the Kasner exponent p is equal to zero and for another -to unity. So we denote P (1,0) the case with p H = 1, p h = 0 and P (0,1) the case with p H = 0, p h = 1.
So that in EGB for P (1,0) we have The exponential solutions are denoted as E with subindex indicating its details -E iso is isotropic exponential solution and E 3+D is anisotropic -with different Hubble parameters corresponding to three-and extra-dimensional subspaces. But in practice, in each particular case there are several different anisotropic exponential solutions, so that instead of using E 3+D we use E i where i counts the number of the exponential solution (E 1 , E 2 etc). In case if there are several isotropic exponential solutions, we count them with upper index: E 1 iso , E 2 iso etc. And last but not least regime is what we call "nonstandard singularity" and denote it is as nS.
It is the situation which arise in Lovelock gravity due to its nonlinear nature. Since the equations (7)-(8) are nonlinear with respect to the highest derivative (Ḣ andḣ in our case), when we solve them, the resulting expressions are ratios with polynomials in both numerator and denominator.
So there exist a situation when the denominator is equal to zero for finite values of H andor h.
This situation is singular, as the curvature invariants diverge, but it happening for finite values of H andor h. Tipler [69] call this kind of singularity as "weak" while Kitaura and Wheeler [70,71] -as "type II". Our previous research demonstrate that this kind of singularity is widely spread in EGB cosmology -in particular, in totally anisotropic (Bianchi-I-type) (4 + 1)-dimensional vacuum cosmological model it is the only future asymptote [43].
In this case the equations of motion (7)-(9) take form (H-equation, h-equation, and constraint correspondingly): As we mentioned, the procedure is different from our previous papers so we are going to explain it in detail on the example of D = 1 case. We substitute ρ from (13) into (11) and (12) with use of the equation of state p = ωρ and solve (11)-(12) with respect toḢ andḣ: One can see that in this case we have only one nonstandard singularity at H = ±1/(2 √ −α)) and for negative α only. If we solve (14), we obtain the locations of the exponential solutions, and they are: only for negative α, and so-called "constant-volume solution" (CVS) (see [49]) with For the CVS to exist the radicand should be positive, which happening if either {α > 0, ω > 1/3} or {α < 0, ω < 1/3}. If we substitute the CVS conditions to (13) we shall see that only {α < 0, ω < 1/3} provide positive energy density.
Now with the preliminaries done, let us draw the phase portrait of the system (14). Despite the fact that, according to [45], the fundamental change occurs at ω = 1/3, we notice that at ω = 0 there are also some minor changes which worth mentioning. So that the situation for α > 0 is presented in Fig. 1. There we presented ω < 0 case on (a) panel, 0 < ω < 1/3 on (b) panel and ω > 1/3 on (c) panel. Green line corresponds toḣ = 0 while black toḢ = 0 (see (14)). Dark blue area corresponds to unphysical (H, h) pairs with ρ < 0.
Let us start with ω < 0 case on (a) panel and let us focus on H > 0 half-plane -the situation for H < 0 is obviously time-reversal with respect to H > 0 -we will comment it a bit later. So that for H > 0 one can see that there are two past asymptotes -P (1,0) and P (0,1) , depending on the overall relation between H and h: for H h it is P (1,0) while for H h it is P (0,1) . As of the future asymptote, it is P iso -isotropic power-law expansion. Let us note that the future asymptote is matter-dominated while past -GB-dominated, all well according to [45]. For H < 0, as we mentioned, the situation is time-reversed -so that P iso is past asymptote while P (1,0) and P (0,1) are future asymptote. Let us note that P iso is stable while being future asymptote but is unstable as a past -even a slightest deviation from isotropy turn the evolution towards either P (1,0) or P (0,1) .
So that in this case there are only two cosmological regimes -P (1,0) → P iso and P (0,1) → P iso and neither of them has realistic compactification.
The next case to consider is 0 < ω < 1/3 and it is presented in Fig. 1(b). One can immediately see the changes from ω < 0 case ( Fig. 1(a)) -theḣ = 0 andḢ = 0 curves have different structure and it slightly affects "zeroth" direction for P (1,0) and P (0,1) asymptotes -in ω < 0 case it is p → 0 − 0 (one can see that p → 0 is reached from below) while for 0 < ω < 1/3 it is p → 0 + 0 (so that it is reached from above). Except for this minor feature there are no other differences between this and the previous cases -the regimes are the same, so we just refer to the previous case for the description. The last case for α > 0 is ω > 1/3 and it is presented in Fig. 1(c). One can immediately see the difference with both previous cases -the past asymptote switched to P iso , and this behavior is in agreement with [45]. The future asymptote is the same as in the previous cases -P iso , so that the only cosmological regime in this case is P iso → P iso -anisotropic transition from unstable isotropic power-law regime to a stable one. Let us note that these are different isotropic power-law solutions -the past asymptote is matter-dominated GB while future -matter-dominated GR so that they have different p i and so different expansion rates. This regime is unique for the matter casein our previous studies dedicated to vacuum and Λ-term cases we never encountered regime of this type.
One can also notice exponential CVS in the ρ < 0 domain where green (ḣ = 0) and black (Ḣ = 0) lines intersect. This behavior is in exact agreement with described above and we expect CVS with positive energy density in {α < 0, ω < 1/3} case.
One cannot miss that the past asymptote changed at ω = 1/3 while future asymptote remains the same. The reason behind it is simple -the future asymptote in this case is governed by GR and for GR the qualitative change similar to that is happening at ω = 1 (Jacobs solution [66]).
And since we limit physical equation of state within ω ∈ (−1, 1], the entire allowed range is within ω 1 and so has the same regime P iso . Now let us turn to α < 0 cases. They have richer dynamics and contain, as we demonstrated earlier, nonstandard singularities and exponential solutions. Due to this richness we decided to present them in different figures. We start with α < 0, ω < 0 in it we can see that there is a fine structure for high h within H < H nS as well as in the vicinity of E iso -we presented these areas in detail on panels (b) and (c) respectively. Apart from them, we detect nS → P iso from h < 0, H < H nS area. Majority of the h > 0, H < H nS trajectories are P (0,1) → P iso , but at high h there are P (0,1) → nS transitions, we denoted the corresponding area to blue in Fig. 2(b). So that within H < H nS there are only three regimes: nS → P iso , P (0,1) → nS and P (0,1) → P iso .
The H > H nS area is shown in detail in Fig. 2(c). We can see that at high h the regime is nS → E iso , while the situation at low and negative h is more complicated. At high H the past asymptote is P (1,0) and it gives rise to two regimes -P (1,0) → E iso and P (1,0) → nS, the latter is colored in blue on panel (c). The final regime is another nS → E iso , which originate from the same nS but at h → −∞. The boundary between P (1,0) → E iso and nS → E iso is denoted by light blue line in Fig. 2(c).
Before moving on to the next case, it is important to mention the CVS solution -we theoretically predicted its existence but it has not participated in the regimes description. In reality it exists at h < 0, H > H nS where green (ḣ = 0) and black (Ḣ = 0) lines intersect. For ω < 0 the intersection angle is very small so on a chosen scales it is very hard to locate it -in the next case the situation improves and we are going to describe CVS in detail.
Now let us move on to α < 0, 1/3 > ω > 0 which we presented in Fig. 3(a) So that the regimes are (see Fig. 3(b)): nS → E iso (formally two different, as there are two nS at h → ±∞), P (1,0) → E iso and P (1,0) → nS (the latter is represented as blue area) and nS → nS as a small part of blue area, roughly separated by black line. One can see that of all regimes only P (0,1) → P iso and P (1,0) → E iso are nonsingular, but neither of them have realistic compactification.
In this case we can clearly see CVS, it is represented as E CV S . One can also see that despite existent, it is neither past nor future asymptote for any regime. From Fig. 3(c) one can clearly see why is it -it is a pole on the phase portrait, referring to the marginal stability of the solution.
It is expectant -while studying the stability of the exponential solutions in [52] (and later it was generalized in [53]), the stability analysis works for H i = 0. On the contrary, CVS has, as its name assumes, H i = 0, making the standard stability analysis unviable, which could be treated as a marginal stability, and now we exactly demonstrate it. For vacuum and Λ-term cases in EGB as well as in the cubic Lovelock gravity we also have CVS, but due to the lesser number of the degrees of freedom, CVS in these cases have directional stability (i.e. stable from H i > 0 direction and unstable from H i < 0 direction, see [58][59][60][61][62]).
Finally let us consider α < 0, ω > 1/3 case, which is presented in Fig. 3(d). The dynamics in this case is much more simple then in previous cases -for both H < H nS the only regime is nS → P iso for both h > 0 and h < 0 while for H > H nS the only regime is nS → E iso for both h > 0 and h < 0 So that there are neither nonsingular nor realistic regimes in this case.
To conclude D = 1 case with matter in form of perfect fluid, we have three nonsingular regimes for α > 0: P (1,0) → P iso , P (0,1) → P iso and P iso → P iso , but neither of them have realistic compactification. Let us note that there are no other regimes apart from the just mentioned, so that all possible regimes in α > 0 case are nonsingular. The situation with α < 0 is more complicated and we have P (0,1) → P iso and P (1,0) → E iso as nonsingular regimes; they exist for ω < 1/3 and for ω > 1/3 there are no nonsingular regimes. Nevertheless, neither of nonsingular regimes for α < 0 have realistic compactification either, making D = 1 case degenerative with respect to realistic compactifications.

IV. D = 2
In this case the equations of motion (7)-(9) take form (H-equation, h-equation, and constraint correspondingly): 6Ḣ + 12H 2 + 2ḣ + 2h 2 + 6Hh + 8α 3(Ḣ + H 2 )(H 2 + 2Hh)+ Similar to the previously described procedure, we substitute ρ from (17) into (15) and (16), use the equation of state p = ωρ and solve (15)-(16) with respect toḢ andḣ: The positions of the nonstandard singularities are determined by zeros of Q and one can see that now, unlike D = 1 case, they are not just straight lines. One can also see that Q is quadratic with respect to h so one can easily solve and plot it. Calculating the discriminant of Q with respect to h gives us From it one can easily see that for α > 0 discriminant is positive if H 2 < 1/(4α) while for negative α the discriminant is always negative -so that for negative α there are no nonstandard singularities -very unexpected result, considering that for D = 1 nonstandard singularities exist only for negative α.
The exponential solutions are defined by P 1 = P 2 ≡ 0 and are as follows: trivial solution , which exist only for negative α, CVS which we describe below and one more solution with and from the requirement of the positivity it exists in two domains: (α > 0, ω < 1/3) and (α < 0, ω > 1/3). Substituting the solution into the integral of energy gives us so that only (α > 0, ω < 1/3) domain gives CVS with positive energy density.
With the preliminaries done, let us present the resulting phase portraits. From the structure of the nonstandard singularities and exponential solutions one can see that D = 2 case is more complicated than the previous D = 1 and so the structure of the regimes. We presented α > 0 cases in Figs. 4-6.
The first case to consider is α > 0, ω < 0, presented in Fig. 4. On (a) panel we presented general view of the phase portrait. One can see the difference from D = 1 cases -now we have two branches and they are separated by ρ < 0 unphysical region, colored in dark blue. As the dynamics is more complicated, it is natural to focus on different regions of the map and mark different features on them. This way we focus on h > 0 half-plane on (b) panel. One can immediately see the difference with D = 1 case -in the latter we have P (0,1) as a past asymptote at h → +∞, H → 0 while now in D = 2 it is replaced by nonstandard singularity. So that the regimes are nS → P iso for h H and P (1,0) → P iso for h H; the latter is the same as in D = 1 case.
On (c) panel we presented the phase portrait focused on h < 0 half-plane. Despite the fact that the very "upper" part (above the ρ < 0 separation region) formally has h < 0, dynamically it belongs to h > 0, so we omit it from the description here. The h < 0 half-plane also has "double" nonstandard singularity, as well as K 3 as high-energy past asymptote. The dynamics in the vicinity has the same property as in D = 1 -it is a pole and so is not participating in the dynamical evolution. So that the regimes are (see Fig. 4(d)): K 3 → E 3+2 and K 3 → nS, separated by dashed blue line, both originate from K 3 -in D = 1 there is no such high-energy regime. The remaining two regimes originate from nS: nS → nS (located between two nS) and nS → E 3+2 . One can see that of these regimes K 3 → E 3+2 is both nonsingular and has realistic compactification. But its area of emergence is quite small with respect to the total initial conditions area.
So that in α > 0, ω < 0 case we have a number of regimes with two of them -P (1,0) → P iso and K 3 → E 3+2 being nonsingular and the latter in addition has realistic compactification.
The next case is α > 0, 1/3 > ω > 0, presented in Fig. 5. The general panels layout is the same as in Fig  To conclude, in α > 0, 1/3 > ω > 0 case we have the same regimes as in ω < 0, but the measure if the physically realistic K 3 → E 3+2 regime is growing with increase of ω until it reach its maximum at ω → 1/3 − 0.
The final case for α > 0 is ω > 1/3, presented in Fig. 6. Similar to the previous cases, the panels layout is as follows: on (a) panel we presented general view of the phase portrait, focus on h > 0 is done on (b) panel and focus on h < 0 is on (c) panel; unlike previous cases there is no "finestructure" of the regime in the vicinity of E 3+2 . Similar to D = 1, at ω = 1/3 we have a change of past asymptote -now it is P iso , so that on h > 0 half-plane the main regime is P iso → P isoanisotropic transition between two different isotropic power-law regimes -the same regime we have in D = 1, ω > 1/3 case. In the h < 0 half-plane we also have change of past asymptote -now it is solely nonstandard singularity, so that the only regimes for h < 0 are nS → nS and nS → E 3+2 (see Fig. 6). To conclude D = 2 α > 0 regimes, we have two nonsingular regimes P (1,0) → P iso and K 3 → E 3+2 at ω < 1/3 with the latter having realistic compactification; its measure over the total initial conditions space is growing with the increase of ω until reaching its maxima at ω → 1/3 − 0. For ω > 1/3 the only non-singular regime is P iso → P iso . Now let us turn our attention to the α < 0 regimes, which are a bit more difficult to attend.
Let us have a look on (17) and calculate its discriminant with respect to h: where x = H 2 . One can check that for α > 0 (and ρ 0 of course) D h > 0 always, so that ∀ρ > 0 ∃ h 1,2 which are solutions of (17). On the contrary, for α < 0 the situation is different: ∀ H ∃ ρ 0 such as for ρ > ρ 0 we have D h < 0. This means that for α < 0 there exists maximal value for density we presented the plot of this relation in Fig. 7(a). From the graph and from the expression one can see that for H 2 < −1/(36α) the maximal density formally is negative, but (20) nevertheless is positive, so that the limit is absent. For comparison reasons we presented α > 0 case in Fig. 7(b)one can easily see that ρ is unlimited there, as there are neither upper limit for h 1 nor lower limit for h 2 . The same situation is for α < 0 and H 2 < −1/(36α), presented in Fig. 7 asymptotes could be seen from all panels of Fig. 8 -they are isotropic exponential solution E iso and isotropic power-law expansion P iso . The separation of different regimes over (H, h) plane could be seen in Fig. 8(c) -they are P (0,1) → E iso , P (0,1) → P iso (coming from h → +∞ on H > 0) and One can note in Fig. 8(c), that apart from two described crossings ofḣ = 0 (green curve) withḢ = 0 (black curve) -P iso and E iso -there is one more more inbetween. As we mentioned earlier, these crossing give rise to exponential solutions -P iso , located at (0, 0) is trivial solution while isotropic exponential solution is nontrivial one. So that formally the third crossing also could be obtained under certain manipulation (it is located at exponential solution in the general sense. One can also note that it happening at ρ = ρ 0 which makes it unstable, so that it closer to CVS from α > 0 then to "normal" exponential solutions. The situation with α < 0, ω > 0 is not much different from the described above ω < 0. We presented the corresponding phase portraits in Fig. 9: 1/3 > ω > 0 case with focus on h > 0 we Dashed blue line corresponds to the ρ 0 from (21). Dark blue area corresponds to the unphysical ρ < 0 initial conditions. Green curve points location ofḣ = 0 while black -Ḣ = 0 (see the text for more details).
presented on (a) panel, ω > 1/3 case with focus on h > 0 we presented on (b) panel and ω > 1/3 case with focus on h < 0 we presented on (c) panel. Comparing Figs. 9(a, b) with Fig. 8(a) shows us little difference and these differences do not change the regimes; the same could be stated about comparison of Fig. 9(c) with Fig. 8(b); the situation in the central region -between P iso and E iso big difference lies in the presence of CVS in Fig. 9(c) -but, as it happening in ρ < 0 area (as we found above), it has no impact on the physical regimes.
One can also see that only one of them, K 3 → E 3+2 , has realistic compactification. It exists for α > 0, ω < 1/3 and the measure of the initial conditions leading to this regime is increasing with growth of ω, reaching its maximum at ω = 1/3 − 0.

V. DISCUSSIONS
Let us summarize the results obtained through the paper. We have considered D = 1, 2 (the number of extra dimensions) flat cosmological models with the spatial section to be a product of three-and extra-dimensional subspaces. As a source we consider perfect fluid, which exist in the entire space (not only in the three-dimensional subspace). We derived analytically the locations of the exponential solutions and nonstandard singularities and plot phase portraits of the evolutionary trajectories to find all possible regimes for the entire range of initial conditions and parameters.
The results for D = 1 demonstrate that for α > 0 the non-singular regimes are P (1,0) → P iso , P (0,1) → P iso and P iso → P iso , but neither of them have realistic compactification; there are no nonstandard singularities for α > 0. For α < 0 situation changes -there are nonstandard singularities and so singular regimes emerge. Nonsingular regimes for α < 0 include P (0,1) → P iso and P (1,0) → E iso and both of them exist only for ω < 1/3; for ω > 1/3 all regimes are singular (have nonstandard singularity as either past or future asymptote or both). But despite the fact that there is a number of non-singular regimes, neither of them has realistic compactification, making D = 1 case degenerative in this sense.
The results for D = 2 somewhat opposite to D = 1: for instance, for D = 1 nonstandard singularities exist only for α < 0 while for D = 2 they exist only for α > 0. Also, in D = 1 case nonstandard singularities are located on a fixed positions (straight lines) while in D = 2 they located along nontrivial curves. The α > 0 cases has a lot of singular regimes, but there are three nonsingular: P (1,0) → P iso and K 3 → E 3+2 for ω < 1/3 and P iso → P iso for ω > 1/3. Among them K 3 → E 3+2 has realistic compactification -the only regime with realistic compactification discovered through the course of this paper. For α < 0 all the regimes are nonsingular and they are the same for both ω > 1/3 and ω < 1/3: P (0,1) → E iso , P (0,1) → P iso , P (1,0) → E iso and K 3 → P iso ; one can see that neither of them has realistic compactification.
In the course of the investigation we have faced a number of observations and issues which require additional discussion and we are going to do it now. First of all, one can note that in D = 1 and α > 0, D = 2 there is a change of regimes while we cross ω = 1/3, and mainly this involve past asymptote regimes. The reason behind it was several times mentioned in the course of the paper and now we want to explain it in detail. In [45] we have demonstrated that ω = 1/3 is a critical value for the equation of state which separate two dynamically distinct cases in (pure) GB cosmology in the presence of the matter in form of perfect fluid -the cases of GB dominance and matter dominance. In particular, for ω < 1/3 the past asymptote (high-energy) is GB-dominated while the future asymptote (low-energy) changes to matter-dominated. For ω > 1/3 the situation is opposite -the past asymptote (high-energy) is matter-dominated while the future asymptote (low-energy) is GB-dominated. In [45] we clearly demonstrated that for ω > 1/3 initial singularity is isotropic one (P iso ) while for ω < 1/3 it is standard K 3 . The same way, for ω > 1/3 late-time regime is K 3 while for ω < 1/3 is isotropic power-law regime P iso . But this situation cannot be directly applied to our case -we have EGB gravity so that the low-energy regime is governed by GR -in vacuum the regime should be K 1 but in the presence of matter it is P iso . And these two P iso are actually different -the former of the mentioned is P GB iso and it is governed by GB in the presence of matter while the latter is P GR iso -isotropic power-law expansion governed by GR in the presence of matter. Also, for GR, as we mentioned in [45], the critical equation of state in ω = 1 which corresponds to Jacobs solution [66], and since we do not consider ω > 1, for all ω under consideration the late-time asymptote (if nonsingular!) is P GR iso . From this explanation one can understand why only past (high-energy) asymptote is affected by the change of ω ≷ 1/3. And well according to the results from [45], for ω > 1/3 the past asymptote becomes P iso in D = 1 and α > 0, D = 2 cases. However, for α < 0, D = 2 it is not the same and the next observation clarify the reason behind it.
The mentioned α < 0, D = 2 case has another interesting feature -the maximal value for the density for a open range of the initial conditions (H 2 > −1/(36α)). In all the other cases -entire D = 1 as well as α > 0, D = 2 -the value for the density is unbounded. In a sense, the situation is similar to the Λ-term cases [58,61,62], when the solution does not exist in a certain region of (H, Λ) for a given α. The reason both in the Λ-term cases and here is the same -certain quadratic equation has negative discriminant in this certain region. But in our case, as now we have more degrees of freedom, the situation is more complicated, and so the analogue is not direct.
Still, the maximal density exist for all H 2 > −1/(36α), reaching ρ 0 → ∞ at H → ±∞, so that formally at H → ∞ the maximal density limit is lifted. But in reality it is not -indeed, any small violation from H = ∞ gives rise to ρ 0 < ∞ and so the matter-dominated regime would be destroyed. That is the reason we do not see it -due to limit on the maximal density it cannot be reached in EGB gravity. Still, in pure GB it formally exists (see [45]).
In the course of study we detected what we call constant-volume solution (CVS) -anisotropic exponential regime whose dimensions are expanding in a way so that the volume of the space is kept constant. In [49] we investigated this regime in great detail for the totally anisotropic flat (Bianchi-I-type) metrics. In the course of study we encountered CVS in Λ-term EGB case (see [58,61,62]), and there this regime is accessible (so that it could be initial or final asymptote for a trajectory) and has directional stability (see [62] for details). Unlikely the Λ-term case, in our current investigation CVS is inaccessible -it is a pole-like singular point on the phase space (see e.g. Fig. 3(c)); though, the solution itself formally exist and we derived it for each particular case. The reason behind this difference lies in different number of independent variables between vacuum and Λ-term cases (one independent variable) on one hand and the case with matter in form of perfect fluid (two independent variables) on the other hand. It seems that this additional bound in vacuum and Λ-term cases plays the crucial role, but formal in-depth investigation of this feature is required and we are going to perform it shortly.
Finally we can compare our results with those obtained in EGB vacuum and Λ-term cases.
Indeed, formally assuming ρ → 0 we should obtain the regimes from the vacuum case, but there could be deviations, considering possible changes as ω ≷ 1/3. Also, for ω → −1 (assuming it is not singular!) we expect to recover Λ-term regimes, but only for Λ > 0, as with ρ > 0 we cannot recover AdS regimes. To perform the comparison, we use our previously obtained results for vacuum [57] and Λ-term [61,62] regimes as well as their revision in [58].
So in D = 1 the vacuum regimes are [57,58] P (1,0) → K 1 for α > 0 and P (1,0) → E iso , nS → E iso and nS → K 1 for α < 0, starting from large H. Our results for D = 1, α > 0 are presented in Fig. 1 and one can see that ρ → 0 and ω < 1/3 regimes would be P (1,0) → P iso -the same past asymptote as in the vacuum case but different future asymptote. Indeed, for exact ρ ≡ 0 we would have K 1 but even tiny nonzero amount of matter destroy this regime turning it into P iso eventually. above, different from the low-energy P iso ) -neither of the regimes are the same as in the vacuum case. For α < 0 and ω < 1/3 the situation is the same -high-energy regimes in our case are the same as in the vacuum case while the low-energy regime is replaced with P iso (see Figs. 2 and 3(a,   b)). For ω > 1/3 the situation is also similar to that at α > 0 -neither of the regimes are the same as in the vacuum case. So that the comparison of the vacuum and perfect fluid regimes for D = 1 reveals that for ω < 1/3 the past asymptote is the same for both cases while the future asymptote is different; for ω > 1/3 both asymptotes are different.
The structure of the regimes for D = 1 Λ-term case is more complicated [58,61]. For α > 0, Λ > 0 the regimes are P (1,0) → E iso and P (0,1) → E iso ; the corresponding regimes could be derived from Fig. 1(a). For ω → −1 the density would become constant and so P iso could not be obtained and would be replaced with E iso . This way, only past asymptotes are the same in both Λ-term and perfect fluid cases -similar to what we had previously from the comparison with vacuum cases.
For α < 0, Λ > 0 the situation is more complicated, as there are two subcases -αΛ ≷ −3/2. For αΛ < −3/2 the Λ-term regimes are P (0,1) → nS and P (1,0) → nS and they both could be obtained in the perfect fluid case (see Fig. 2). For αΛ > −3/2 the Λ-term regimes are P (1,0) → E 1 iso , nS → E 1 iso , nS → E 2 iso and P (0,1) → E 2 iso . First two of them could easily be retrieved (see Fig. 2) but not the remaining two. The situation with them resemble the previous case when P iso was replaced with isotropic exponential solution, so now it is the same but the isotropic exponential regime is different from the used in the first pair. So that for αΛ < −3/2 both regimes are retrieved completely while for αΛ > −3/2 we retrieve all but low-energy regimes -the same as for α > 0.
The cases with D = 2 have more complicated structure; vacuum regimes [57,58] for α > 0 are P (1,0) → K 1 for one and K 3 → E 3+2 ← nS → nS ← K 1 for another branch. Comparing these regimes with those obtained with use of ρ → 0 from Fig. 4(c) we can verify that all except the low-energy regime are the same; the low-energy regime K 1 is replaced with P iso for the same reasons as discussed in the D = 1 case. The regimes in Fig. 5 are also the same but not in Fig. 6, so that we recover all regimes except low-energy one for ω < 1/3 and recover neither for ω > 1/3 -again, in agreement with D = 1 results.
The resulting regimes for α < 0 D = 2 vacuum model are K 3 → K 1 for one and P (1,0) → E iso ← nS with K 1 → nS for another branch. Similarly to α > 0 case, some of the regimes could be easily retrieved: K 3 → K 1 is "transformed" into K 3 → P iso ; P (1,0) → E iso is the same in both. Two remaining regimes are more difficult to link: K 1 is replaced with P iso , but now there is change of direction -the vacuum regime K 1 → nS is replaced with P (0,1) → P iso -as we seen before, K 1 is replaced with P iso , but here the direction of evolution also has been changed; nS is replaced with P (0,1) . The same change happened with the remaining regime: nS → E iso is replaced with P (0,1) → E iso . So that, part of D = 2 vacuum regimes also could be obtained from perfect fluid approach, but with some changes, like mentioned K 1 → nS regime.
Finally let us connect D = 2 Λ-term regimes [58,61] with those obtained from the perfect fluid approach in the current paper. For Λ > 0, α > 0 regimes for D = 2 have "fine structure" [61] -different sequence of regimes depending on αΛ ∈ [15/32; 1/2]. We generally do not recover this fine structure; and the recovered regimes along ρ = const ≡ Λ do not correspond to any particular αΛ from the Λ-term case. We can see the same two nonstandard singularities in both cases, but no exponential solutions in the perfect fluid case. Formally one could retrieve them , but their emergence would be quite nontrivial and cannot be described within perfect fluid approach. This could indicate the fundamental difference between Λ-term and perfect fluid cases.
The regimes for Λ > 0, α < 0 are different for αΛ ≷ −5/6: for αΛ −5/6 there are two (one for exact equality) isotropic exponential solutions as future asymptotes while for αΛ < −5/6 there are anisotropic exponential solutions. The past asymptotes P (1,0) and K 3 are recovered perfectly while nS are replaced with P (0,1) -similar to α > 0 case. Of the future asymptotes, similar to the previous α > 0 D = 2 Λ-term case, anisotropic exponential solutions cannot be recovered -for the same reasons. On the contrary, isotropic exponential solutions could, with the similar technics as in D = 1 Λ-term case (see above).
To conclude, vacuum regimes could be reconstructed from ω < 1/3 perfect fluid regimes with some changes in the future asymptotes; ω > 1/3 regimes have no connections with vacuum regimes.
Λ-term regimes could be reconstructed only partially for D = 1 case; D = 2 case regimes are even less possible to reconstruct, in particular, Λ-term anisotropic exponential solutions could not be obtained from the perfect fluid approach.

VI. CONCLUSIONS
In this paper we have described the regimes which emerge in D = 1, 2 (the number of extra dimensions) flat cosmological models in EGB gravity with the matter source in form of perfect fluid. The results of the investigation suggest that in D = 1 there are no regimes with realistic compactification while in D = 2 there is. It is the transition from high-energy Kasner regime K 3 to exponential regime with expansion of three-and contraction of the remaining two-dimensional subspaces. It is very remarkable fact, as it could lead to a possible explanation of both the Dark Energy problem and features from effective string theory in the Early Universe within the same theory. Truth be told, we observed the late-time acceleration in both Λ-term and vacuum models; and if in the former of them accelerated expansion is not something amazing, in the latter it is, as in the GR accelerated expansion without any matter source is absent. Nevertheless, neither of these two could hardly describe the current state our Universe -we definitely observe existence of the ordinary matter. That is why the results of our current paper are important -we demonstrated that (at least in D = 2) the model with matter in form of perfect fluid and without Λ-term can have accelerated expansion phase at late times.
We shall proceed with the investigation of the same model in EGB but it higher dimensions, as well as to consider higher Lovelock contributions -this would give us answer on the question if this behavior is common and how much is it spread with respect to the parameters and initial conditions.
We have noted an interesting feature -in α < 0, D = 2 case there exists maximal density for the matter, and due to this fact for ω > 1/3 the past asymptote is GB-dominated, opposite to all other cases, where for ω > 1/3 the past asymptote is matter-dominated. In higher-dimensional EGB cosmologies the structure of the branches would be even more complicated so that the dynamics could be even more interesting and there could be other unexpected features, which makes this case more interesting to consider.