Cosmological solutions in Einstein–Gauss–Bonnet gravity with static curved extra dimensions

Abstract

In this paper we perform systematic investigation of all possible solutions with static compact extra dimensions and expanding three-dimensional subspace (“our Universe”). Unlike previous papers, we consider extra-dimensional subspace to be constant-curvature manifold with both signs of spatial curvature. We provide a scheme how to build solutions in all possible number of extra dimensions and perform stability analysis for the solutions found. Our study suggests that the solutions with negative spatial curvature of extra dimensions are always stable while those with positive curvature are stable for a narrow range of the parameters and the width of this range shrinks with growth of the number of extra dimensions. This explains why in the previous papers we detected compactification in the case of negative curvature but the case of positive curvature remained undiscovered. Another interesting feature which distinguish cases with positive and negative curvatures is that the latter do not coexist with maximally-symmetric solutions (leading to “geometric frustration” of a sort) while the former could – this difference is noted and discussed.

1 Introduction

Einstein–Gauss–Bonnet Gravity (EGB) is the simplest example of a larger family of gravity theories known as Lovelock Gravities [1]. Lovelock gravities are characterized by the fact that their actions possess higher power curvature terms but whose variation lead to equations of motion which remain of second order derivative in the metric. Lovelock gravities are therefore the most natural generalization of General Relativity to higher space-time dimensions. EGB gravity whose action has additional term quadratic in the curvature with respect to general relativity exists for space-time dimensions $d\ge 5$ (in four dimensions the quadratic Gauss Bonnet term does not affect the equations of motion being topological). It is the most studied Lovelock gravity in literature as it has many of the features of more generic Lovelock gravities but keeping a relatively simple action. Moreover EGB gravity can also be seen as the low energy limit of certain string theories [2].

EGB gravity, whose action holds $\Lambda$-term, Einstein–Hilbert term and quadratic Gauss-Bonnet term, has the remarkable feature that it can possess up to two independent maximally symmetric solutions (even with different sign of the curvature scale). Indeed, in opposition to General Relativity, the curvature scale of the maximally symmetric solutions is not determined only by the $\Lambda$-term but is a function of all three couplings of the theory. This can be immediately seen by plugging the ansatz of a maximally symmetric space-time in the equations of motion of EGB gravity. One gets a quadratic equation in the scale and therefore up to two independent maximally symmetric solutions. Indeed also for the EGB black hole solutions one gets in general two branches with different asymptotic behavior [3]. However the discriminant of the quadratic equation for the curvature scale can also be negative for a range of values of the coupling constants. In this case there exist no maximally symmetric space-time solution at all. This situation remained almost unexplored in literature due to the fact that it is difficult to give reasonable asymptotic falloff conditions in such cases. However this situation of non existing maximally symmetric solutions can be of interest in the context of dynamical compactification in cosmology. Indeed one may wonder why a space-time should tend to a compactified manifold with three large dimensions and D compact dimensions instead of tending to an isotropic space-time which seems more natural. The non-existence of a maximally-symmetric solution would force the space-time to search for a less symmetric configuration and give a very natural explanation for compactification. This has been explored, to the best knowledge of the authors for the first time in [4,5,6]. The situation where the EGB theory does not admit maximally symmetric solutions was called in these papers “geometric frustration”. The term “geometric frustration” is normally used in the context of condensed matter physics when a system can not take the configuration of minimal energy due to topological obstructions. In the previously cited papers, for simplicity it was assumed that space-time had the structure of a warped product of a 4-dimensional flat FRW space-time with a D dimensional constant curvature space with an independent scale factor. In the case that curvature of the compact dimensions is negative and moreover the couplings of the EGB are chosen from the open region of couplings space where geometric frustration occurs it was shown that realistic compactification scenarios exist where the scale factor of the extra dimensions tends to a constant (i.e., stabilizes). No realistic compactification scenario was found for the case when the extra dimensions have positive curvature or when there is no geometric frustration. Considering that the cited papers were focused on a generic number D of extra dimensions and especially to the large D limit it is still possible that some cases with positive curvature of the extra dimensions for some particular value of D have just been overseen.

It is interesting to note that the presence of higher-order curvature terms in the Lagrangian is one of the features of string-inspired theories. Historically, Scherk and Schwarz [7] were the first to demonstrate the presence of the $R^2$ and $R_{\mu \nu } R^{\mu \nu }$ terms in the Lagrangian of the Virasoro–Shapiro model [8, 9]. A presence of curvature-squared term of the $R^{\mu \nu \lambda \rho } R_{\mu \nu \lambda \rho }$ types was demonstrated [10] for the low-energy limit of the $E_8 \times E_8$ heterotic superstring theory [11] to match the kinetic term for the Yang–Mills field. Later it was demonstrated [2] that the only combination of quadratic terms that leads to a ghost-free nontrivial gravitation interaction is the Gauss–Bonnet (GB) term:

$$\begin{aligned} L_{GB} = L_2 = R_{\mu \nu \lambda \rho } R^{\mu \nu \lambda \rho } - 4 R_{\mu \nu } R^{\mu \nu } + R^2. \end{aligned}$$

This term, initially discovered by Lanczos [12, 13] (therefore it is sometimes referred to as the Lanczos term) is an Euler topological invariant in (3+1)-dimensional space-time, but not in (4+1) and higher dimensions. Zumino [14] extended Zwiebach’s result on higher-than-squared curvature terms, supporting the idea that the low-energy limit of the unified theory should have a Lagrangian density as a sum of contributions of different powers of curvature. In this regard the Einstein–Gauss–Bonnet (EGB) gravity could be seen as a subcase of more general Lovelock gravity [1], but in the current paper we restrain ourselves with only quadratic corrections and so to the EGB case.

Generally speaking, all extra-dimensional theories have one thing in common – we need to explain where additional dimensions are “hiding” – since we do not sense them, at least with the current level of experiments. One of the possible ways to hide extra dimensions and to recover four-dimensional physics, is to build a so-called “spontaneous compactification” solution. Exact static solutions with the metric set as a cross product of a (3+1)-dimensional manifold and a constant curvature “inner space”, were found for the first time in [15], where (3+1)-dimensional manifold being Minkowski (the generalization for a constant curvature Lorentzian manifold was done in [16]). In the context of cosmology, it is more useful to consider spontaneous compactification with the four-dimensional part given by a Friedmann–Robertson–Walker (FRW) metric. In this case it is also natural to consider the size of the extra dimensions being time dependent rather than static. Indeed, in [17] it was exactly demonstrated that in order to have more realistic model one needs to consider the dynamical evolution of the extra-dimensional scale factor. In [16], the equations of motion with time-dependent scale factors were written down for arbitrary Lovelock order in the special case of spatially flat metric (the results were further proven in [18]). The results of [16] were further analyzed for the special case of 10 space-time dimensions in [19]. In [20], the dynamical compactification was studied with the use of Hamiltonian formalism. More recently, searches for spontaneous compactifications were performed in [21], where the dynamical compactification of the (5+1) Einstein–Gauss–Bonnet model was considered; in [22, 23] with different metric Ansätzen for scale factors corresponding to (3+1)- and extra-dimensional parts; and in [4,5,6] (mentioned above), where general (e.g., without any Ansätz) scale factors and curved manifolds were considered. Also, apart from cosmology, the recent analysis was focused on properties of black holes in Gauss-Bonnet [3, 24,25,26,27] and Lovelock [28,29,30,31,32] gravities, features of gravitational collapse in these theories [33,34,35], general features of spherical-symmetric solutions [36], and many others.

When we are looking for exact cosmological solutions, two main ansätzen are involved – power-law and exponential. One of first approaches to power-law solutions in EGB gravity performed in [16, 37] and more recently they were studied in [18, 38,39,40,41]. One of the first approaches to the exponential solutions was done in [42], the recent works include [43,44,45]. We separately described the exponential solutions with variable [46] and constant [47] volume; let us note [48] 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. General scheme for finding exponential solutions in arbitrary dimensions and with arbitrary Lovelock contributions taken into account described in [49]. Deeper investigation revealed that not all of the solutions found in [49] could be called “stable” [50]; see also [44] for more general approach to the stability of exponential solutions in EGB gravity.

The simplest case – when the spatial section is the product of spatially-flat three- and extra-dimensional subspaces, systematic study of all possible regimes in EGB and partially in cubic Einstein-Lovelock gravity was performed in [51,52,53,54,55,56]. In particular, for vacuum EGB case it was done in [51] and reanalyzed in [52]. We also added cubic Lovelock term and analyzed the resulting vacuum Einstein–Lovelock cosmology in [53, 54]. Similar analysis for EGB model with $\Lambda$-term was performed in [55, 56] and reanalyzed in [52]. All these studies demonstrate that there are exponential regimes with expanding three and contracting extra dimensions and they are not suppressed. On the contrary, it is relatively difficult to reach power-law regime naturally.

In the studies described above we made two important assumptions – both subspaces (three- and extra-dimensional) were considered to be spatially flat and isotropic. But neither of these conditions could be called “natural”, so it is interesting to investigate what happens if we left these conditions? In the Friedmann cosmology spatial curvature plays important role, for example, positive curvature changes the possibility to reach inflationary asymptotic [57, 58]. As we previously mentioned, in EGB gravity the influence of the spatial curvature was studied in [4, 5], where we described “geometric frustration” regime and further investigated it in [6].

We addressed effects of both curvature and anisotropy earlier in [59]. Particularly, we considered initially totally anisotropic (Bianchi-I-type) $(5+1)$- and $(6+1)$-dimensional models and numerically studied their evolution. The former of them has only one stable anisotropic exponential solution – with expanding three and contracting two dimensions – and it is the only dynamical attractor of the system. On the contrary, the latter has two possibilities – expanding three and contracting three or expanding four and contracting two dimensions, and depending on the initial conditions we could end up in both of the possibilities. So that if (in)appropriate exponential solution exists and stable, initially anisotropic Universe could end up with “wrong” compactification.

We can also note that apart from vacuum and $\Lambda$-term models we considered models with perfect fluid as a source: initially we considered them in [43], some deeper studies of (4+1)-dimensional Bianchi-I case was done in [40] and deeper investigation of power-law regimes in pure GB gravity in [41]. Systematic study for all D was started in [60] for low D and is currently continued for high D cases.

The aim of this paper is actually to check again the compactification, and to do it also for smaller values of D and see if for some particular value there exist compactification with stabilized extra dimensions. The interest in checking again the positive curvature in extra dimensions is that from a point of view of Kaluza-Klein theory the positive curvature of extra dimensions allows to include non-abelian interactions. In order to do so we check separately the existence of solutions and then their stability due to the fact that if solution exist but is unstable it is useless from a practical point of view. The stability issue of the solutions was not addressed in the previous papers.

The structure of the manuscript is as follows: first we write down equations of motion and then consider several particular cases which differs from each other (for several low D – number of extra dimensions – equations of motion simplify via dropping some of the terms), ending with general case (which has all possible terms). For each of these cases we obtain a solution, write down perturbed equations and solve them around the found solution. After all cases described, we summarize the results, discuss them and draw conclusions.

2 Equations of motion

We start with the standard Einstein–Gauss–Bonnet Lagrangian in the cosmological background (see, e.g., [18])

$$\begin{aligned} {{\mathcal {L}}} = R + \alpha {{\mathcal {L}}}_{GB} - \Lambda , \end{aligned}$$

(1)

where R is the Ricci scalar, $\Lambda$ is $\Lambda$-term (or boundary term) and ${{\mathcal {L}}}_{GB}$,

$$\begin{aligned} {{\mathcal {L}}}_{GB} = R_{\mu \nu \alpha \beta } R^{\mu \nu \alpha \beta } - 4 R_{\mu \nu } R^{\mu \nu } + R^2 \end{aligned}$$

(2)

is the Gauss–Bonnet Lagrangian while $\alpha$ is its coupling constant. We consider space-time to be warped product of Lorenzian and constant curvature manifolds. The former is $(3+1)$-dimensional while the latter is D-dimensional. Then its metric could be written as

$$\begin{aligned} g_{\mu \nu }= & {} \mathop {\mathrm{diag}}\nolimits \Bigg \{ -1, a^2(t), a^2(t), a(t)^2, b(t)^2, b(t)^2 \chi ^2(x_4) \nonumber \\&\ldots , b^2(t) \displaystyle {\prod \limits _{i=4}^{D-4}\chi ^2(x_i)} \Bigg \}, \end{aligned}$$

(3)

where $\chi (x) = \sin (x)$ for positive and $\chi (x) = \sinh (x)$ for negative curvature of extra dimensions.

Now we modify the metric above by inserting lapse function N(t) ($-1 \mapsto -N^2(t)$ and substitute the modified metric into the Lagrangian (2), we calculate the action and vary it with respect N(t), a(t) and b(t) to obtain equations of motion: constraint and dynamical equations for a(t) and b(t), respectively:

$$\begin{aligned}&D(D-1) \left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) + 6DHH_b + 6H^2\nonumber \\&\quad + \alpha \left( D(D-1)(D-2)(D-3)\right. \nonumber \\&\quad \left. \times \left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) ^2 + 24DH^3H_b \right. \nonumber \\&\quad \left. + 12D(D-1)H^2(H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}) + 24D(D-1)H^2H_b^2\right. \nonumber \\&\quad \left. + 12D(D-1)(D-2)H^2\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right) = \Lambda , \nonumber \\&D(D-1)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) + 4(\dot{H} + H^2)\nonumber \\&\quad + 2D(\dot{H}_b + H_b^2) + 4DHH_b + 2H^2 \nonumber \\&\quad + \alpha \left( D(D-1)(D-2)(D-3)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) ^2\right. \nonumber \\&\quad \left. + 16D(\dot{H} + H^2)HH_b\right. \nonumber \\&\quad \left. + 8D(D-1)(\dot{H} + H^2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\quad \left. + 8D(\dot{H}_b + H_b^2)H^2 + 16D(D-1)(\dot{H}_b + H_b^2)HH_b\right. \nonumber \\&\quad \;\, \left. + 4D(D-1)(D-2)(\dot{H}_b + H_b^2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\quad \left. + 4D(D-1)H^2\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\left. \quad +8D(D-1)(\dot{H} + H^2)H_b^2\right. \nonumber \\&\quad \left. + 8D(D-1)(D-2)HH_b\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right) = \Lambda , \nonumber \\&(D-1)(D-2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) + 6(\dot{H} + H^2)\nonumber \\&\quad + 2(D-1)(\dot{H}_b + H_b^2) + 6(D-1)HH_b + 6H^2 \nonumber \\&\quad + \alpha \left( (D-1)(D-2)(D-3)(D-4)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) ^2\right. \nonumber \\&\quad \left. + 24(\dot{H} + H^2)H^2 \right. \nonumber \\&\quad \left. + 12(D-1)(D-2)(\dot{H} + H^2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\left. \quad + 48(D-1)(\dot{H} + H^2)HH_b + 24(D-1)(\dot{H}_b + H_b^2)H^2 \right. \nonumber \\&\quad \left. + 4(D-1)(D-2)(D-3)(\dot{H}_b + H_b^2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\quad \left. + 24(D-1)(D-2)(\dot{H}_b + H_b^2)HH_b \right. \nonumber \\&\quad \left. + 12(D-1)(D-2)(\dot{H} + H^2)\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\quad \left. + 24(D-1)H^3H_b \right. \nonumber \\&\quad \left. + 12(D-1)(D-2)(D-3)HH_b\left( H_b^2 + \displaystyle \frac{\gamma _D}{b(t)^2}\right) \right. \nonumber \\&\quad \left. + 24(D-1)(D-2)H^2H_b^2 \right) = \Lambda , \end{aligned}$$

(4)

where $H \equiv \dot{a}(t)/a(t)$ is the Hubble parameter associated with “ordinary” space, $H_b \equiv \dot{b}(t)/b(t)$ is the Hubble parameter associated with extra dimensions and $\gamma _D$ is normalized curvature of extra dimensions.

Equation that defines maximally symmetric solutions reads

$$\begin{aligned}&(D+3)(D+2)(D+1)D \alpha ^2H^4\nonumber \\&\quad +(D+3)(D+2)\alpha H^2-\xi =0, \end{aligned}$$

(5)

and cosmologically it corresponds to isotropic solution. It has real solutions iff

$$\begin{aligned} \xi \geqslant -\frac{(D+2)(D+3)}{4D(D+1)}. \end{aligned}$$

(6)

where $\xi =\alpha \Lambda$.

3 (3+2)-dimensional case with curvature

For $D=2$ system (4) takes form

$$\begin{aligned}&2\displaystyle \frac{(\gamma _D + \dot{b}(t)^2)}{b(t)^2} + 4\displaystyle \frac{\ddot{a}(t)}{a(t)} + 4\displaystyle \frac{\ddot{b}(t)}{b(t)} + 8\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)} + 2\displaystyle \frac{\dot{a}(t)^2}{a(t)}\nonumber \\&\quad + \alpha \left( 16\displaystyle \frac{\ddot{a}(t)(\gamma _D + \dot{b}(t)^2)}{a(t) b(t)^2} \right. \nonumber \\&\quad + \left. 32\displaystyle \frac{\ddot{a}(t) \dot{a}(t) \dot{b}(t)}{a(t)^2 b(t)} + 16\displaystyle \frac{\ddot{b}(t)\dot{a}(t)^2}{a(t)^2 b(t)} + 8\displaystyle \frac{\dot{a}(t)^2 (\gamma _D + \dot{b}(t)^2)}{a(t)^2 b(t)^2} \right. \nonumber \\&\quad + \left. 32\displaystyle \frac{\ddot{b}(t) \dot{a}(t) \dot{b}(t)}{a(t) b(t)^2} + 16\displaystyle \frac{\dot{a}(t)^2 \dot{b}(t)^2}{a(t)^2 b(t)^2} \right) = \Lambda , \nonumber \\&6\displaystyle \frac{\ddot{a}(t)}{a(t)} + 2\displaystyle \frac{\ddot{b}(t)}{b(t)} + 6\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)} + 6\displaystyle \frac{\dot{a}(t)^2}{a(t)} + \alpha \left( 24\displaystyle \frac{\ddot{a}(t) \dot{a}(t)^2}{a(t)^3} \right. \nonumber \\&\quad + \left. 48\displaystyle \frac{\ddot{a}(t) \dot{a}(t) \dot{b}(t)}{a(t)^2 b(t)} + 24\displaystyle \frac{\ddot{b}(t)\dot{a}(t)^2}{a(t)^2 b(t)} + 24\displaystyle \frac{\dot{a}(t)^3 \dot{b}(t)}{a(t)^3 b(t)} \right) = \Lambda , \nonumber \\ \end{aligned}$$

(7)

complimented with a constraint equation

$$\begin{aligned}&2\displaystyle \frac{(\gamma _D + \dot{b}(t)^2)}{b(t)^2} + 12\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)} + 6\displaystyle \frac{\dot{a}(t)^2}{a(t)}\nonumber \\&\quad + \alpha \left( 24\displaystyle \frac{\dot{a}(t)^2 (\gamma _D + \dot{b}(t)^2)}{a(t)^2 b(t)^2}\right. \nonumber \\&\quad \left. + 48\displaystyle \frac{\dot{a}(t)^3 \dot{b}(t)}{a(t)^3 b(t)} + 48 \displaystyle \frac{\dot{a}(t)^2 \dot{b}(t)^2}{a(t)^2 b(t)^2} \right) = \Lambda . \end{aligned}$$

(8)

We rewrote the system in terms of scale factors, as we are going to look for solutions with “stabilized extra dimensions” (so that the size of extra dimensions naturally becomes constant with respect to time). On the same time we require that three-dimensional subspace to expand with acceleration (after all, we are looking for a solution which could describe observed Universe), then, the conditions for such solution to exist are $a(t) = \exp (H_0t)$, $b(t) = b_0 \equiv \mathop {\mathrm{const}}\nolimits$, and the system (7)–(8) takes a form:

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{2\gamma _D}{b_0^2} + 6H_0^2 + \displaystyle \frac{24\gamma _D \alpha H_0^2}{b_0^2} = \Lambda , \;\; 12H_0^2 + 24\alpha H_0^4 = \Lambda . \end{array} \nonumber \\ \end{aligned}$$

(9)

One can see that three dynamical equations shrink to two – one of the variables becomes static ($b(t) = b_0 \equiv \mathop {\mathrm{const}}\nolimits$), so there is no “dynamical” equation which corresponds to it anymore. One can also note that we keep $\gamma _D$ arbitrary, unlike analysis for “geometric frustration” regime [4,5,6]. Now choosing new variables $x=1/b_0^2$ and $y=H_0^2$, we can rewrite (9) as

$$\begin{aligned} \begin{array}{l} 2z + 6y + 24\alpha zy = \Lambda , \quad 12y + 24\alpha y^2 = \Lambda , \end{array} \end{aligned}$$

(10)

where we absorbed $\gamma _D$ into x: $z = \gamma _D x$. Formally we now can get a direct solution: solve second of (10) w.r.t. y, substitute resulting $y_\pm$ into first of (10) and get $z_\pm$ for each branch. But to keep analysis consistent with the following sections, dedicated to higher-dimensional cases, where the equations will be cubic (for $(3+3)$) or even quartic (for higher-dimensional cases), we use same approach as we will be using further. Namely, we introduce $\xi = \alpha \Lambda$ and use scaling of the variables with respect to each other:

$$\begin{aligned} \begin{array}{l} \Lambda = \xi /\alpha , ~~~ b_0^2 = \zeta \alpha , ~~~ H_0^2 = \theta /\alpha . \end{array} \end{aligned}$$

(11)

With these redefinitions the system (10) takes a form

$$\begin{aligned} \begin{array}{l} 24\gamma _D\theta + 6\theta \zeta - \xi \zeta + 2\gamma _D =0 , \\ 12\theta + 24\theta ^2 = \xi . \end{array} \end{aligned}$$

(12)

We can immediately solve the second of it w.r.t. $\xi$:

$$\begin{aligned} \begin{array}{l} \xi = 12\theta (2\theta +1), \end{array} \end{aligned}$$

(13)

and substitute it into the first of (12) to get $\zeta$:

$$\begin{aligned} \begin{array}{l} \zeta = \displaystyle \frac{\gamma _D (12\theta +1)}{3\theta (4\theta +1)}. \end{array} \end{aligned}$$

(14)

Now the initial system is brought to just 1-parametric solution for $\xi$ (13) and $\zeta$ (14) with $\theta$ as a parameter. Let us investigate existence of the solutions for each particular combination of $\{\alpha ,\,\gamma _D\}$; the solutions are illustrated in Fig. 1.

Fig. 1

Graphs illustrating the behavior of derived functions $\xi$ (13) and $\zeta$ (14) for different cases in $(3+2)$-dimensional model: $\alpha > 0$, $\gamma _D > 0$ case: $\zeta (\theta )$ in a and $\xi (\theta )$ in b; $\alpha < 0$, $\gamma _D > 0$ case: $\zeta (\theta )$ in c and $\xi (\theta )$ in d; $\zeta (\theta )$ for $\alpha > 0$, $\gamma _D < 0$ case presented in e while $\zeta (\theta )$ for $\alpha < 0$, $\gamma _D < 0$ case presented in f (see the text for more details)

Case 1  $\alpha > 0$, $\gamma _D > 0$. From definitions (11) $\alpha > 0$ means $\theta > 0$ and $\zeta > 0$; from (14) one can see that it is always fulfilled (see Fig. 1a) and from (13) one can see that $\xi > 0$ for this case (see Fig. 1b).

Case 2  $\alpha < 0$, $\gamma _D > 0$. Now $\alpha < 0$ and it means $\theta < 0$ and $\zeta < 0$; from (14) one can see that it is fulfilled for $\theta \in (-\infty ; -1/4)\cup (-1/12; 0)$ (see Fig. 1c). From (13) one can see that $\xi > -3/2$ in this case (see Fig. 1d).

Case 3  $\alpha > 0$, $\gamma _D < 0$. As in Case 1, $\alpha > 0$ means $\theta > 0$ and $\zeta > 0$, but now (14) one can see that it is never fulfilled, which means there ar no solutions for this case (see Fig. 1e).

Case 4  $\alpha < 0$, $\gamma _D < 0$. In a way this case is “opposite” to the Case 2, as now $\theta < 0$ and $\zeta < 0$ are fulfilled for $\theta \in (-1/4; -1/12)$ (see Fig. 1f). The corresponding range for $\xi$: $\xi \in (-3/2; -5/6)$ and the graph is the same as in Case 2 (see Fig. 1d).

3.1 Linear stability of the solutions

To address the linear stability, we perturb the system (7) around the solution with stabilized extra dimensions ($a(t) = \exp (H_0t)$, $b(t) = b_0 \equiv \mathop {\mathrm{const}}\nolimits$). For simplicity, we rewrite the system (7) back in terms of Hubble parameter $H(t) = \dot{a}(t)/a(t)$ (additionally, this effectively diminish number of degrees of freedom by one which is crucial for our task), and we perturb the system around solution $H(t) = H_0 + \delta H(t)$, $b(t) = b_0 + \delta b(t)$ with $H_0$ and $b_0$ governed by (9). The resulting system of perturbed equations takes a form

$$\begin{aligned}&4b_0 (1+4\alpha H_0^2) \ddot{\delta }b(t) + 8b_0 H_0 (1+4\alpha H_0^2) {\dot{\delta }} b(t)\nonumber \\&\quad + 2b_0 (6H_0^2 - \Lambda ) \delta b(t) + ( 16\alpha \gamma _D + 4b_0^2 ) {\dot{\delta }} H(t) + \nonumber \\&\quad + 12H_0 (4\alpha \gamma _D + b_0^2) \delta H(t) = 0, \nonumber \\&2( 1+12\alpha H_0^2 ) \ddot{\delta }b(t) + 6H_0( 1+12\alpha H_0^2 ) {\dot{\delta }} b(t)\nonumber \\&\quad + (12H_0^2( 1+12\alpha H_0^2 ) - \Lambda ) \delta b(t) + \nonumber \\&\quad + 6b_0(1+4\alpha H_0^2) {\dot{\delta }} H(t) + 24H_0b_0(1+4\alpha H_0^2) \delta H(t) = 0. \nonumber \\ \end{aligned}$$

(15)

To find the solution of the system in the exponential form we transform it into “normal modes” with redefinition of the second derivative ${\dot{\delta }} b = \delta y$; then the system (15) could be replaced with

$$\begin{aligned} \begin{array}{l} \begin{pmatrix} {\dot{\delta }} y\\ {\dot{\delta }} H \\ {\dot{\delta }} b \end{pmatrix} = M \begin{pmatrix} \delta y\\ \delta H \\ \delta b \end{pmatrix} \end{array} \end{aligned}$$

(16)

with M being matrix made of corresponding coefficients. Then the solutions of the (16) could be written in the exponential form with the exponents being eigenvalues of M. Then, for solution to be stable, all these exponents should simultaneously have negative real parts. With use of (11), (13) and (14) these exponents could be rewritten in terms of only $\theta$ for each choice of $\alpha$ and $\gamma _D$.

Fig. 2

Graphs illustrating stability of the solutions for different cases in $(3+2)$-dimensional model: $\alpha > 0$, $\gamma _D > 0$ case in a; $\alpha < 0$, $\gamma _D > 0$ and $\alpha < 0$, $\gamma _D < 0$ cases in b. Different colors correspond to different (three) branches of the eigenvalues (see the text for more details)

Our analysis suggests that for all cases where solutions exist, they are unstable: for Case 1 ($\alpha > 0$, $\gamma _D > 0$) one of the exponents always has positive real part (see Fig. 2a); same situation is for Case 2 ($\alpha < 0$, $\gamma _D > 0$) and Case 4 ($\alpha < 0$, $\gamma _D < 0$), presented in Fig. 2b. For Case 3 ($\alpha > 0$, $\gamma _D < 0$) there are no solutions, as we obtained earlier.

Overall, we can see that in ($3+2$)-dimensional model with curved 2-dimensional extra-dimensional subspace, there are no stable solutions with stabilizing extra dimensions, regardless of the curvature of extra-dimensional part.

4 (3+3)-dimensional case with curvature

The system of dynamical equations (4) for $(3+3)$-dimensional case reads

$$\begin{aligned}&6\displaystyle \frac{(\gamma _D + \dot{b}(t)^2)}{b(t)^2} + 4\displaystyle \frac{\ddot{a}(t)}{a(t)} + 6\displaystyle \frac{\ddot{b}(t)}{b(t)} + 12\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)} + 2\displaystyle \frac{\dot{a}(t)^2}{a(t)}\nonumber \\&\quad + \alpha \left( 48\displaystyle \frac{\ddot{a}(t)(\gamma _D + \dot{b}(t)^2)}{a(t) b(t)^2} + \right. \nonumber \\&\quad + \left. 48\displaystyle \frac{\ddot{a}(t) \dot{a}(t) \dot{b}(t)}{a(t)^2 b(t)} + 24\displaystyle \frac{\ddot{b}(t)\dot{a}(t)^2}{a(t)^2 b(t)}\right. \nonumber \\&\quad \left. + 24\displaystyle \frac{\dot{a}(t)^2 (\gamma _D + \dot{b}(t)^2)}{a(t)^2 b(t)^2} + 24\displaystyle \frac{\ddot{b}(t) (\gamma _D + \dot{b}(t)^2)}{b(t)^3} \right. \nonumber \\&\quad + \left. 96\displaystyle \frac{\ddot{b}(t) \dot{a}(t) \dot{b}(t)}{a(t) b(t)^2} + 48\displaystyle \frac{\dot{a}(t) \dot{b}(t) (\gamma _D + \dot{b}(t)^2)}{a(t) b(t)^3}\right. \nonumber \\&\quad \left. + 48\displaystyle \frac{\dot{a}(t)^2 \dot{b}(t)^2}{a(t)^2 b(t)^2} \right) = \Lambda , \nonumber \\&2\displaystyle \frac{(\gamma _D + \dot{b}(t)^2)}{b(t)^2} + 6\displaystyle \frac{\ddot{a}(t)}{a(t)} + 4\displaystyle \frac{\ddot{b}(t)}{b(t)} + 12\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)}\nonumber \\&\quad + 6\displaystyle \frac{\dot{a}(t)^2}{a(t)} + \alpha \left( 24\displaystyle \frac{\ddot{a}(t) \dot{a}(t)^2}{a(t)^3} \right. \nonumber \\&\quad + \left. 24\displaystyle \frac{\ddot{a}(t)(\gamma _D + \dot{b}(t)^2)}{a(t) b(t)^2} + 96\displaystyle \frac{\ddot{a}(t) \dot{a}(t) \dot{b}(t)}{a(t)^2 b(t)}\right. \nonumber \\&\quad \left. + 48\displaystyle \frac{\ddot{b}(t)\dot{a}(t)^2}{a(t)^2 b(t)} + 24\displaystyle \frac{\dot{a}(t)^2 (\gamma _D + \dot{b}(t)^2)}{a(t)^2 b(t)^2} \right. \nonumber \\&\quad + \left. 48\displaystyle \frac{\dot{a}(t)^3 \dot{b}(t)}{a(t)^3 b(t)} + 48\displaystyle \frac{\ddot{b}(t) \dot{a}(t) \dot{b}(t)}{a(t) b(t)^2}\right. \nonumber \\&\quad \left. + 48\displaystyle \frac{\dot{a}(t)^2 \dot{b}(t)^2}{a(t)^2 b(t)^2} \right) = \Lambda , \end{aligned}$$

(17)

complimented with a constraint equation

$$\begin{aligned}&6\displaystyle \frac{(\gamma _D + \dot{b}(t)^2)}{b(t)^2} + 18\displaystyle \frac{\dot{a}(t) \dot{b}(t)}{a(t) b(t)} + 6\displaystyle \frac{\dot{a}(t)^2}{a(t)}\nonumber \\&\quad + \alpha \left( 72\displaystyle \frac{\dot{a}(t)^2 (\gamma _D + \dot{b}(t)^2)}{a(t)^2 b(t)^2} \right. \nonumber \\&\quad + \left. 72\displaystyle \frac{\dot{a}(t)^3 \dot{b}(t)}{a(t)^3 b(t)} + 72\displaystyle \frac{\dot{a}(t) \dot{b}(t) (\gamma _D + \dot{b}(t)^2)}{a(t) b(t)^3}\right. \nonumber \\&\quad \left. + 144 \displaystyle \frac{\dot{a}(t)^2 \dot{b}(t)^2}{a(t)^2 b(t)^2} \right) = \Lambda . \end{aligned}$$

(18)

The overall procedure is quite similar to the previously described $(3+2)$-dimensional case. The original system (17)–(18) could be brought to the following form under “stable compactification” requirement ($a(t) = \exp (H_0t)$, $b(t) = b_0 \equiv \mathop {\mathrm{const}}\nolimits$):

$$\begin{aligned}&\displaystyle \frac{6\gamma _D}{b_0^2} + 6H_0^2 + \displaystyle \frac{72\gamma _D \alpha H_0^2}{b_0^2} = \Lambda , \nonumber \\&\displaystyle \frac{2\gamma _D}{b_0^2} + 12H_0^2 + 24\alpha H_0^4 + \displaystyle \frac{48\gamma _D \alpha H_0^2}{b_0^2} = \Lambda ; \end{aligned}$$

(19)

where we kept $\gamma _D$ arbitrary. Choosing new variables $x=1/b_0^2$ and $y=H_0^2$, expressing one of new variables from first of (19) and substituting it into the second of (19), we can arrive to a pair of cubic equations:

$$\begin{aligned}&432\alpha ^2 y^3 + 180\alpha y^2 + (15-6\xi )y = \Lambda , \nonumber \\&864\alpha ^2 z^3 + 72\alpha (2\xi +5) z^2 + 30z = \Lambda (2\xi + 3); \end{aligned}$$

(20)

where we absorbed $\gamma _D$ into x: $z = \gamma _D x$ and used standard notation $\xi = \alpha \Lambda$. From the definitions used, we should have solutions within $x>0$, $y>0$, and the analysis in these coordinates becomes quite cumbersome (though not impossible). To simplify things we are going to use different approach – same as for $(3+2)$-dimensional case. Namely, we use the same redefinitions (11); substituting them into (19), the system takes a form

$$\begin{aligned}&72\gamma _D\theta + 6\theta \zeta - \xi \zeta + 6\gamma _D = 0; \nonumber \\&24\theta ^2\zeta + 48\gamma _D\theta + 12\theta \zeta - \zeta \xi + 2\gamma _D = 0. \end{aligned}$$

(21)

And this system has 1-parametric family of solutions:

$$\begin{aligned} \begin{array}{l} \xi = \displaystyle \frac{3\theta (144\theta ^2 + 60\theta + 5)}{6\theta + 1}, ~~\zeta = \displaystyle \frac{2}{3} \displaystyle \frac{\gamma _D (6\theta +1)}{\theta (4\theta +1)}. \end{array} \end{aligned}$$

(22)

So that given $\alpha$ and $\gamma _D$ we can build all possible solutions for all possible $\theta$. Let us analyze possible solutions and areas of their definitions. To start with, let us notice that from (11) it is clear that for $\alpha > 0$ we should have $\zeta > 0$ and $\theta > 0$ (as both $b_0^2$ and $H_0^2$ cannot be negative) while for $\alpha < 0$ we should have $\zeta < 0$ and $\theta < 0$. Let us consider all four possible combinations of signs for $\alpha$ and $\gamma _D$ separately.

Fig. 3

Graphs illustrating the behavior of derived functions (22) for different cases in $(3+3)$-dimensional model: $\alpha > 0$, $\gamma _D > 0$ case: $\zeta (\theta )$ in a and $\xi (\theta )$ in b; $\alpha < 0$, $\gamma _D > 0$ case: $\zeta (\theta )$ in c and $\xi (\theta )$ in d; $\zeta (\theta )$ for $\alpha > 0$, $\gamma _D < 0$ case presented in e while $\zeta (\theta )$ for $\alpha < 0$, $\gamma _D < 0$ case presented in f (see the text for more details)

Case 1  $\alpha > 0$, $\gamma _D > 0$. In that case $\theta > 0$ and substituting it into $\zeta$ we can see that $\zeta > 0$ for all $\theta > 0$ (see Fig. 3a). Substituting $\theta > 0$ into $\xi$ we can see that $\xi > 0$ (see Fig. 3b). So that for $\alpha > 0$ and $\xi > 0$ there always exist solution with $\gamma _D > 0$.

Case 2  $\alpha < 0$, $\gamma _D > 0$. In this case, since $\alpha < 0$, we should consider only $\theta < 0$ and the solution should have $\zeta < 0$, which is satisfied in two regions (see Fig. 3c): $\theta < -1/4$ and $-1/6< \theta < 0$. Comparing with Fig. 3d and performing some basic math, we can learn that within the first range ($\theta < -1/4$) we have $\xi > 0$ at $\theta < \theta _1$, and $\xi$ is negative at $\theta \in (\theta _1, -1/4)$, with $\xi (-1/4) = -3/2$. Within the second range ($-1/6< \theta < 0$), $\xi > 0$ for $\theta \in (-1/6, \theta _2)$ and is negative for $\theta \in (\theta _2, 0)$. Within the second range we can spot minimum value for $\xi$ which happening at $\theta _{min} = \root 3 \of {109 + 6\sqrt{330}}/72 + 1/(72\root 3 \of {109 + 6\sqrt{330}}) - 11/72 \approx -0.0669$; $\xi (\theta _{min}) \approx -0.5467$. All quoted specific values for $\theta$ are roots or singularities coming from (22); $\theta _1$ and $\theta _2$ ($\theta _1 < \theta _2$) are roots of $(144\theta ^2 + 60\theta + 5)=0$: $\theta _{1,\, 2} = (-5\pm \sqrt{5})/24 \approx \{-0.302,\,-0.115\}$.

To conclude, for $\alpha < 0$ we can have $\gamma _D > 0$ solutions for both $\xi < 0$ and $\xi > 0$.

Case 3  $\alpha > 0$, $\gamma _D < 0$. In this case we should have $\theta > 0$ and after plotting $\zeta (\theta )$ in Fig. 3e, we see that $\zeta < 0$ everywhere which means that there are no solutions of this kind.

Case 4  $\alpha < 0$, $\gamma _D < 0$. In this case we have $\theta < 0$, after plotting $\zeta (\theta )$ in Fig. 3f we see that it is “opposite” to $\alpha < 0$, $\gamma _D > 0$ case, again we should have $\zeta < 0$ and the area of existence is complimentary to $\alpha < 0$, $\gamma _D > 0$ case – now it is $\theta \in (-1/4, -1/6)$. Since $\xi (\theta )$ does not depend on $\gamma _D$, we can use Fig. 3d again to find that $\xi < -3/2$ in this case. So that $\gamma _D < 0$ solution could exist only for $\alpha < 0$ and they could be found within $\theta \in (-1/4, -1/6)$ range.

Now let us analyze linear stability of the obtained solutions.

4.1 Linear stability of the solutions

Similarly to the previous case, to study the linear stability, we perturb the system (17) around the solution with stabilized extra dimensions. Since the curvature of the internal subspace is zeroth, we can rewrite the system (17) in terms of Hubble parameter $H(t) = \dot{a}(t)/a(t)$ (effectively this diminish number of degrees of freedom by one which is crucial for our task), then we perturb the system around solution $H(t) = H_0 + \delta H(t)$, $b(t) = b_0 + \delta b(t)$ with $H_0$ and $b_0$ governed by (19) and pervious section describe a way how to find them. The resulting system of perturbed equations takes a form

$$\begin{aligned}&(24\alpha H_0^2 b_0^2 + 24\alpha \gamma _D + 6b_0^2) \ddot{\delta }b(t)\nonumber \\&\quad + (48\alpha H_0^3 b_0^2 + 48\alpha H_0\gamma _D + 12 H_0 b_0^2 ) {\dot{\delta }} b(t) \nonumber \\&\quad + ( 72 \alpha H_0^2 \gamma _D + 18 H_0^2 b_0^2 - 3\Lambda b_0^2 + 6\gamma _D ) \delta b(t)\nonumber \\&\quad + ( 48\alpha b_0\gamma _D + 4 b_0^3 ) {\dot{\delta }} H(t) \nonumber \\&\quad + ( 144\alpha H_0 b_0 \gamma _D + 12H_0 b_0^3 ) \delta H(t) = 0, \nonumber \\&( 48\alpha H_0^2 b_0 + 4b_0 ) \ddot{\delta }b(t) + ( 144\alpha H_0^3 b_0 + 12H_0 b_0 ) {\dot{\delta }} b(t)\nonumber \\&\quad + (48\alpha H_0^4 b_0 + 24H_0^2 b_0 - 2\Lambda b_0) \delta b(t) \nonumber \\&\quad + (24\alpha H_0^2 b_0^2 + 24\alpha \gamma _D + 6 b_0^2) {\dot{\delta }} H(t)\nonumber \\&\quad + (96\alpha H_0^3 b_0^2 + 96\alpha H_0\gamma _D + 24H_0 b_0^2) \delta H(t) = 0. \end{aligned}$$

(23)

We again use normal modes, and with use of (11) and (22) the exponents from the eigenvalues could be rewritten in terms of only $\theta$ for each choice of $\alpha$ and $\gamma _D$. So that we examine the exponents for all cases and present our analysis in Fig. 4.

Fig. 4

Graphs illustrating stability of the solutions for different cases in $(3+3)$-dimensional model: $\alpha > 0$, $\gamma _D > 0$ case in a panel; $\alpha < 0$, $\gamma _D > 0$ and $\alpha < 0$, $\gamma _D < 0$ cases in b (“large-scale”) and c (“fine structure”). Different colors correspond to different (three) branches of the eigenvalues (see the text for more details)

Case 1  $\alpha > 0$, $\gamma _D > 0$. Two out of three exponents are real and negative while the third is real and positive for all $\theta$, making this solution unstable (see Fig. 4a). We also remind the reader that for Case 3 there are no solutions.

Cases 2 and 4  $\alpha < 0$ for $\gamma _D > 0$ and $\gamma _D < 0$ respectively – have the same structure of stability areas so we report them together. First branch is stable for $\theta \in (-\infty , \theta _1) \cup (-1/4, \theta _4)$, second branch has negative real part of the exponent everywhere in $\theta < 0$ and the third within $\theta \in (-\infty , \theta _2) \cup (\theta _5, 0)$. Here $\theta _1$ and $\theta _2$ are the same as defined above while $\theta _3 < \theta _4$ are roots of $(55296\theta ^4 + 44352\theta ^3 + 12768\theta ^2 - 1480\theta + 55) = 0$ – radicand from the expression and $\theta _5 = -1/12$. Since we call solution as stable if all three exponents are negative, the overall stability range is the overlap of all three branches; the resulting range is $\theta \in (-\infty , \theta _1) \cup (-1/4, \theta _2) \cup (\theta _5, \theta _4)$. Comparing these ranges with those of existence, we can see that Case 4 solutions are stable on the entire area of definition ($\theta \in (-1/4, -1/6)$). On the other hand, Case 2 solutions overlap with stability range for $\theta \in (-\infty , \theta _1) \cup (-1/6, \theta _2) \cup (\theta _5, \theta _4)$. The situation is illustrated in Figs. 4b, c, where on (b) panel we presented “large-scale” picture while on (c) panel we focus on the range near zero, where there is “fine-structure” of the solutions.

Let us make a note on $\xi$, as it defines $\Lambda$: for Case 2 both ranges $\theta \in (-\infty , \theta _1)$ and $\theta \in (-1/6, \theta _2)$ give $\xi \in (0, +\infty )$ which, combined with $\alpha < 0$, gives us $\Lambda < 0$ – so that these Case 2 solutions exist for $\alpha < 0$, $\Lambda < 0$. The last range $\theta \in (\theta _5, \theta _4)$ gives $\xi \in (-0.5448, -0.5)$ – tiny but negative range, resulting in $\Lambda > 0$.

For Case 4 we have $\theta \in (-1/4, -1/6)$ which results in $\xi \in (-\infty , -3/2)$, so that $\Lambda > 0$.

5 (3+4)-dimensional solution

For $(3+4)$-dimensional case the resulting system for stabilized extra dimensions would be 4th order polynomial, making it even harder to solve explicitly then $(3+3)$-dimensional case. So that we use same technic as for $(3+3)$ dimensions, namely, we use (11) for the resulting system and obtain (after neglecting denominator):

$$\begin{aligned} \begin{array}{l} 144\gamma _D \theta \zeta + 6\theta \zeta ^2 - \xi \zeta ^2 + 24\gamma _D^2 + 12\gamma _D \zeta = 0, \\ 24\theta ^2 \zeta + 144\gamma _D \theta + 12\theta \zeta - \xi \zeta + 6\gamma _D = 0. \end{array} \end{aligned}$$

(24)

We solve the second of (24) with respect to $\xi$ and substitute it into the first of (24); since we use normalization $\gamma _D = \pm 1$, obviously $\gamma _D^2 = 1$, so the resulting equation takes a form

$$\begin{aligned} \begin{array}{l} 24\zeta ^2\theta ^2 + 6\zeta ^2\theta - 6\zeta \gamma _D - 24 = 0, \end{array} \end{aligned}$$

(25)

which has a solution

$$\begin{aligned} \begin{array}{l} \zeta _{\pm } = \displaystyle \frac{3\gamma _D \pm 3 |8\theta + 1|}{6\theta (4\theta + 1)}. \end{array} \end{aligned}$$

(26)

Now we can substitute it into previously found expression for $\xi$ from the second of (24):

$$\begin{aligned} \xi _{\pm }= & {} \displaystyle \frac{12\theta }{|8\theta +1|\pm \gamma _D} \big ( \pm 96\gamma _D \theta ^2 + 2\theta |8\theta +1| \pm 30\gamma _D \theta \nonumber \\&+ |8\theta +1| \pm 2\gamma _D \big ). \end{aligned}$$

(27)

Equations (26) and (27), together with definitions (11), completely determine two branches of 1-parametric solution with stabilized extra dimensions in $(3+4)$-dimensional model. Let us analyze when these solutions exist. As in $(3+3)$-dimensional case, we consider all possible combinations of $\alpha$ and $\gamma _D$ separately.

Fig. 5

Graphs illustrating the behavior of derived functions (22) for different cases in $(3+4)$-dimensional model: $\alpha > 0$, $\gamma _D > 0$ case: $\zeta _\pm (\theta )$ in a ($\zeta _+$ in red and $\zeta _-$ in blue) and $\xi _+(\theta )$ in b; $\alpha < 0$, $\gamma _D > 0$ case: $\zeta _\pm (\theta )$ in c ($\zeta _+$ in red and $\zeta _-$ in blue), $\xi _-(\theta )$ in d and $\xi _+(\theta )$ in e; $\alpha > 0$, $\gamma _D < 0$ case: $\zeta _\pm (\theta )$ in f panel ($\zeta _+$ in red and $\zeta _-$ in blue) and $\xi _+(\theta )$ in g; $\alpha < 0$, $\gamma _D < 0$ case: $\zeta _\pm (\theta )$ in h ($\zeta _+$ in red and $\zeta _-$ in blue) and $\xi _-(\theta )$ in i; (see the text for more details)

Case 1  $\alpha > 0$, $\gamma _D > 0$. In that case $\theta > 0$ and substituting it into $\zeta _\pm$ we can see that $\zeta _+ > 0$ everywhere in $\theta > 0$ (colored in red in Fig. 5a) while $\zeta _- < 0$ (colored in blue in Fig. 5a). So we can conclude that only $\zeta _+$ is viable and plot corresponding $\xi _+$ in Fig. 5b; from it one can see that $\xi _+ > 0$ which means $\Lambda > 0$ (since $\alpha > 0$).

Case 2  $\alpha < 0$, $\gamma _D > 0$. In that case $\theta < 0$ and substituting it into $\zeta _\pm$ we can see that $\zeta _- < 0$ everywhere in $\theta < 0$ (colored in blue in Fig. 5c) while $\zeta _+ < 0$ only for $\theta > -1/4$ (colored in red in Fig. 5c). For “–” branch we plot $\xi _-$ in Fig. 5d and one can see that it is always positive (and so that $\Lambda < 0$) while for “+” branch $\xi _+$ is plotted in Fig. 5e and we can see that $\xi$ could be both positive and negative there. Minimal possible value is $\xi _+(-1/4) = -3/2$; $\xi _+(\theta )$ hits zero at $\theta _1 = -(5+\sqrt{5})/40 \approx -0.1809$ and $\theta _2 = -3/28 \approx -0.1071$ (and at $\theta = 0$); it has maximum at $\theta =-1/8$ with $\xi _+(-1/8) = 3/8$ and finally it has local minimum at $\theta _{min} = -3/56 \approx -0.05357$ with $\xi _+(\theta _{min}) = -27/56 \approx -0.48214$. So that this case has a variety of parameters where solutions could possibly exist.

Case 3  $\alpha > 0$, $\gamma _D < 0$. In that case $\theta > 0$ and substituting it into $\zeta _\pm$ we can see that $\zeta _+ > 0$ everywhere in $\theta > 0$ (colored in red in Fig. 5f) while $\zeta _- < 0$ (colored in blue in Fig. 5f). So we can conclude that only $\zeta _+$ is viable and plot corresponding $\xi _+$ in Fig. 5g; from it one can see that $\xi _+ < 0$ which means $\Lambda < 0$ (since $\alpha > 0$). One cannot miss familiarity between Cases 1 and 3 – “mirror-like” behavior of $\zeta _\pm$ branches (compare Fig. 5a, f) but since they have different sign for $\gamma _D$, the resulting $\xi$ also have different sign (compare Fig. 5b, g).

Case 4  $\alpha < 0$, $\gamma _D < 0$. In that case $\theta < 0$ and substituting it into $\zeta _\pm$ we can see that $\zeta _+ > 0$ everywhere in $\theta < 0$ (colored in red in Fig. 5h) while $\zeta _- < 0$ only for $\theta < -1/4$ (colored in blue in Fig. 5h). So that the only viable branch is $\zeta _-$ and only for $\theta < -1/4$. The resulting $\xi _-(\theta )$ graph is presented in Fig. 5i – one can clearly see that $\xi < -3/2$ and so $\Lambda > -3/(2\alpha )$. Again, one cannot miss familiarity between Cases 2 and 4 – their $\zeta _\pm$ curves are “mirror-symmetric” with respect to $\zeta =0$ (compare Fig. 5c, h) but unlike Cases 1 and 3 described above, Cases 2 and 4 have more complicated structure which results in quite different branch/range combination for each of them. Also, since now we have only one branch, for $\xi$ we also have only single choice.

Overall, one can see much more abundant potential dynamics then for $(3+3)$-dimensional case; let us find out if these solutions are stable or not.

5.1 Linear stability of (3+4)-dimensional solutions

The general scheme for finding stable solutions is exactly the same as in $(3+3)$-dimensional case – we perturb the system of dynamical equations around the solution, separate perturbations, find normal modes and investigate when they all simultaneously have negative real parts. Skipping technical details, we report the results for all cases separately; the results are presented in Fig. 6.

Fig. 6

Graphs illustrating stability of the solutions for different cases in $(3+4)$-dimensional model: $\alpha > 0$, $\gamma _D < 0$ case in a; $\alpha < 0$, $\gamma _D > 0$ with $\zeta _-$ branch in b (“large-scale”) and c (“fine structure”) panels; $\alpha < 0$, $\gamma _D > 0$ with $\zeta _+$ branch as well as $\alpha < 0$, $\gamma _D < 0$ with $\zeta _-$ branch in d (“fine structure”). Different colors correspond to different (three) branches of the eigenvalues (see the text for more details)

Case 1  In this case one of the modes is always real and positive while two others are always real and negative – overall it means that the solutions are unstable for this case. The situation resemble same Case 1 from $(3+3)$-dimensional case so Fig. 4a would serve as a good illustration.

Case 2  As previously described, here we have possible solutions on both $\zeta _\pm$ branches. For $\zeta _-$ branch one of the modes always has positive real part which makes this branch unstable. Situation is illustrated in Fig. 6b, c, with (b) panel presents “large-scale“ picture and (c) panel – “fine-tuning” at small values of $\theta$ where the structure is complicated. On the contrary, $\zeta _+$ has all three modes with negative real parts at $\theta \in (-\infty , -1/4) \cup (-1/5, -1/8) \cup (-1/16, -3/56)$ (see Fig. 6d). If we unite it with the region of existence described earlier, we will have to drop the first range, resulting in $\theta \in (-1/5, -1/8) \cup (-1/16, -3/56)$ for range for existence and stability of Case 2 solution.

Case 3  Here we have $\zeta _+$ and all three modes are real and negative everywhere on $\theta > 0$, making entire range of definition for solutions within this class stable (see Fig. 6a).

Case 4  Similarly to the previous case, all three modes have negative real parts everywhere in the domain of definition (now, for $\theta < -1/4$) making them stable (see Fig. 6d). It is interesting to note that $\zeta _+$ for Case 2 and $\zeta _-$ for Case 4 have exactly the same exponents in the stability analysis – it seems that sign swaps for $\alpha$ and $\gamma _D$ exactly cancel each other.

To conclude, solutions with $\gamma _D > 0$ exist and stable only for $\alpha < 0$ (Case 2), for $\theta \in (-1/5, -1/8) \cup (-1/16, -3/56)$, which corresponds to $\xi \in (-27/56, -15/32) \cup (-0.3, 3/8)$ (so that for both positive and negative $\Lambda$). Solutions with $\gamma _D < 0$ exist and stable for both $\alpha > 0$, $\Lambda > 0$ (Case 3) and $\alpha < 0$, $\Lambda > 0$ ($\xi < -3/2$) (Case 4).

6 General $(3+D)$-dimensional case

The procedure for the general case is exactly the same as for previously described $(3+3)$- and $(3+4)$-dimensional cases. So with use of (11), the resulting system for stabilized extra dimensions reads

$$\begin{aligned}&D(D-1)(D-2)(D-3)\gamma _D^2 + D\zeta \gamma _D (12\theta +1)(D-1)\nonumber \\&\quad + \zeta ^2 (6\theta - \xi ) = 0; \nonumber \\&D(D-1)(D-2)(D-3)(D-4) \gamma _D^2\nonumber \\&\quad + \zeta \gamma _D (D-1)(D-2)(24\theta +1)\nonumber \\&\quad + \zeta ^2(24\theta ^2 + 12\theta - \xi ) = 0. \end{aligned}$$

(28)

One can immediately see that due to the multiplier $(D-4)$ in the second equation, we cannot obtain $D=4$ case as a subcase of general one. Following the procedure, we express $\xi$ from one of the equations and substitute it into another one. The resulting equation is quadratic with respect to $\zeta$ but, unlike previous $(3+3)$- and $(3+4)$-dimensional cases its solutions need radicals to be written explicitly, nevertheless, they still could be written in a closed form:

$$\begin{aligned} \xi= & {} \displaystyle \frac{12\theta (1+2\theta )\zeta ^2 + \gamma _D(D-1)(D-2)(24\theta +1)\zeta + \gamma _D^2(D-1)(D-2)(D-3)(D-4)}{\zeta ^2}, \nonumber \\ \zeta _\pm= & {} \displaystyle \frac{-\gamma _D (D-1) (6D\theta - 24\theta - 1 \pm \sqrt{\mathcal {D}}) }{6\theta (4\theta +1)},~\text{ where } \nonumber \\ \mathcal {D}= & {} (D-1) \left( (D-1)\gamma _D^2 (6D\theta - 24\theta - 1)^2+ 24\theta (D-2)(D-3)(4\theta +1) \right) . \end{aligned}$$

(29)

Similarly the the previous cases, we consider separately four cases with different signs for $\alpha$ and $\gamma _D$. But unlike previous cases, now everything depends not only on our parameter $\theta$, but also on D – number of extra dimensions. In the previous cases it was easy to plot the resulting curves and find roots/divergences, now it is a bit more complicated.

Case 1 ($\alpha > 0$, $\gamma _D > 0$) and Case 3 ($\alpha > 0$, $\gamma _D < 0$) are quite similar: in both we have $\alpha > 0$ so that we can consider only $\theta > 0$ and for solution to exist we should have $\zeta > 0$. There, we can see from (29) that $\mathcal {D}^2 > (D+1)^2(6D\theta - 24\theta - 1)^2$ which makes $\zeta _+ > 0$ and $\zeta _- < 0$ for both signs of $\gamma _D$. So that for cases 1 and 3 we have a solution for the entire $\theta > 0$ for $\zeta _+$. When transforming to $\xi$, it corresponds to $\xi > 0$ for Case 1 and $\xi < -D(D-1)/(4(D-2)(D-3))$ for Case 3; one can see that the latter is different from $(3+4)$-dimensional case while the formed is the same.

Similar analysis could be used to Case 2 ($\alpha < 0$, $\gamma _D > 0$) and Case 4 ($\alpha < 0$, $\gamma _D < 0$). Now we have $\alpha < 0$ so that we consider only $\theta < 0$ and should have $\zeta < 0$ for the solution to exist. Again, comparing discriminant and the free term, we can conclude that Case 4 solutions exist for $\zeta _-$ and $\theta < -1/4$ while Case 2 solutions exist for all $\theta < 0$ on $\zeta _-$ and $0> \theta > -1/4$ on $\zeta _+$. One can see that existence regions for these two cases are exactly the same as for $(3+4)$-dimensional case.

6.1 Linear stability of general $(3+D)$-dimensional solutions

In the previously described cases for given $\alpha$ and $\gamma _D$ we could build a solution for any $\theta$ from its region of existence. This general case is different – $\theta$ is still our parameter but we additionally have number of extra dimensions D; also, solutions for $\zeta _\pm$ (29) no more have simple form – they have radicals which further complicate the analysis. Still, we can formally follow all the steps – perturb the equations of motion around the solution, substitute (11), calculate eigenvalues of the matrix for perturbation equations and the result will depend on $\theta$ with D and with $\alpha$ and $\gamma _D$ as parameters. As the resulting eigenvalues being a roots of cubic equation, their implicit analysis is quite troublesome, nevertheless one can demonstrate that they do not have extremum as a function of D for $D > 4$. With this knowledge at hand we investigate several particular D cases (we used $D = 5,\, 7,\, 10,\, 12$) and, seing no qualitative difference, decide that this is common behavior for general D. This way we demonstrate that in Case 1 one of the modes is always positive while two others always negative, making solutions from this case unstable – the same situation as in previous cases; Fig. 4a would serve as a good illustration for this case. On the other hand, solutions for Case 3 are stable, as all three modes are negative for $\theta > 0$ and $\zeta _+$. This is the same as in $(3+4)$-dimensional case, so that Fig. 6a is a good illustration. Case 4 also proves to be stable for $\zeta _-$ and $\theta < -1/4$ so that Fig. 6d is an example. Finally, Case 2 is the most complicated of all four cases - well, just like in $(3+3)$- and $(3+4)$-dimensional cases as well. For $\zeta _-$ branch, one of the modes always has positive real part, making it unstable – just like in $(3+4)$-dimensional case (see Fig. 6b, c). For $\zeta _+$ branch the situation is as follows: the range of definition is $0> \theta > -1/4$; out of three modes, one has negative real part for $\theta _1< \theta < \theta _4$, another mode has negative real part for $\theta \in (-1/4, \theta _2) \cup (\theta _3, 0)$ and final mode always negative. Here $-1/4< \theta _1< \theta _2< \theta _3< \theta _4 < 0$ are defined as follows: $\theta _1$ is very lengthy radicand coming from denominator of expression for mode; $\theta _3 = - 1/(4D)$ is another zero of the same denominator; $\theta _2$ and $\theta _4$ are roots of

$$\begin{aligned}&48D^3 (5D-6) (D+1) (3D^2 - 19D + 32)\theta ^4\nonumber \\&\quad + 24D^2(31D^4 - 185D^3 + 266D^2 - 272D - 480)\theta ^3 \nonumber \\&\quad + 12D (21D^4 - 103D^3 + 108D^2 - 268D - 384) \theta ^2\nonumber \\&\quad + 2(D-2)(2D-3)(9D^2 - 32D + 48)\theta \nonumber \\&\quad + (D-1)(D+2)(2D-3) = 0. \end{aligned}$$

(30)

For a solution to be stable, all three modes should be negative; so the resulting range is intersection of all three; since $-1/4< \theta _1< \theta _2< \theta _3< \theta _4 < 0$, the resulting range is $\theta \in (\theta _1, \theta _2)\cup (\theta _3, \theta _4)$ – see Fig. 6d as an example. Of course, all of these $\theta$’s scales with D and in Fig. 7 we presented them for some values of D – the ranges for $\theta$ are presented in (a) panel; if we convert these ranges into $\xi$, the resulting areas are presented in Fig. 7b. One can see from Fig. 7 that the area is decreasing with growth of D, making smaller the measure of the parameters when solutions in Case 2 could exist.

Fig. 7

Behavior of allowed regions for Case 2 solutions with varying D: regions for $\theta$ in a and for $\xi$ in b (see the text for more details)

7 Summary

We realize that the manuscript so far was quite technical and decided to summarize relevant results in a separate section. Before turning to summarizing the results, there is one more important note we want to make. Through this paper we report existence and stability regions mainly using $\theta$, and it is done for the reason. Indeed, as we described earlier, $\theta$ is the only “dynamical” parameter of the theory (“dynamical” here means that this parameter vary from one exact solution to another with $\alpha$, $\gamma _D$ and D remain unchanged). If we have a look on graphs of $\xi$ and $\zeta$ as a functions of $\theta$ we will see that the same value of $\xi$ and $\zeta$ often achieved at different values of $\theta$: this corresponds to the situation when there are multiple solutions for the same $\alpha$ and $\xi$. For instance, for $(3+3)$-dimensional case there could be up to three (and for general $(3+D)$ up to four) distinct solutions (with different $H_0$ and $b_0$ each) for a given $\alpha$ and $\xi$. And there will be distinct $\theta$ which corresponds to each of them. So that if we are interested in finding all possible solutions, we are interested in range of $\theta$; if we are interested in when at least one solution exist, range of $\xi$ would suffice.

With these notes taken, let us move to summarizing the results.

7.1 Negative curvature of extra dimensions

In fact, this case was studied quite well in [4, 5] and here it is just one of two subcases of the general $\gamma _D$. Our analysis suggest that the solutions exist in $(3+2)$- and $(3+3)$-dimensional cases only for $\alpha < 0$ while for $(3+4)$- and higher-dimensional cases – for both signs of $\alpha$. For $(3+2)$-dimensional case these solutions are unstable (see Fig. 2b) while for higher dimensions the situation becomes more complicated. The solutions with $\alpha < 0$ exist and stable within some ranges of $\theta$: for $(3+3)$ it is $\theta \in (-1/4, -1/6)$ which corresponds to $\xi < -3/2$ (see Fig. 4c); for $(3+4)$ it is $\theta < -1/4$ (and again $\xi < -3/2$); exactly the same result we obtain for general $(3+D)$ case. Let us note that despite limited range (in $(3+3)$-dimensional case it is even limited from both above and below), it is the entire range of where solution exist. So that we can conclude that for $\alpha < 0$ all existing solutions are stable.

Now let us turn to $\alpha > 0$ solutions – they exist only starting from $(3+4)$ and in higher-dimensional cases. Our analysis proves that $(3+4)$- and general $(3+D)$-dimensional cases have exactly the same description: it is $\zeta _+$ branch which exists (see Fig. 5f) and is stable (see Fig. 6a) everywhere within range of definition. So that not only for $\alpha < 0$, but for $\alpha > 0$ all existing solutions are stable as well. So that we can conclude that all existing solutions with negative curvature and $D>2$ are always stable.

7.2 Positive curvature of extra dimensions

Let us consider separately cases $\alpha > 0$ and $\alpha < 0$, as we did in the negative curvature section. The solutions for $\alpha > 0$ always exist (starting from $(3+2)$ and in any higher dimensions) but are always unstable – it is what we called Case 1 and one of the nodes always has positive real part in all D cases. On the contrary, starting from $(3+3)$-dimensional case, $\alpha < 0$ solutions have some range of stability: it is $\theta \in (-\infty , \theta _1) \cup (-1/6, \theta _2) \cup (\theta _5, \theta _4)$, and all $\theta$’s are defined in section dedicated to $(3+3)$-dimensional model; let us note that the entire $\xi > 0$ is covered within this range of $\theta$. For $(3+4)$-dimensional case the situation is as follows: there are two branches, $\zeta _\pm$, of them $\zeta _-$ defined everywhere at $\theta < 0$ (see Fig. 5c) but is stable nowhere (see Fig. 6b, c) while $\zeta _+$ is defined at $0> \theta > -1/4$ and is stable within $\theta \in (-1/5, -1/8) \cup (-1/16, -3/56)$, which corresponds to $\xi \in (-27/56, -15/32) \cup (-0.3, 3/8)$. Let us note the difference between this and $(3+3)$-dimensional cases – in the latter stability range is $\xi > 0$ – it is unbounded but single-signed. Finally general $(3+D)$-dimensional case follow $(3+4)$-dimensional: for $\alpha < 0$ we have unstable $\zeta _-$ branch and stability of $\zeta _+$ branch within two intervals; dependence of these intervals on D is depicted in Fig. 7a while corresponding range for $\xi$ – in Fig. 7b.

7.3 Concluding remarks

We summarize all existing and stability ranges in Table 1. Please note that we united results for different branches (when apply), as we are interested in conditions when solutions exist and are stable, leaving technical details to the main paper. From Table 1 we can clearly see the difference between the cases with positive and negative curvature: solutions with negative curvature are always stable while solutions with positive curvature are stable only within some range of $\xi$ and the size of this range (“measure of stable trajectories”) strongly decreases with growth of D, as indicated in Fig. 7b. The only exceptions are $(3+3)$-dimensional case, for which solutions are stable for all $(\alpha< 0, \Lambda < 0)$ and $(3+4)$-dimensional case, where solutions are stable for $\alpha < 0$, $\xi \in (-27/56, -15/32) \cup (-0.3, 3/8)$.

This can explain why we have not detected this sort of behavior back in [4, 5] – we studied the behavior numerically, and were interested in large-D asymptotic – in this case the behavior of models with negative curvature strongly favors detectability (the always stable – regardless of the parameters) while to find stable solution with positive curvature we need to fine-tune $\xi$, and with growth of D the size of range for $\xi$ decreases, making it practically impossible to obtain “correct” value from this range via random pick.

Table 1 Conditions for existing and stability of compactified solutions in different number of extra dimensions D

Finally, let us attend one more interesting feature of the obtained solutions. Beforehand, studying the case with negative curvature in [4, 5] we noticed that maximally-symmetric solutions do not coexist with compactified ones. In the case of positive curvature they do. This statement is important so we want to elaborate a bit. If we add criterium for maximally-symmetric solution (6) to stability regions on Fig. 7b, we obtain situation illustrated in Fig. 8. In there we took Fig. 7b and added exact existence separatrix for maximally-symmetric solution (6) as red dashed line. Since existence criteria for maximally-symmetric solutions is $\xi \geqslant \xi _{iso}$, then one can see that both regions satisfy it. Unfortunately, the expressions for three out of four $\xi (D)$ curves are quite cumbersome, so we are left with only one curve – $\xi _3(D)$ (which is the upper boundary of the bottom (smaller) strip). It corresponds to $\theta _3 = -1/(4D)$ (see appropriate section) and being substituted to $\xi$ in Eq. (29) (alongside with $\zeta _+$ from the same equation) gives us $\xi _3 = -(D^2+2D-9)/(4D(D-2))$. If we compare it with $\xi _{iso}$ from (6) we can see that

$$\begin{aligned} \begin{array}{l} \xi _3 - \xi _{iso} = \displaystyle \frac{3(D-1)}{4D(D+1)(D-2)} > 0, \end{array} \end{aligned}$$

(31)

so that $\xi _3 > \xi _{iso}$ always and we can claim that for $D>2$ there always exist stable compactified solutions with positive curvature which coexist with maximally-symmetric solutions – this situation is different from the case with negative curvature of extra dimensions. Figure 8 suggests that the statement is even stronger – that all stable compactified solutions with positive curvature coexist with maximally-symmetric solutions. To check this we can numerically calculate both $\xi _4$ (lower boundary of the bottom strip) for large values of D and compare it with $\xi _{iso}$ to verify that $\xi _{iso} < \xi _4$. For instance, for $D=100$ we have $\xi _4(D=100) = -0.2599803$ while $\xi _{iso}(D=100) = -0.2600495$, which support our statement in its stronger form; higher values for D support it as well.

Fig. 8

Allowed regions for $\xi$ as a functions of the number of extra dimensions (in blue; same as Fig. 7b) and separatrix for existence of maximally-symmetric solutions (red dashed line) (see the text for more details)

8 Conclusions

In this paper we have considered a compactification scenario where stabilisation of extra dimensions occurs due to presence the Gauss-Bonnet term and non-zero spatial curvature. The sign of spatial curvature can be positive or negative, in the latter case we need additional factorisation to get a compact inner space. Scale factor of the extra dimension space as well as the effective cosmological constant in three-dimensional “big” submanifold are obtained from the system of two algebraic equations. Roots of these equations depend upon coupling constants of the theory, and some combinations of coupling constants can lead either to absence of roots or to negative roots without an adequate physical meaning.

The first goal of the present paper is to describe combinations of coupling constants which allows a physically acceptable compactification solutions. We show that a product of constants $\xi =\alpha \Lambda$ characterises the existence conditions completely if supplemented by the information about signs of $\alpha$ and the spatial curvature. From our results summarised in Table 1 we can see that the existence conditions for both signs of the spatial curvature are not very restrictive. We note also that in the case of negative spatial curvature the compactification solution exist only in those range of $\xi$ where the maximally symmetric solution does not exist. This proves the “geometrical frustration” hypothesis for negatively curved inner spaces (see [4] for details). We can note that for $\alpha >0$ the compactification solution exists exactly for the range of $\xi$ where maximally symmetric solution does not exist, while for negative $\alpha$ there is a range of $\xi$ where neither compactification nor maximally symmetric solution exist. However, compactification solutions with positively curved inner spaces can co-exist with the maximally symmetric solution.

The second goal of our paper is to study stability of compactification solutions with respect to homogeneous perturbations of the metric. The first results obtained are the same for both signs of the spatial curvature – any compactification solution with only $D=2$ extra dimensions is unstable. This explains why such solutions have never been found in actual numerical integrations of equations of motion. The situation with bigger number of extra dimensions is, however, quite different for negatively and positively curved extra spaces. For the negative curvature extra space, any solution with $D>2$ extra dimensions is stable. Note also, that for $D>3$ both signs of the constant $\alpha$ are possible for the solution to exist (and, than, to be stable).

For the positive curvature case, on the contrary, stability conditions impose severe restrictions for possible set of the coupling constant of the theory under study. First of all, there are no stable compactification solutions with a positive $\alpha$. Second, for negative $\alpha$ a solution is stable only in rather narrow interval of $\xi$, and the width of this interval decreases with increasing D (see Fig. 7b).

Let us also note that for a flat internal space ($\gamma _D = 0$) the solutions under consideration do exist for $\alpha < 0$ and $\Lambda = - 3/(2 \alpha ) > 0$ and are stable for $D > 2$ [61] and $D=2$ [62].

We find that the range of $\xi$ which allows stable compactification solutions is located in zone where maximally symmetric solution also exist. The actual fate of a particular trajectory may depend upon the initial conditions, and this question needs further investigations. In general, positive curvature case seems to be more physically relevant since it leads to compactness of the inner space directly, while addition factorisation is needed for the negative curvature case. Our results indicate, however, that this advantage is in some sense compensated by the fine-tuning of coupling constants needed for the compactification solution to be stable. Since this fine-tuning is rather serious for the case of large number of dimensions in the inner space, further Lovelock terms might change these conclusions, and this is the matter of a separate investigation which we plan to provide in future.