Realistic compactification in spatially flat vacuum cosmological models in cubic Lovelock gravity: Low-dimensional case

In this paper we begin to perform systematical investigation of all possible regimes in spatially flat vacuum cosmological models in cubic Lovelock gravity. We consider the spatial section to be a product of three- and extra-dimensional isotropic subspaces, with the former considered to be our Universe. As the equations of motion are different for $D=3, 4, 5$ and general $D \geqslant 6$ cases, we considered them all separately. Due to the quite large amount different subcases, in the current paper we consider only $D=3, 4$ cases. For each $D$ case we found values for $\alpha$ (Gauss-Bonnet coupling) and $\beta$ (cubic Lovelock coupling) which separate different dynamical cases, all isotropic and anisotropic exponential solutions, and study the dynamics in each region to find all possible regimes for all possible initial conditions and any values of $\alpha$ and $\beta$. The results suggest that in both $D$ cases the regimes with realistic compactification originate from so-called"generalized Taub"solution. The endpoint of the compactification regimes is either anisotropic exponential (for $\alpha>0$, $\mu \equiv \beta/\alpha^2<\mu_1$ (including entire $\beta<0$)) or standard low-energy Kasner regime (for $\alpha>0$, $\mu>\mu_1$), as it is compactification regime, both endpoints have expanding three and contracting extra dimensions. There are two unexpected observations among the results -- all realistic compactification regimes exist only for $\alpha>0$ and there is no smooth transition between high-energy and low-energy Kasner regimes, the latter with realistic compactification.


I. INTRODUCTION
It is interesting to note that the idea of extra dimensions is older then General Relativity (GR) itself. Indeed, the first ever extra-dimensional model was constructed by Nordström in 1914 [1], and it unified Nordström's second gravity theory [2] with Maxwell's electromagnetism. After Einstein proposed GR [3], it still took years before it was accepted: during the solar eclipse of 1919, the bending of light near the Sun was measured and the deflection angle was in perfect agreement with GR, while Nordström's theory, as most of the scalar gravity theories, predicted a zeroth deflection angle.
Though, the idea of extra dimensions was not forgotten -in 1919 Kaluza proposed [4] a similar model but based on GR: in his model five-dimensional Einstein equations could be decomposed into 4D Einstein equations and Maxwell's electromagnetism. For such decomposition to exist, the extra dimensions should be "curled" or compactified into a circle and "cylindrical conditions" should be imposed. The work by Kaluza was followed by Klein who proposed [5,6] a nice quantum mechanical interpretation of this extra dimension and so the theory, called Kaluza-Klein after its founders, was finalized. It is interesting to note that their theory unified all known interactions at that time. As a time flew, more interactions were known and it became clear that to unify all of them, more extra dimensions are needed. At present, one of the promising theories to unify all interactions is M/string theory.
One of the distinguishing features of M/string theories is the presence of the curvature-squared corrections in the Lagrangian. Scherk and Schwarz [7] demonstrated the presence of the R 2 and R µν R µν terms in the Lagrangian of the Virasoro-Shapiro model [8,9]; presence of the term of R µνλρ R µνλρ type was found in [10] for the low-energy limit of the E 8 × E 8 heterotic superstring theory [11] to match the kinetic term of the Yang-Mills field. Later it was demonstrated [12] that the only combination of quadratic terms that leads to a ghost-free nontrivial gravitation interaction is the Gauss-Bonnet (GB) term: This term, first discovered by Lanczos [13,14] (and so sometimes it is referred to as the Lanczos term), is an Euler topological invariant in (3+1)-dimensional space-time, but in (4+1) and higher dimensions it gives nontrivial contribution to the equations of motion. Zumino [15] extended Zwiebach's result on higher-than-squared curvature terms, supporting the idea that the low-energy limit of the unified theory might have a Lagrangian density as a sum of contributions of different powers of curvature. The sum of all possible Euler topological invariants, which give nontrivial contribution to the equations of motion in a particular number of space-time dimensions, form more general Lovelock gravity [16].
When one hears of the extra spatial dimensions, the natural question arises -where they are?
Our everyday experience clearly indicates there are three spatial dimensions, and experiments in physics and theory support this (for instance, in Newtonian gravity in more then three space dimensions there are no stable orbits, while we clearly see they are). The string theorists working with extra dimensions came with an answer -the extra spatial dimensions are compact -they are compactified on a very small scale, so small that we cannot sense them with our level of equipment.
But with that answer, another natural question comes to mind -how come that they are compact?
The answer to this question is not that simple. One of the ways to hide extra dimensions and to recover four-dimensional physics, is to consider so-called "spontaneous compactification" solution.
Exact static solutions with the metric as a cross product of a (3+1)-dimensional Minkowski spacetime and a constant curvature "inner space", were found for the first time in [17] (the generalization for a constant curvature Lorentzian manifold was done in [18]). For cosmology, it is more useful to consider the four-dimensional part given by a Friedmann-Robertson-Walker metric, and the size of extra dimensions to be time-dependent rather then static. In [19] it was demonstrated that in order to have a more realistic model one needs to consider the dynamical evolution of the extradimensional scale factor. In [18], the equations of motion with time-dependent scale factors were written for arbitrary Lovelock order in the special case of a spatially flat metric (the results were further proven and extended in [20]). The results of [18] were further analyzed for the special case of 10 space-time dimensions in [21]. In [22], the dynamical compactification was studied with use of the Hamiltonian formalism. More recently, searches for spontaneous compactifications were made in [23], where the dynamical compactification of the (5+1) Einstein-Gauss-Bonnet (EGB) model was considered; in [24,25] with different metric Ansätze for scale factors corresponding to (3+1)-and extra-dimensional parts. Also, apart from the cosmology, the recent analysis has focused on properties of black holes in Gauss-Bonnet [26][27][28][29]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.
If we want to find exact cosmological solutions, the most common Ansatz used for the scale factor is exponential or power law. Exact solutions with exponents for both the (3+1)-and extra-dimensional scale factors were studied for the first time in [41], and exponentially increasing (3+1)-dimensional scale factor and an exponentially shrinking extra-dimensional scale factor were described. Power-law solutions have been considered in [18,42] and more recently in [20,[43][44][45][46] so that by now there is an almost complete description of the solutions of this kind (see also [47] for comments regarding physical branches of the power-law solutions). Solutions with exponential scale factors [48] have also been studied in detail, namely, the models with both variable [49] and constant [50] volume; the general scheme for finding anisotropic exponential solutions in EGB was developed and generalized for general Lovelock gravity of any order and in any dimensions [51].
The stability of the exponential solutions was addressed in [52] (see also [53] for stability of general exponential solutions in EGB gravity), and it was demonstrated that only a handful of the solutions found and described in [51] could be called "stable", while the most of them are either unstable or have neutral/marginal stability.
In order to find all possible cosmological regimes in Einstein-Gauss-Bonnet gravity, one needs to go beyond an exponential or power-law Ansatz and keep the scale factor generic. We are especially interested in models that allow dynamical compactification, so that we consider the metric as the product of a spatially three-dimensional and extra-dimensional parts. In that case the three-dimensional part is "our Universe" and we expect for this part to expand while the extradimensional part should be suppressed in size with respect to the three-dimensional one. In [54] we demonstrated the there exist the phenomenologically sensible regime when the curvature of the extra dimensions is negative and the Einstein-Gauss-Bonnet theory does not admit a maximally symmetric solution. In this case both the three-dimensional Hubble parameter and the extradimensional scale factor asymptotically tend to the constant values. In [55] we performed a detailed analysis of the cosmological dynamics in this model with generic couplings. Recent analysis of this model [56] revealed that, with an additional constraint on couplings, Friedmann-type late-time behavior could be restored.
With the exponential and power-law solutions described in the mentioned above papers, another natural question arise -could these solutions describe realistic compactification or are they just solutions with no connection to the reality? To answer this question, we have considered the cosmological model in EGB gravity with the spatial part being the product of three-and extra dimensional parts with both subspaces being spatially flat. As both subspaces are spatially flat, the equations of motion could be rewritten in terms of Hubble parameters and then they become first order differential equations and could be analytically analyzed to find all possible regimes, asymptotes, exponential and power-law solutions. For vacuum EGB model 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.
The same analysis but for EGB model with Λ-term was performed in [59,60] and reanalyzed in [58]. The results suggest that the only realistic regime is the transition from high-energy GB Kasner to anisotropic exponential solution, it requires D 2, see [58][59][60] 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 these studies 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 [61]. The initial anisotropy could affect the results greatly -indeed, say, in vacuum (4 + 1)-dimensional EGB gravity with Bianchi-I-type metric (all the directions are independent) the only future asymptote is nonstandard singularity [44]. Our analysis [61] suggest that the transition from Gauss-Bonnet Kasner regime to anisotropic exponential expansion (with expanding three and contracting extra dimensions) is stable with respect to breaking the symmetry within both three-and 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 [61] for details).
The effect of the spatial curvature on the cosmological dynamics could be dramatic -say, positive curvature changes the inflationary asymptotics [62,63]. In the case of EGB gravity the influence of the spatial curvature 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 could be called a spiritual successor of [57][58][59][60] -now we are performing the same analysis but for cubic Lovelock gravity. In this paper we consider only D = 3, 4 (the number of the extra spatial dimensions) for vacuum case, other D cases, as well as Λ-term case and possible influence of anisotropy, spatial curvature and different kinds of matter source are to be considered in the papers to follow.
The manuscript is structured as follows: first we introduce Lovelock gravity and derive the equations of motion in the general form for the spatially-flat (Bianchi-I-type) metrics. Then we add our Ansatz and write down simplified equations. After that we describe the scheme we are going to use to analyze the particular cases. Then we consider particular cases with D = 3 and D = 4; for the former, we are going to describe the scheme step-by-step. In each section, dedicated to the particular case, we describe it and briefly summarize its features. Finally we summarize both cases, discuss their differences and similarities. After that we compare the dynamics in this cubic Lovelock with the dynamics in quadratic Lovelock (Einstein-Gauss-Bonnet) case, described in [57,58]. At last, we draw conclusions and formulate perspective directions for further investigations.

II. EQUATIONS OF MOTION
Lovelock gravity [16] has the following structure: its Lagrangian is constructed from terms where δ i 1 i 2 ...i 2n j 1 j 2 ...j 2n is the generalized Kronecker delta of the order 2n. One can verify that L n is Euler invariant in D < 2n spatial dimensions and so it would not give nontrivial contribution into the equations of motion. So that the Lagrangian density for any given D spatial dimensions is sum of all Lovelock invariants (1) upto n = D 2 which give nontrivial contributions into equations of motion: where g is the determinant of metric tensor, c n are coupling constants of the order of Planck length in 2n dimensions and summation over all n in consideration is assumed. The metric ansatz has the form g µν = diag{−1, a 2 1 (t), a 2 2 (t), . . . , a 2 n (t)}.
As we mentioned earlier, we are interested in the dynamics in cubic Lovelock gravity, so we consider n up to three (n = 0 is boundary term while n = 1 is Einstein-Hilbert, n = 2 is Gauss-Bonnet and n = 3 is cubic Lovelock contributions). Substituting metric (3) into the Lagrangian and following the usual procedure gives us the equations of motion: As mentioned in the Introduction, we want to investigate the particular case with the scale factors split into two parts -separately three dimensions (three-dimensional isotropic subspace), the dynamical equation that corresponds to h, and the constraint equation, Looking at (6) In previous papers, dedicated to study cosmological dynamics in EGB gravity, we performed analysis in all dimensions, sensitive to EGB case [57][58][59][60]. In the cubic Lovelock, the structure of the equations of motion is different in D = 3, 4, 5 and in the general D 6 cases. Also, since the current paper is dedicated to the vacuum case, we have Λ ≡ 0.

III. GENERAL SCHEME
The procedure of the analysis is exactly the same as described in our pervious papers [57][58][59][60] and is as follows: • we solve (8) with respect to H -one can see that (8)  If we take the discriminant of (8) with respect to H, and then its discriminant with respect to h, we obtain critical values for (α, β) which separate qualitatively different cases; • we find analitically isotropic exponential solutions: to do this we substituteḢ =ḣ ≡ 0 as well as h = H into (6)-(8); the system simplifies into a single equation and we solve it, finding not only roots but also the ranges of (α, β) where they exist; • we find analitically anisotropic exponential solutions: to do this we substitutė H =ḣ ≡ 0 into (6)-(8); the system could be brought to two equations: bi-six order poly-nomial in h with powers of α and β as coefficients and H = H(h, α, β). Both of them are usually higher-order with respect to their arguments so retrieving the solutions in radicand is impossible. But if we consider the discriminant of the former of them, the resulting equation gives us critical values for (α, β) which separate areas with different number of roots; • altogether first three items provides us with a set of critical values for (α, β) which separate domains with different dynamics; • we solve (6)-(7) with respect toḢ andḣ; • we substitute obtained H i curves intoḢ andḣ and obtain the latter as a single-variable functions:Ḣ(h) andḣ(h); • the obtainedḢ(h) andḣ(h) expressions are analyzed for all possible domains in (α, β) space to obtain all possible regimes; • obtained exponential regimes are compared with exact isotropic and anisotropic solutions (see [51]) to find the nature of the exponential regimes in question; • power-law regimes are analyzed in terms of Kasner exponents (p i = −H 2 i /Ḣ i ) to verify that low-energy power-law regimes are standard Kasner regimes with p i = p 2 i = 1 while high-energy power-law regimes are Lovelock Kasner regimes with p i = (2n − 1) = 5.
The described above scheme allows us to completely describe all existing regimes for a given set of the parameters (α, β). In the first case to consider, D = 3, we are going to completely describe all steps in detail.
Before continuing with the particular cases, it is necessary to 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 term (see, e.g., [20]). This way, for Einstein-Hilbert it is n = 1 and so p i = 1 (see [64]) and the corresponding regime is K 1 , which is usual low-energy regime in vacuum EGB case (see [58,60]) and we expect it to play the same role here. For Gauss-Bonnet it is n = 2 and so p i = 3 and the regime K 3 is typical high-energy regime for vacuum EGB case (again, see [58,60]). Finally, for cubic Lovelock it is n = 3 and so p i = 5 and the regime K 5 is expected to be typical high-energy regime for this case.
Another power-law regime is what we call "generalized Taub" (see [65] for the original solution).
It is the regime which was mistakenly taken for K 3 in [60], but in [58] it was corrected and explained (they both have p i = 3 which causes misinterpreting). It is a situation when for one of the subspaces the Kasner exponent p ≡ −H 2 /Ḣ 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.
We denote the exponential solutions 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. 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. The final regime is what we call "nonstandard singularity" and we denote is as nS. It is the situation which arise in Lovelock gravity due to its nonlinear nature. Since the equations (6)- (7) 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 then the denominator is equal to zero for finite values of H andor h. This situation is singular, as the curvature invariants diverges, but it happening for finite values of H andor h.
Our previous research demonstrate that this kind of singularity is wide spread in EGB cosmology -in particular, in totally anisotropic (Bianchi-I-type) (4 + 1)-dimensional vacuum cosmological model it is the only future asymptote [44].
As mentioned above, for this case we describe each step with details while in the cases to follow we skip the details. Finding the discriminant of the discriminant of (11) with respect to H gives us two critical values for µ ≡ β/α 2 ; µ 1 = 1/6, µ 2 = 3/2. These are two values for µ which qualitatively change the behavior of the H(h) curves.
To summarize, for α < 0 we have exponential solutions only for µ 3/2. roots ζ 1,2 < 0 and so h = ± αζ 1,2 ; for µ = 1/6 additional double root is added to the above roots (making total three different roots), for µ 3 > µ > 1/6 the double root from before is splitted making four different roots, at µ = µ 3 pairs of roots coinside leaving only two different roots and for µ > µ 3 roots degrade into nonstandard singularities leaving us with no roots at all. At µ 3/2 we have roots again, but they are negative, so that the roots for h = ± √ αζ are imagenary. Concluding, for α > 0 we have exponential solutions only for µ µ 3 .
Now we can solve (11) with respect to H and plot the resulting curves in Fig. 1. There red curve corresponds to H 1 , blue to H 2 and green to H 3 . The panels layout is as follows: One can see that on (a) and (d) panels (and so β < 0 and arbitrary α) we have one isotropic exponential  Before proceeding with the description of the regimes, one more important note should be taken. As could be seen from Fig. 1, some of the branches are discontinued-for instance, let us consider H 3 branch from Fig. 1(a) (α < 0, β < 0 case). One can see that it describe "internal" part of the loop but starting from some h > 0 it "jumps" into another branch of evolution. Obviously, this cannot happen in real physical evolution, so that it is the description which we use allows such "jumps", while the real physical evolution, for instance, for hyperbolic-like curve in the first quadrant of Fig. 1(a), is described partially by H 2 and partially by H 3 . Then to recover the real physical evolution, we "glue" the appropriate parts of the (Ḣ(h),ḣ(h)) curves. As a way of example, for the above mentioned hyperbolic-like curve from the first quadrant of Fig. 1(a), we plot (Ḣ(h),ḣ(h)) for H 2 branch in Fig. 6(a) (see Fig. 2(b)), for H 3 branch-in Fig. 6(b) (see Fig. 2(c)).
Then we notice the discontinuity in the Fig. 6(b)-it corresponds to the "jump" from the "inner" loop-like evolution curve to the curve under consideration at some h = h cr . So that for h < h cr we use H 2 part and for h > h cr -the H 3 part; the resulting glued curve is presented in Fig. 6(c).
One can verify that the resulting curve is free of any "jumps" and so represents physical evolution curve.
Exactly the same could be done for the analysis in terms of the Kasner exponents p i ≡ −H 2 i /Ḣ i , and the example for the same curve is presented in Figs. 6(d)-(f). Again, in Fig. 6(d) we presented Kasner exponents (red for p H , blue for p h and green for p i ) for H 2 branch (see Fig. 4(b)), in The resulting regimes are: for α < 0, β < 0 case (Figs. 1(a), 2(a)-(c) and 4(a)-(c)) we have P 1,0 ↔ E iso and P 0,1 ↔ E iso on hyperbolic-like curve and K 1 ↔ nS as well as nS ↔ nS on the inner loop-like curve. The regimes P 1,0 and P 0,1 are the regime with power-law behavior for one of the Hubble parameters with p i = 1 and "effective" p j = 0 for another, hence the designation: for P 1,0 we have p H = 1, p h = 0 while for P 0,1 we have p H = 0, p h = 1. Let us note that for D = 3 there is no real difference between P 1,0 and P 0,1 (since both subspaces are three-dimensional) but for future cases there is, so we treat these two cases separately. Overall, in the α < 0, β < 0 case there are no regimes with realistic compactifications.
One cannot miss the moment when H 1 and H 2 branches "touch" each other at µ = 3/2; at that point two isotropic exponential solutions E 1,2 iso coincide. For µ > 3/2 the branches "detach" each other, forming banana-like curves as presented in Fig. 1(c). The isotropic exponential solutions degrade into anisotropic ones, located on these banana-like curves. Similar to the previous case, the asymptotes for H and h are ± −α/10β and the corresponding regimes are nonstandard singularities. Then, combining all above mentioned with the analysis of the (Ḣ(h),ḣ(h)) curves allows us to conclude the regimes: P 0,1 ↔ E 1 , P 1,0 ↔ E 2 and nS ↔ E 1,2 on banana-like curves and nS ↔ K 1 on central curve. So that in this case we have P 1,0 , P 0,1 ↔ E 1,2 as a nonsingular regime, but E 1,2 are located in the first and third quadrants, meaning either (H > 0, h > 0) or (H < 0, h < 0), so that both three-and extra-dimensional spaces are expanding or contracting, which contradict our choice (H > 0, h < 0), so that we discard this regime and cannot call it realistic.
We go on with α > 0, β < 0 ( Fig. 1(d), 3(a)-(c) and 5(a)-(c)) case. Similar to two previous cases, there are three limiting values for H as h → ±∞ (as well as for h as H → ±∞) and they are the same as for the previous case (± −α/10β), and, also similar to the previous cases, the corresponding regimes are nonstandard singularities. The anisotropic exponential solutions are located on edge-shaped curves in second and fourth quadrants while isotropic-on hyperbola-like curve in first and third quadrants. Combining all these with the analysis of the (Ḣ(h),ḣ(h)) curves allows us to draw the regimes: E 2 ↔ P 1,0 , E 1 ↔ P 0,1 and E 1,2 ↔ nS on the edge-shaped curves, K 1 ↔ nS and nS ↔ nS on central loop-like curve and nS ↔ E iso on hyperbola-like curve.
In this case we have two interesting features which worth mentioning-first, we finally have realistic compactification-indeed, P 1,0 , P 0,1 → E 1,2 are realistic compactifications, as both E 1,2 have different signs for H and h 1 . The second feature is inability to reach E iso from the standard singularity-indeed, isotropic exponential solution is connected only to nonstandard singularity. This feature does not affect realistic compactification abundance, but we note it for completeness; also, in the Gauss-Bonnet case [57][58][59][60] we never experienced such situation, neither in vacuum nor in Λ-term cases.
The difference, as one can see from Figs. 1(e, f) lies in circle-like curve in the second and fourth quadrants. The regime which is common for both cases is P 1,0 ↔ K 1 , and it is another example of the realistic compactification. The regimes within the circle-like curve in the second and fourth quadrants are subject to the "fine-structure" (see the description of the anisotropic exponential solutions abundance above) and are presented in Fig. 7. There we presented the following α > 0,  coinside with our description provided above.
For visualisation purposes we plot all discovered regimes on H(h) curves, the results are presented in Fig. 8. The panels layout follow that of Fig. 1: β > 0 and µ < 0.3 on (e) panel; red curve corresponds to H 1 , blue -to H 2 and green -to H 3 . The arrows represent the evolution according with respect to grow of the cosmic time t.
This concludes our consideration of the first and the simplest D = 3 case. We concluded that there are two sets of realistic regimes: P 1,0 , P 0,1 → E 1,2 for α > 0, β < 0 and P 1,0 , P 0,1 → K 1 for α > 0, β > 0; interesting that both of them require α > 0. Let us note that only half of these regimes have H > 0, h < 0 -the other half have H < 0, h > 0, but since this particular case has D = 3 -the number of extra dimensions coincide with the number of "our" spatial dimensions, we can "switch" between them, which effectively doubles the number of regimes. One also can notice 6Ḣ + 12H 2 + 6ḣ + 12h 2 + 18Hh + 8α 3Ḣ(H 2 + 6Hh + 3h 2 ) + 3ḣ(3H 2 + 6Hh+ (17) First we find the discriminant of (17) with respect to H and then its discriminant with respect to h. The resulting equation is 15th order polynomial with respect to µ = β/α 2 and it has the following roots: single roots To find isotropic exponential solutions, we substituteḢ =ḣ ≡ 0 as well as h = H into (15)-(17), the system is simplified to a single equation Analyzing (18) we find out that there is one root if β < 0 (regardless of α) and two roots if α < 0, 5α 2 /6 β > 0; in all other cases there are no isotropic exponential solutions.
One can see that in D = 4 case, unlike D = 3, we have "proper" cubic-Lovelock Kasner regime with p = 5, as predicted. But as E 1 is located along H 1 , it has either (H < 0, h > 0) or (H > 0, h < 0). Only one of them is of interest to us (H > 0, h < 0), but this regime is reachable only as past asymptote: E 1 → K 5 ; the E 1 from K 5 → E 1 has (H < 0, h > 0). So that the regime which could give us realistic compactification, cannot be reached from K 5 , despite it exists (but unstable).
Singular regimes for α < 0, µ < 0 case include E iso ↔ nS, E 2 ↔ nS, nS ↔ nS and K 1 ↔ nS along H 2 − H 3 curve and nS ↔ E 1 , nS ↔ nS and K 1 ↔ nS along H 1 curve. So that another anisotropic exponential solution E 2 is located between nonstandard singularities and cannot be reached from the initial cosmological singularity.
Next case is α < 0, β > 0, µ < 5/6, which is presented in Fig. 9(b). There the only nonsingular regimes are K 5 ↔ K 1 (along H 3 , green line) and P 1,0 ↔ E 1 iso . Similar to D = 3, H 1 and H 2 branches have fixed limit at h → ±∞, the corresponding regime is nonstandard singularity. This makes E 2 iso to be connected only with nS; apart from these regimes there is also K 1 ↔ nS. So that in this case we have K 5 → K 1 which is the transition from high-energy to low-energy Kasner, similar to the situation described in the Einstein-Gauss-Bonnet case [57]. But similar to the just described above the situation with the anisotropic exponential solution, K 1 in this case has either (H < 0, h > 0) or (H > 0, h < 0). And again, similar to the previously described case, for K 5 → K 1 we have (H < 0, h > 0) while for K 1 → K 5 we have (H > 0, h < 0). So that the Kasner regime which is suitable for us (three expanding and D contracting dimensions) is unstable as t → ∞. We shall discuss it more in the proper section of the manuscript. Let us also note that µ = 5/6 is the limiting case of (b) panel where the isotropic exponential solutions coincide.
The nonstandard singularities in this case have power-law behavior, that is why they are denoted as nS/P -they have singularity at H = ± −10α/β/30 and they approach it in a power-law manner in a finite time. So that the singularity remains nonstandard, yet, unlike other appearances, in this particular case its behavior is know, so that we specify it and denote as nS/P . Now let us consider α < 0, β > 0, µ > 5/6 case, which is depicted in Fig. 9(c). Similar to the previous case, as h → ±∞ the values of H tend to the same constants and the corresponding regimes are again nS/P for the same reasons as in the previous µ < 5/6 case. Again, similar to the previous case, we have K 5 along H 3 (green line), but now isotropic exponential solutions "turned" into anisotropic, located on "banana"-shaped curves. So, similar to the previous case, the we have K 5 ↔ K 1 nonsingular regime and it is, again, have H < 0 for K 1 so it is invalid for our purpose of finding proper compactification. But unlike the previous case, we also have P 1,0 ↔ E 2 -transition to the anisotropic exponential solution-another nonsingular regime. But this anisotropic solution is located in the first and third quadrants and so it has either (H > 0, h > 0) or (H < 0, h < 0)so both three and extra dimensions are either expanding or contracting and so we cannot explain compactification with this regime. The singular regimes for this case include E 1 ↔ nS/P , E 2 ↔ nS and K 1 ↔ nS.
So that in these three cases which combine entire α < 0 domain there is no realistic compactification. Let us continue with α > 0 regimes. First one to consider is α > 0, β < 0, presented in Fig. 9(d). There we have isotropic exponential solution on H 2 (blue hyperbola-like curve) and K 5 on H 1 (red line), but unlike previous cases there is no K 5 ↔ K 1 transition, as there is anisotropic exponential solution on H 1 as well, so we have K 5 ↔ E 1 instead. And again, similar to one of the previous cases, this exponential solution has H < 0 on K 5 → E 1 and H > 0 on E 1 → K 5 , so that we cannot describe realistic compactification with this regime. Similar to D = 3 case the isotropic exponential solution is "surrounded" by nonstandard singularities -the regimes involve E iso are E iso ↔ nS and E iso ↔ nS/P . Another non-singular regime is P 1,0 ↔ E 2 which takes place on H 3 "edge"-shaped part, and in this case E 2 has H > 0, h < 0 and is stable past-time asymptote, so that it could describe the compactification. The singular regimes in this case include E 1 ↔ nS, E 2 ↔ nS, K 1 ↔ nS/P and K 1 ↔ nS. To conclude, this case provides us with P 1,0 → E 2 which could describe realistic compactification.
The last case to consider is α > 0, β > 0, and it has wide variety of the regimes, as well as fine structure, similar to the previous D = 3 case. But unlike it, now all the regimes from the fine structure are connected with the initial singularity, making them more physical. We presented two representative cases -µ < µ 1 and µ > µ 4 = 1/6 in Figs 9(e) and 9(f) respectively, while all cases inbetween -in Fig. 10.
Let us start the description with growth of µ.
The first case is µ < µ 1 , presented in Fig. 9(e). The regimes along H 3 (green curve) in the second quadrant are: and green -to H 3 , in accordance with the designation in Fig. 9) (see the text for more details). like figure formed from H 1 and H 2 branches in the center are (starting from red H 1 branch): where first K 1 belong to H 1 branch while the last -to H 2 . One can see that there are two anisotropic exponential solutions along this curve, but both of them are unstable. So that in this case only P 1,0 → E 2 is compactification regime. The next case is µ = µ 1 -in that case H 2 (blue curve) and H 3 (green curve) "touch" each other, but the regimes remain unchanged. With further increase of µ the H 2 and H 3 branches "detouch" each other but in a way presented in Fig. 10(a), where we presented finestructure (and so only the second quadrant) for µ 3 > µ > µ 1 . One can see that now H 3 is "cut" into two physical branches -one of them together with a piece from H 2 is forming P 1,0 → K 1 regime (see also Fig. 9(f) where it seen clearly). Another part of H 3 is forming with remaining parts of H 2 and H 1 the second physical branch with the regimes (starting from K 1 on H 1 ): One can see there four different anisotropic exponential solution (similar as in D = 3 case, see Fig. 7(c)) but now located along the same physical branch. Three out of them are unstable (E 1 , E 2 and E 4 ) and the remaining stable E 3 is "surrounded" by nS and cannot be a part of compactification scenario. To conclude, the only compactification regime here is P 1,0 → K 1 .
The further growth of µ to µ = µ 3 is presented in Fig. 10(b) -the situation and the regimes are almost the same with the difference that now an additional exponential solution (E 4 by designations on Fig. 10(b)) emerged "between" two nonstandard singularities on the upper part of H 3 . Since it emerged between nS, no new physically significant regimes are formed. With increase of µ this exponential solution, which has H i = 0 (constant-volume exponential solution, see [50]) at µ = µ 3 , is split into two for µ 8 > µ > µ 3 -one of them (E 2 ) is stable while another (E 3 ) is unstable (see Fig. 10(c) for illustration). But again, both of them are surrounded by nonstandard singularities so no new physically sensitive regimes appear.
The next qualitative change of the situation happening at µ = µ 8 , displayed in Fig. 10 no changes of the realistic regimes occur. With further growth of µ to µ 4 > µ > µ 8 these newly formed exponential solutions turn to nonstandard singularities (see Fig. 10(e)); the remaining two exponential solutions (E 1 and E 2 on the H 3 ) merge with nonstandard singularity between them into new exponential solution at µ = µ 4 (see Fig. 10(f)). Finally, at µ > µ 4 ( Fig. 9(f)) no exponential solution remains and there are only nS ↔ nS and K 5 ↔ nS along H 2 − H 3 physical branch. So that for the entire µ > µ 1 we have only P 1,0 → K 1 as realistic compactification regime; all the regimes within fine structure are non-realistic.
To conclude, in D = 4 we encounter proper cubic Lovelock Kasner regime K 5 with p = 5, but there are no realistic compactification regimes originating from K 5 . Instead, there are P 1,0 → E 3+4 for α > 0, µ µ 1 and P 1,0 → K 1 for α > 0, µ > µ 1 . So that for entire α > 0 we have either of these regimes, and, similar to the previous cases, have no realistic regimes for α < 0.

VI. DISCUSSIONS
In the current paper we have analyzed the cosmological dynamics of the cubic Lovelock gravity, with Einstein-Hilbert and Gauss-Bonnet terms present as well. We have chosen a setup with a topology being a product of two isotropic subspaces -three-dimensions, representing our Universe, and D-dimensional, representing extra dimensions. Both subspaces are flat, which simplifies our equations of motion and makes it possible to analyze them analytically. In a sense, it is a logical continuation of [57][58][59][60], where we considered the same problem but in EGB gravity -vacuum case in [57,58] and Λ-term case in [58][59][60]. In [58] we reviewed all the results for EGB from [57,59,60] and changed the visualization of the regimes -in the original papers [57,59,60] we use tables to list of all the regimes, and this way sometimes is not easy to read. On the contrary, in [58] we put all the regimes on H(h) curves and added arrows to demonstrate t → ∞ directed evolution. In the current paper we decided to keep visualization from [58].
First of all, let us summarize the results, as they are scattered over mini-conclusions in each particular section. The fist one is D = 3 and it has interesting feature -since the equations of motion are cubic in both H and h, there could be up to three branches of the solutions. On the other hand, it is the lowest possible dimension for cubic Lovelock gravity, so there are no Kasner solutions (see [47]). Then the only possibility is what we call "generalized Taub" solution -the situation when the expansion in each direction is characterized by Kasner exponent p i = −H 2 i /Ḣ i equal to either 1 or 0; so that for our topology it is either P 1,0 (p H = 1, p h = 0 -expansion of the three-dimensional subspace and "static" extra dimensions) or P 0,1 (p H = 0, p h = 1 -expansion of the extra-dimensional subspace and "static" three dimensions). Then the remaining brancheswhich cannot be connected to either P 1,0 or P 0,1 , form closed evolution curves for (α < 0, β < 0) (see Fig. 8(a)) or (α > 0, β > 0) (see Figs. 8(e) and 7); for (α < 0, β > 0) and (α > 0, β < 0) they encounter nonstandard singularities (see Figs. 8(b)-(d)). The realistic compactification regimes are P 1,0 /P 0,1 → E 3+3 for (α > 0, β < 0) and P 1,0 /P 0,1 → K 1 for (α > 0, β > 0); let us note that both of the regimes exist only for α > 0.
The above-mentioned "generalized Taub" solution deserves some additional comments. Formally it fits "generalized Milne" solution -the second branch of the power-law solutions in Lovelock gravity (see [20] for details), but it is only formal -it fits only because it is degenerative. As it was demonstrated in [47], strict "generalized Milne" cannot exist in pure highest-order Lovelock gravity, as it leads to degeneracy in the equations of motion. But if additional (lower-order) Lovelock contributions are involved, it could restore this branch of solutions, but it was never demonstrated exactly before. So that on the particular example of [3 + D] spatial splitting we demonstrated this possibility. Still, a little is known about this regime and it deserves additional investigation in the separate papers.
When we consider this "generalized Taub" solution as a past asymptote -and this is the case for all possible realistic compactification models in D = 3, 4 -it feels a bit unnatural. Indeedthe P 1,0 regime imply H → ∞ and h → 0 as t → 0 (by "0" we mean here initial cosmological singularity), so that we initially have "burst"-like expansion of three-dimensional subspace while the extra-dimensional subspace is almost static. In addition to the feeling of unnaturalness, it is a question if this regime could be reached from totally anisotropic space, in a manner it was done in [61] for EGB case. So that it gives additional reason to deeply investigate this regime and we are going to do it in the nearest time.
The results of our analysis suggest that the variety and abundance of the regimes is closer to Λ-term EGB, rather then vacuum EGB models. The reasons for that are not clear, but we expect that number of the free parameters plays a role here. Indeed, for vacuum EGB model there is only one parameter -α, Gauss-Bonnet coupling, while for Λ-term EGB and vacuum cubic Lovelock there are two -α and Λ for the former and α and β (cubic Lovelock coupling) for the latter. In that case the dynamics of the Λ-term cubic Lovelock gravity would be even more interesting and we are going to consider this case shortly.

VII. CONCLUSIONS
This concludes our study of the cosmological models in vacuum cubic Lovelock gravity with D = 3, 4 extra dimensions. We have found that in both cases there are regimes with successful compactification, but all of them originate from "generalized Taub" solution; for the future asymptote we have either Kasner regime or anisotropic exponential solution.
Apart from the regimes with successful compactification, we described and plotted on H(h) curves all possible regimes for all initial conditions and all structurally different cases. The variety and abundance of the regimes exceed even Λ-term EGB case, featuring transition between two anisotropic exponential solutions and transition between two different "generalized Taub" solutions.
There are two interesting observations which require additional investigation, as both are quite unexpected. First of them is that all of the regimes with realistic compactification have α > 0 requirement. This is totally unexpected, as in both vacuum and Λ-term EGB cases we have viable compactifications for both signs of α. For the Λ-term case the joint analysis of our cosmological bounds and those coming from AdS/CFT and other considerations allows us to conclude α > 0 (see [58,59]), but for that we involved external (to our results) analysis. In this case without any external bounds we already have realistic compactification only for α > 0.
The second observation is that there is no K 5 → K 1 transition with realistic compactification.
In EGB vacuum case [57,58] we have the transitions of this kind, so it was natural to assume that in higher-order Lovelock gravity they also present, but our investigation reveals that they are not.
There is K 5 → K 1 transition, but with contracting three and expanding extra dimensions, so it formally exist, but with no compactification scenario. As both of these observations are unexpected and in disagreement with what we have learned from study of EGB case, this is a good direction for further improvement of our understanding of Lovelock gravity.