Nonsingular cosmological models with strong gravity in the past

In scalar-tensor Horndeski theories, nonsingular cosmological models - bounce and genesis - are problematic because of potential ghost and/or gradient instabilities. One way to get around this obstacle is to send the effective Planck mass to zero in the asymptotic past ("strong gravity in the past"). One may suspect that this feature is a signal of a strong coupling problem at early times. However, the classical treatment of the cosmological background is legitimate, provided that the strong coupling energy scale remains at all times much higher than the scale associated with the classical evolution. We construct various models of this sort, namely (i) bouncing Universe which proceeds through inflationary epoch to kination (expansion within general relativity, driven by massless scalar field); (ii) bouncing Universe with kination stage immediately after bounce; (iii) combination of genesis and bounce, with the Universe starting from flat space-time, then contracting and bouncing to the expansion epoch; (iv)"standard"genesis evading the strong coupling problem in the past. All these models are stable, and perturbations about the backgrounds are not superluminal.


Introduction
Nonsingular cosmological models bouncing cosmology and genesis from Minkowski space are of continuous interest as alternatives to or completions of inflation. Provided that the spatial curvature is negligible, a prerequisite for the construction of these models is the stable violation of the null energy condition (and, more generally, null convergence condition). It is known since 2010 [1,2,3] that the latter feature can exist in Horndeski theories [4] (for reviews see, e.g., Refs. [5,6]). These are scalar-tensor modifications of gravity, with the Lagrangians containing second derivatives of both the metric and scalar field, and yet with the second-order equations of motion. Indeed, within Horndeski theories, numerous explicit examples of stable early genesis [1,7,8,9,10,11,12] and bouncing [13,14,15,16,17,15,18,19,20,21] stages were constructed.
However, within Horndeski theory, these cosmologies typically suffer from either singularities or gradient and/or ghost instabilities at some earlier or later stage (possibly well after the initial genesis epoch and, likewise, well before or well after the bounce). In earlier papers, this property was observed explicitly in most cases where the evolution was followed by sufficiently distant past and future [10,12,15,18,19,20,21] (see Ref. [22] for the discussion of other problematic properties of the model of Ref. [21]). Later on, the problem has been formulated as a no-go theorem [23,24]. Namely, in the unitary gauge and in the spatially flat Friedmann-Lemaître-Robertson-Walker background ds 2 = dt 2 − a 2 (t)δ ij dx i dx j , the quadratic actions for tensor (transverse traceless) perturbation h ij and scalar perturbation ζ have the forms The theorem states that if in a Horndeski theory the background is nonsingular during the entire evolution −∞ < t < +∞, the coefficient G T is strictly positive at all times, and the following two integrals are divergent at lower and upper limits, respectively, then F S < 0 and/or F T < 0 in some time interval, i.e., there exists either ghost or gradient instability (or both).
As a digression, we emphasize that like most cosmology model builders, we stick to the study of homogeneous and isotropic backgrounds and their stability against linearized perturbations. Like most others, we are confident that the standard (3+1) decomposition with algebraic gauge conditions (unitary gauge in our case) is adequate for this particular purpose: the wave equations derived from the actions (1) are manifestly strongly hyperbolic for positive G T , F T , G S , F S and manifestly elliptic for negative sound speeds squared (but we tend to agree with Ref. [27] that the algebraic gauges, including unitary, may not be convenient for analyzing the evolution at a fully nonlinear level).
We leave aside the issue of stability at the nonlinear level and, even more so, the issue of well posedness of general backgrounds in Horndeski theories; the latter issues are discussed, e.g., in Refs. [27,28,29,30]. In this regard, positivity of G T , F T , G S , and F S is necessary, albeit possibly not a sufficient condition for a healthy cosmological model.
One way to deal with instabilities implied by the no-go theorem is to arrange for (or merely declare) a sufficiently low energy scale of the UV completion and make sure that the unstable modes (with energies below this scale) do not have enough time to develop [10,31,32,33]. Another is to get around these instabilities altogether [34,35,36,37,38] by making use of beyond Horndeski [39,40] or more general degenerate higher-order scalar-tensor theories [41,42], which, however, have problems with superluminality [43]. In this paper we follow yet another route [24], namely, we stick to the Horndeski theory and ensure that the coefficients G T , F T , G S , and F S ("effective Planck masses" squared) in the quadratic actions for perturbations (1) sufficiently rapidly decay as one goes backwards in time to t → −∞, so that the integral in the left-hand side of Eq. (2a) is actually convergent. In this way we relax the assumption of the no-go theorem and construct Horndeski models without gradient or ghost instabilities (at the linearized level). A peculiarity of this case is that the gravitational and scalar interactions are strong at early times, which, among other things, signalizes a potentially strong coupling problem. 1 For brevity, we refer to this property as "strong gravity in the past." In the latter class of models, the fact that the effective Planck masses tend to zero as t → −∞ does not necessarily mean that the classical field theory treatment of the (homogeneous and isotropic) cosmological evolution is not legitimate at early times [45,46]. Indeed, the classical analysis is valid, provided that the quantum strong coupling energy scale E strong stays well above the energy scale of the classical evolution E class (for power-law evolution, the latter is E class ∼ |t| −1 as t → −∞). This issue was considered in Refs. [45,46,47] in the framework of the class of models suggested in Ref. [24]. Using the dimensional analysis, it was shown, order by order in perturbation theory, that there actually exists a region in the parameter space where the classical field theory treatment is legitimate. 2 It is worth noting, though, that the parameters of the concrete model given in Ref. [24] do not belong to this region.
Interestingly, the result of the all-order analysis of Ref. [47] coincides with the result of Ref. [46] obtained by studying the cubic order only: both lead to the same constraints on the parameters. 3 1 There exists even more radical proposal that the effective Planck masses squared G T and F T vanish at some finite time t 0 , i.e., [44]. It remains to be seen whether or not models of this sort are tractable within classical field theory. 2 We note in passing that another model with vector field and power-law background solution was constructed in Ref. [48]; it describes stable early genesis which is legitimately treated within classical field theory. 3 A would-be caveat in the analysis of Ref. [47] is that it did not give explicit comparison of strong coupling energy scales emerging at different orders of perturbation theory. The subtlety is that when the strong coupling scale inferred from a higher order term is below that Explicit Horndeski models with strong gravity in the past have not been constructed so far. It is the purpose of this paper to fill this gap: we introduce several Horndeski cosmologies of this sort, which are stable at all times; we emphasize that we always work in the Jordan frame. We ensure that these models are free of the strong coupling problem, i.e., E class E strong at all times, even though E class → 0 as t → −∞. Our cosmologies are complete in the sense that at late times the Universe expands in a standard way: at large positive t, the models turn into general relativity with a conventional massless scalar field that drives the expansion. This is kination epoch which is assumed to end up with reheating through, say, one of the mechanisms of Refs. [50,51]. The least straightforward part of our construction is to ensure the (linear) stability of the solutions during the entire evolution. We also make sure that the speed of the perturbations about our backgrounds does not exceed the speed of light. So, our cosmologies are exotic but healthy (modulo possible pathologies at nonlinear level).
The first model, elaborated in greater detail, is the bouncing Universe. In the asymptotic past the Universe contracts with the power-law behavior of the scale factor, then the contraction terminates and expansion begins (bounce).
Depending on the choice of the Lagrangian, the expansion epoch may or may not pass through the inflationary stage; we give examples of both scenarios. As described above, we follow the evolution up to the kination epoch.
Another model is a combination of genesis and bounce: the Universe starts from the flat space-time, then contracts, passes through the bounce and then evolves in the same way as in the first example. For completeness, we design yet another "standard" genesis model (in which the Universe expands from the beginning), which satisfies the condition [47] of the absence of strong coupling and thus improves on the model of Ref. [24].
This paper is organized as follows. We introduce our subclass of models from the Horndeski class in Sec. 2, where we also discuss general properties of these models. Bouncing Universes are constructed in Sec. 3. Models with genesis are presented in Sec. 4. We conclude in Sec. 5. In Appendix A we derive the condition of the absence of strong coupling at early times in the models of Sec. 3, while in Appendix B we give details of our numerical treatment of the models of Sec. 4.

Generalities
In this paper we consider a subclass of the Horndeski theories. The general form of the Lagrangian for this subclass is coming from lower order ones, the naive estimate for the strong coupling scale may break down [49]. The dominance of the cubic order shows that this is not the case in models we consider.
It is convenient for our purposes to work in the Arnowitt-Deser-Misner (ADM) formalism. The ADM form of the metric is where γ ij is the three-dimensional metric, N is the lapse function and N i = γ ij N j is the shift function. In ADM terms, the action for the Horndeski theory subclass (3) has the form where The relationship between the two formalisms is established by choosing the equal-time slices as slices of constant φ and defining the time coordinate in such a way that φ(t) is a prescribed monotonous function,φ > 0 (as an example, it is convenient to choose at large positive times e φ = t). This gives .
Then one has [40,52,53] where It is worth noting that general relativity (GR) description of gravity is restored for B 4 = −A 4 = M 2 P /2 = const, where M P is the reduced Planck mass, which we set equal to 1 in what follows. Note also that the transition from the "covariant" formulation (3) to ADM action (4) is not unique: it depends on the choice of the function φ(t). Thus, one can impose additional constraints on the functions A 2 , A 3 , B 4 , or, in other words, on the background solution. We will use this freedom in Sec. 3.2.1.
The equations of motion for homogeneous, isotropic, and spatially flat background are obtained by setting N = N (t), γ ij = a 2 (t)δ ij and varying the action (4) with the respect to N (t) and a(t). They read [54] where H =ȧ/(aN ) is the Hubble parameter. To perform the stability analysis, one writes where a(t) and N 0 (t) are background solutions, ∂ i N T i = 0 and Note that the ADM formulation automatically implies the unitary gauge, δφ = 0. The residual gauge freedom is fixed by setting Y = 0 and W T i = 0, so that the spatial part of the metric reads γ ij = a 2 e 2ζ (e h ) ij .
Variables α, β, and N T i enter the action without temporal derivatives; the dynamical degrees of freedom are ζ and transverse traceless h ij , i.e., scalar and tensor perturbations.
In what follows we omit subscript 0 in the notation for the background lapse function. Then the quadratic action for tensor perturbations reads [12] where Likewise, the quadratic action for scalar perturbation ζ is [12] S ss = dtd 3 xN a 3 G Sζ where To avoid ghost and gradient instabilities, one requires that We also require that the speed of perturbations does not exceed the speed of light, It has been argued that the latter conditions are necessary for the existence of the UV completion [55,56].
It is worth noting that under rescaling of the Lagrangian functions with constant λ, one has A 4 (t, N ) → A 4 (λt, N ), and solutions to equations of motion (10) scale as while the coefficients of the quadratic action transform as In particular, stability and subluminality conditions remain intact under the transformation (16). The scaling property of (16), (17), and (18) implies, in particular, that the overall time scale of evolution can be chosen at one's will, so that at epochs described by GR it is safely longer than the Planck time.

Ansatz
Our purpose in this paper is to design the functions A 2 , A 3 , A 4 , and B 4 in such a way that the model admits a cosmological solution of interest. To this end, we do not need to work in complete generality. To construct bouncing cosmologies in this Section, we make use of the following Ansatz : where µ > 0 is a time-independent parameter which may be different for different cosmologies, 4 f = f (t) is a positive function of time which is extracted as a prefactor in (19a) and (19b) for convenience. Functions a 2 , a 3 are given by Thus, our Ansatz generalizes Ref. [24] and involves four arbitrary functions of time f (t), x(t), v(t), and y(t). The construction of concrete cosmological models boils down to the design of these functions.
Given the Ansatz (19), the background equations of motion (10) are while the functions entering (13) are given by Thus We always choose the functions f (t), . . . y(t) in such a way that inequalities (14) and (15) are satisfied.

Bounce followed by inflation
In this Section, we construct a linearly stable bounce solution which evolves through the following stages: • contraction, with power-law behavior of the scale factor • bounce • inflation, with the Hubble parameter almost constant in time • kination, with the Horndeski field reduced to a massless scalar field.
Unlike at contraction and bounce, gravity at inflation and kination is described by conventional GR, and the expansion is driven by the scalar field.
To build the model, we make use of the following approach. We design the functions f (t), x(t), v(t), and y(t) in (19) and (20) for the contraction, inflation, and kination epochs separately and then construct smooth interpolations between these epochs. One of these interpolating stages involves bounce, and we have to figure out the conditions for its realization. Needless to say, we have to ensure stability and absence of superluminality throughout the whole evolution. Clearly, the construction involves a lot of guesswork, some of which is sketched in what follows. The existence of a consistent solution is ultimately proven by a numerical example.

Early times: the Universe contracts
We begin with the earliest epoch, i.e., large negative times. We require the power-law contraction with constant lapse function, Here we set N = 1 by making use of the ambiguity of transition from covariant to ADM formalism pointed out in Sec. 2 after Eq. (9). As we discussed there, this is equivalent to imposing a constraint on the functions A 2 , . . . , A 4 , or, in other words, on f (t), . . . , y(t). We will encounter this constraint in due course, see Eq. (26).
The desired behavior (23) is achieved by choosing where c, x 0 , v 0 , and y 0 are constant parameters. Our next purpose is to find the complete set of constraints on these parameters. There are several sources of these constraints.
(i) We have to ensure that the background equations (21) are satisfied (with N = 1). Making use of (23) and (24) one finds that the background equations reduce to algebraic equations This set of equations determines the Hubble coefficient χ and also constrains the values of c, x 0 , v 0 , and y 0 . The latter constraint is precisely the one that ensures N = 1. For an appropriate root 5 of (25), the constraint can be written as follows: Then the Hubble coefficient is given by So, the first set of constraints on the parameters defining the model at early times is that v 0 is not arbitrary but is given by Eq. (26), it must be real (argument of square root must be positive), and the Hubble parameter given by (27) must be positive, (ii) The second set of constraints comes from the stability requirement (14) and the absence of superluminal propagation (15). We make use of (22) and write Thus, the constraints F T , G T > 0 are satisfied automatically and c T ≡ 1, while the constraints F S , G S > 0, c 2 S ≤ 1 are nontrivial (but time independent).
(iii) One more constraint comes from the desire to get around the no-go theorem of Refs. [23,24]. A necessary condition for having a stable bouncing solution in the Horndeski theory with GR asymptotics as t → +∞, is [24] i.e., this integral must be convergent in the lower limit of integration. At large negative times we have so the convergence of the integral requires (iv) Yet another constraint is obtained by requiring that despite the fact that F T , F S , G T , and G S (effective Planck masses squared) tend to zero as t → −∞, the background evolution can be described classically at early times. We consider this issue in Appendix A along the lines of Ref. [47]. The outcome is simple: the classical treatment of early time evolution is legitimate provided that Note that this constraint together with (29) implies that i.e., the contraction velocity |ȧ| increases.
In our framework this phenomenon manifests itself in the behavior of superhorizon tensor (and also scalar) modes in the contracting Universe 6 : in the BKL case, one of the two solutions for a superhorizon mode of given conformal momentum grows as t increases and diverges in the formal limit t → 0 (while another solution stays constant in time).
This means that the Universe becomes strongly anisotropic and inhomogeneous at late times, which is undesirable (see, e.g., Ref. [60] for discussion). To avoid BKL, one makes sure that the time-dependent superhorizon solution decays, instead of growing, as t increases towards zero. In our framework, the equation of motion for superhorizon perturbation is obtained from (11) with spatial derivatives neglected, One of its solutions is constant in time, while another is It decays as t increases towards zero for This constraint ensures also that the BKL phenomenon is absent for scalar perturbations. Given that µ < 1, it is weaker than (29).
Thus, the parameter µ in the Lagrangian and the Hubble coefficient χ must belong to the intervals [constraints (i), (iii), and (iv)] To find the allowed range of parameters entering the Lagrangian, we note that in accordance with Eqs. (16), (17), and (18), both equations of motion for background and constraints coming from the absence of instability and super- (26), (27) and (28). The parameter v 0 is given by Eq. (26) in terms of other parameters. So, it is sufficient to determine the allowed range of x 0 , y 0 for a given value of µ and one value of c. With the reduced Planck mass set equal to 1, c is roughly the inverse characteristic time scale in Planck units, so it should be small. 7 In our numerical example below we set c = 4 · 10 −3 , and here we stick to this choice. The allowed ranges of x 0 and y 0 for several values of µ are shown in Fig. 1. In fact, the allowed range is not empty in the entire interval 1/2 < µ < 1; this is illustrated in Fig. 1.
It is useful to note that the asymptotics of v(t) and x(t) may be chosen in such a way that This is possible for all allowed values of µ. In what follows we consider this case only.

Inflation after bounce
As outlined in the beginning of this Section, our next step is to describe the inflationary stage, and then discuss the transition from contraction to inflation through bounce. Models of inflation in the Horndeski theory have been proposed in Refs. [3,54,10,12,61,62]; here we employ the construction similar to Ref. [24]. Namely, exact exponential The white black-framed area shows the allowed range of parameters x 0 and y 0 , where all constraints of this Section are satisfied. We set µ = 0.6 in the left panel (Fig. 1a), µ = 0.8 in the central panel (Fig. 1b) and µ = 0.95 in the right panel (Fig. 1c); c = 4 · 10 −3 everywhere.
expansion occurs when the functions in the Lagrangian take constant values, where the choice (33a) is made to restore GR already at inflation (recall that G 4 = −A 4 = f /2). With the Ansatz and we denote the (time-independent) solution to these equations by H = H 1 and N = N 1 . We require H 1 > 0, parameters using (34): Let us turn to the requirements of background stability and subluminal propagation of perturbations. Using (22) and (35), we arrive at the constraints These constraints can be written in a simple form: which also leads to It is worth noting that the value of y at inflation has the opposite sign to its value at contraction: and y 0 < 0 (see Fig. 1), respectively. On the contrary, the values of x and v have the same signs at these two stages, see Eqs. (32) and (37).
Now, let us comment on the possible behavior of the functions x(t), v(t), y(t) and f (t) near the bounce. Since x(t) and v(t) do not change signs during the transition from contraction to inflation, it is natural to take them monotonously changing from x 0 to x 1 and from v 0 to v 1 , respectively. The function f (t) flattens out, withḟ < 0 both at contraction and bounce. Near the bounce, at t ≈ t b , we have H(t) ≈ 0. Assuming that the functions x(t) < 0 and v(t) > 0 vary slowly in comparison with f (t) in the vicinity of the bounce, we obtain from the first equation of motion, Eq. (21a) with H(t b ) = 0, that We note in passing that with our choice of signs of x and v at contraction, inflation and bounce ( the argument of square root is positive. Then, we take the time derivative of Eq. (21a) and solve it together with Eq. (21b) forḢ andṄ to finḋ , where we again neglectẋ(t),v(t), andẏ(t) in comparison withḟ (t). SinceḢ(t b ) > 0, we see that y(t b ) may be negative, but only slightly: .
With this qualification, most of the smooth functions interpolating between (24) and (33) indeed give rise to the bouncing solution. We present a numerical example in Sec. 3.2.4.

Kination epoch after transit from inflation
To describe the final kination epoch with GR and free massless scalar field, we make use of the covariant formalism with the Lagrangian (3). It is convenient to use the freedom of field redefinition and choose the background field φ as follows: This choice corresponds to the Lagrangian Indeed, it is straightforward to check that the scalar field equation and Friedmann equation have the solution (39) with a = const · t 1/3 , N = const, and H = (3tN ) −1 . Note that during the transition from inflation to kination, the coefficient of X in the Lagrangian changes sign (it changes from x 1 < 0 at inflation to 2/3 at kination). This is in complete accordance with [50,51].
Towards the kination epoch, other terms in the scalar field Lagrangian should tend sufficiently rapidly to zero.
This can be achieved, e.g., by requiring the following asymptotics of functions in the covariant Lagrangian (3) at large φ (distant future; recall that GR is restored already at inflation, G 4 = 1/2): where ω 2 (φ) and ω 3 (φ) are damping factors, which suppress higher order terms. They can be chosen rather arbitrarily.
We choose them as follows: The reason for this choice is that we obtain simple ADM functions A 2 and A 3 using conversion formulas (6)- (9): In fact, we have to generalize these expressions by introducing time shift, t → t − t * , where t * is a new parameter. The point is that once contraction is described literally by formulas of Sec. 3.2.1, inflation begins soon after t = 0 and ends at some much later time t e . The time shift (with t * of order of t e but somewhat smaller than t e ) has to be introduced to account for the fact that kination begins around t e = 0. Thus, the asymptotic behavior of x(t), v(t), and y(t) at large t is (recall that f (t) = 1 at inflation and later) By choosing the functions x(t), v(t), and y(t) in such a way that they interpolate between constant values x 1 , v 1 , and y 1 at inflation and functions (41) at late times, we obtain a smooth transition from inflation to kination. The issues of stability and subluminality are, however, tricky at transition epochs; designing a completely stable and subluminal model requires considerable trial and error effort.

Numerical example
Here we present a concrete model which proves by example that there exists stable and subluminal cosmology with the desired properties listed in the beginning of this Section. We emphasize again that by rescaling the functions in the Lagrangian in accordance with Eq. (16), one can make all time scales like c −1 , inverse inflationary Hubble parameter H −1 1 , etc., arbitrarily long, much longer than the Planck time. This observation applies also to other models considered in this paper.
We choose the parameter µ near the center of allowed interval 1/2 < µ < 1: As we already mentioned, at the contracting stage we choose, quite arbitrarily, The parameters x 0 and y 0 are then chosen from the allowed region shown in Fig. 1b; by trial and error we find convenient values, consistent with the sign choice (32): The value of v 0 and the Hubble coefficient χ at the contraction stage are found from (26) and (27): v 0 = 5.19 · 10 −6 , χ = 0.01 .
This completes the description of the contraction stage.
We would like to have bounce at some time before t = 0; we request [although we do not have to do so in view of scaling (16)] that the characteristic time scales are large compared to 1 (i.e., Planck time). We begin with the function f (t) which should interpolate between f = −ct at contraction and f = 1 at inflation. A simple choice is The parameter s is the inverse time scale of the transition, and we set We wish the bounce to occur roughly at t ∼ −s −1 , so the maximum value of |H| at contraction is estimated as Let us now turn to the inflationary stage and transition to it. A simple Ansatz for the inflationary Hubble parameter H 1 is that it is comparable to the maximum value |H| max at contraction. There is no reason to think that the lapse function at inflation is considerably different from 1. We choose H 1 = 3.7 · 10 −5 , N 1 = 0.82 .
We also have to specify the value of y 1 at inflation. To this end, we introduce a simple Ansatz of proportionality between x(t) and v(t): so that v(t)/x(t) is time independent at the transition from contraction to inflation. Note that with the numerical values (43) and (44), the estimate (38) is consistent with our choice N ∼ 1. Equation (47) implies then (35) gives and numerically This set of parameters is consistent with the stability and subluminality constraints (36). The values of x 1 and v 1 are found from (35): Note that |x 1 |, v 1 , and y 1 are much smaller than |x 0 |, v 0 , and |y 0 |, respectively. This has to do with two properties of the contracting stage which distinguish it from inflation. First, the constraints shown in Fig. 1b are consistent with fairly large values of |y 0 |, and it is indeed quite large in our example; on the contrary, inflationary y 1 is bounded by H 1 , see (36b). Second, equation of motion (21b) contains a term proportional toḟ y 0 = −cy 0 , which is not so small at contraction and drives |x 0 | and v 0 to fairly large values; this term vanishes at inflation. We note in passing that Eq. (21a) is satisfied at the contraction stage due to the partial cancellation between x 0 and 3v 0 .
The transition from contraction through bounce to inflation is described by x(t), v(t), and y(t) smoothly interpolating between x 0 , v 0 , y 0 and x 1 , v 1 , y 1 [and with f (t) given by (45)]. A nontrivial part of the construction is to make sure that stability and subluminality conditions (14) and (15) are satisfied. By trial and error, we find appropriate forms where the functions U x (t) = ln e −1.5·s·(t−80) + e 2 e −1.5·s·(t−80) + e (49a) U y (t) = ln e −3.8·s·(t+180) + e 2 e −3.8·s·(t+180) + e (49b) interpolate between 0 and 1, while v(t) is given by (47), and the parameter s is the same as in (46).
We show the behavior of the Hubble parameter and lapse function at contraction, bounce, and beginning of inflation in Fig. 3. The scalar coefficients F S and G S and scalar sound speed c S are shown in Figs. 4a, 4b and 5a, respectively; the stability and subluminality are explicit. We show tensor coefficient F T for completeness in Fig. 5b (recall that G T = F T and c T = 1 at all times). So, after a rather short transition period, the inflationary stage sets in. Depending on the parameters of the model, it can last for a longer or shorter time. Note that this property may be of interest from a phenomenological viewpoint [62]. We take, quite arbitrarily, the duration of inflation approximately equal to ∆t inf ≈ 1.55 · 10 6 (in Planck units), which corresponds to the number of e-foldings at inflation N e = N 1 H 1 ∆t inf ≈ 46. To have the transition from inflation to kination, we take at late times with the parameter which regulates the duration of inflation equal to t * = 1.5 · 10 6 , and where T = 4.5 · 10 4 , and the function again interpolates between 0 and 1 [the value of parameter s is still given by (46)]. The parameters t * and T are chosen in such a way that the functions x(t), v(t), and y(t) are reasonably smooth in the transition period, and inflation smoothly ends somewhat later than t * . This is illustrated in Fig. 6. At late times we obtain the correct kination behavior given by (41), and the Hubble parameter asymptotes to H = [3(t − t * )N ] −1 . The Hubble parameter and lapse function are shown in Fig. 7, while the scalar coefficient F S and scalar sound speed c S are shown in Fig. 8.
Clearly, the model is stable and subluminal at the transition from inflation to kination. The sound speed tends to 1 rather slowly, since the ratio v(t)/x(t) exhibits slow decay (t − t * ) −3 (and y(t) decays as (t − t * ) −5 ).
To end up this Section, we note that since the duration of inflation is fairly long, the complete expressions for x(t), v(t), and y(t), valid at all times, are obtained by simple superpositions of (48) and (50), e.g., etc. This completes our discussion of the model with bounce, inflation, and kination.

Bounce directly to kination
Contraction and bounce need not necessarily proceed into the inflationary stage: a short transition epoch after bounce may end up directly at kination. In this scenario, the initial stage is described in the same way as in Sec. 3.2.1, whereas the evolution after bounce proceeds as in Sec. 3.2.3. Let us give a numerical example which shows that stable and subluminal cosmology of this sort is indeed possible. We again consider a model with µ = 0.8 and choose the parameters of the contraction stage as in (42), (43), and (44). The function f (t) is again given by (45), so that we restore GR at later times. We would like to approach the behavior (41) soon after bounce and, by trial and error, end up with the following example: where U x (t) and U y (t) are still given by (49), and now v 2 = 1.04 · 10 8 , y 2 = 9.6 · 10 10 .
We show the Hubble parameter and lapse function for this model in Fig. 9 and the scalar coefficient F S and scalar sound speed in Fig. 10. The latter figure illustrates that the model is stable and subluminal at all times.

Ansatz
To illustrate that interesting cosmologies can be obtained within various Ansätze, in this Section we construct genesis models by choosing the functions in the Lagrangian (5) in the following form: The parameter µ > 0 is similar to that in the previous Section, and the parameter δ > 0 is new. The functions A 4 and B 4 depend on N now; functions a 2 (N ), a 3 (t, N ), a 4 (t, N ), and b 4 (t, N ) are chosen as follows: Note that the parameter x is now time independent, unlike in the models of Sec. 3 where functions x(t) and v(t) entering a 2 exhibited step-function behavior.
Given the Ansatz (52), the background equations of motion (10) are The functions entering (13) are given by We will make sure that the functions f (t), y(t) and z(t) are such that inequalities (14) and (15) are satisfied.

Contracting genesis followed by bounce
Here we construct a model with genesis of contracting Universe. This cosmology begins with the flat space-time, then the Universe starts to contract, and the rate of contraction increases. At some moment of time the bounce occurs: the contraction terminates and expansion begins. We consider for definiteness the case in which bounce is followed by inflationary expansion; inflation is assumed to end up as in Sec. 3.2. Alternatively, bounce may lead directly to kination epoch like in Sec. 3.3; we do not elaborate on this possibility.
We begin with early times and consider the following asymptotics motion (54) we arrive at [the fact that N = 1 for this solution is due to the particular choice of the coefficient (−1/3) of N −4 in (53a); this choice replaces in this model the constraint (26)]. For δ > 0, the scale factor tends to a constant as t → −∞: as required for genesis. Contracting genesis occurs for χ > 0.
Let us note that for small δ, early time asymptotics is approached slowly (backwards in time), since the expansion parameter is (−t) −δ : where coefficients χ 1 and N 1 are not particularly small. This complicates the numerical analysis; we explain how we get around this obstacle in Appendix B. The same comment applies to the genesis model of Sec. 4.3.
The asymptotic behavior of coefficients (55) entering the quadratic action for perturbations is The constraints on parameters arise from the same requirements as in Sec. 3.2.1. Let us list them: (i) Contraction at early times: (ii) Stability of background and subluminality of perturbations: (iii) Evading the no-go argument of Ref. [24]: (iv) Absence of strong coupling in the past. This issue has been studied in detail in Ref. [47] with the result µ + 3 2 δ < 1; (v) Belinsky-Khalatnikov-Lifshitz phenomenon. In the same way as in Sec. 3, we obtain for time-dependent superhorizon perturbations They decay as t increases towards zero, provided that 2µ + 1 > δ. This constraint is weaker than (iii).
All these constraints can be satisfied without much of fine-tuning. In view of (61), the constraints give The constraint F S > 0 is satisfied automatically. This set of inequalities can be satisfied for all µ and δ from their allowed range.
Let us turn to inflationary epoch. Like in Sec. 3.2, inflation occurs for time-independent coefficients in the Lagrangian: In this case, equations of motion (54) read and we denote the (time-independent) solution to these equations by H = H 1 and N = N 1 . We require H 1 > 0, In analogy to Sec. 3.2, it is convenient to treat H 1 and N 1 as independent parameters and express y 1 and z 1 through these parameters using (64): We emphasize that unlike in the bounce scenario of Sec. 3, the inflationary epoch in our current model is not described by GR since z 1 = 0. The conditions for the background stability and subluminal propagation of perturbations, Eqs. (55), read These can also be satisfied without much of fine-tuning. As an example, the allowed range of x and H 1 is shown in To see that the contracting genesis stage can consistently pass through bounce to inflationary expansion, we now give an explicit numerical example. As we alluded to above, we choose µ = 0.8, δ = 0.1. The parameter relevant to the contracting genesis stage is chosen as c = 1.7545 · 10 −2 . By trial and error we find convenient values for other Lagrangian parameters at early times, consistent with the system of constraints (59)- (62) and (66): The value of the Hubble coefficient χ at the contraction stage is found from (57), χ = 0.25. Next, we turn to the function f (t) which should interpolate between f = −ct at contraction and f = 1 at inflation. To ensure stability and subluminality at all times, we choose this function in somewhat more complicated form than in (45): The transition from contraction to inflation is described by y(t) and z(t) smoothly interpolating between y 0 , z 0 and y 1 , z 1 . By trial and error, we find appropriate functions where the functions U y (t) = 1 + ln e 3.8·s·(t+150) + e e 3.8·s·(t+150) + e 2 , U z (t) = ln e −5.8·s·(t−605) + e 2 e −5.8·s·(t−605) + e interpolate between 0 and 1. A nontrivial requirement leading to (69) is again stability and subluminality of perturbations at all times.
We show the behavior of the Hubble parameter and lapse function at contraction, bounce and beginning of inflation in Fig. 12. The scalar coefficient F S and scalar sound speed c S are shown in Fig. 13: the stability and subluminality are explicit (although not obvious in Fig. 13a, the coefficient F S is, in fact, strictly positive at all times). We end up this Section by the following remark. As we pointed out above, gravity at the inflationary epoch is [see (68)] and some positive value y 2 . The latter property follows from the expression for F S , which for z = 0 has the following form: where N 2 > 0 and H 2 > 0 are lapse function and Hubble parameter at GR inflation. It is straightforward to design appropriate functions z(t) and y(t) without spoiling the stability and subluminality properties. An example is where y 2 = 2 · 10 −4 , t * * = 1.8 · 10 4 and U z1 (t) = ln e −5.8·s·(t−500) + e 2 e −5.8·s·(t−500) + e , U y1 (t) = 1 + ln e 3.8·s·t + e e 3.8·s·t + e 2 .
The transition from late GR inflation to kination proceeds in the same way as in Sec. 3.2.

Genesis without strong coupling
As we pointed out in the Introduction, the genesis model of Ref. [24] suffers from the strong coupling problem at early times. For completeness, we present here a version of this model which is free of the strong coupling problem. The Lagrangian is the same as in [24] but with different parameters. Namely, the Lagrangian functions are defined by Ansatz (52) with x = const, y = const, and z = 0, i.e., where µ > 0, δ > 0 and As before, we choose the asymptotic behavior f = −ct (c > 0) as t → −∞, and using equation of motion (54) obtain the genesis behavior where ξ is given by The asymptotics of the scalar coefficients and scalar sound speed squared are x + 3c · y · (2µ + δ + 1) 2 , Parameters of this model should obey several constraints. The first one comes from the requirement of evading the no-go argument of Ref. [24]: The second constraint ensures that the classical treatment of early time evolution is legitimate [47]: The third one is the requirement that the Universe expands at early times: Finally, one requires the background stability and subluminal propagation of perturbations. For y > 0 (as needed for healthy inflation, see below), all these constraints are satisfied, provided that x < 0 and x < −3c · y(2µ + δ + 1) .
The transition from the genesis stage to inflation is achieved simply by flattening out the function f (t) to f = 1. For f = 1, equations of motion (54) read and we denote the (time-independent) solution to these equations by H = H 1 and N = N 1 . We require H 1 > 0, The requirements of background stability and subluminal propagation of perturbations, Eqs. (55), read in this case As before, it is convenient to treat H 1 and N 1 as independent parameters and express x and y through these parameters using (76): .

(78b)
Then the constraints (75) and (77) reduce to In accordance with (78), these constraints translate into constraints on x and y. It is straightforward to see that the latter are satisfied, provided that y > 0 and x obeys (75).
Now, let us turn to our numerical example. We choose We obtain the values of x and y by considering the inflationary stage. We choose H 1 and N 1 at inflation as follows:

Conclusion
This paper demonstrates that it is relatively straightforward to construct, within the Horndeski class of scalar-tensor theories, nonsingular cosmological models which are healthy, i.e., free of instabilities and superluminal propagation of perturbations. The price to pay is strong gravity in the past, the property that effective Planck masses tend to zero as t → −∞. We have made sure, however, that the latter property does not spoil the description of the background within classical field theory. In this way we have constructed bouncing models, genesis cosmology, and a combination thereof. These may or may not pass through the inflationary stage, as we explicitly demonstrated in Sec. 3.
In our constructions, we heavily exploited the functional freedom that exists in the Horndeski class of theories.
On the one hand, this freedom is instrumental for designing models with prescribed properties (in other words, for employing the "inverse method" [21]); on the other, it makes the whole approach not so appealing. It is certainly desirable to have a better idea of which Horndeski theories, if any, have a chance to be realistic as low energy effective theories.
written in terms of metric variables, with the Lagrangian (5) (recall that N denotes the background lapse function). Different terms in the Lagrangian (5) contain different powers of the scale factor, which in the contracting Universe nontrivially depends on time, It is therefore convenient to work with physical momenta and frequencies. Assuming that they are higher than E class , i.e., assuming that E strong E class , we can neglect slow dependence of the scale factor on t and at a given time treat a as an (instantaneously) time-independent parameter (with the exception of expressions that involve the Hubble parameter explicitly). Of course, this assumption must be justified a posteriori: the whole analysis is valid if the classical treatment of the background is legitimate, E strong E class . Having this in mind, we introduce physical spatial and temporal coordinatesx ≡ xa, t ≡ tN = t.
Note that N = 1 in our case, but we keep the notationt for concordance with spatial coordinates. Then the derivatives are, respectively,∂ We now rewrite our Lagrangian (5) in terms of physical coordinates. We define and find Finally, we have √ −gd 4 x = (1 + α)e 3ζ dtd 3x . Thus, the action written in terms of physical variables does not contain the scale factor anymore, but otherwise the Lagrangian has the same structure as the original one (5), except for the first term in the right-hand side of (A.4).
The action for perturbations, written in terms of physical variablesx i ,t, is identically the same as the limiting case of the action that we encountered in Ref. [47], where we studied the strong coupling issue in the model with genesis described in Sec. 4.3. Namely, the model with genesis has an additional parameter δ > 0 in (70), while the bouncing model we discuss does not. By direct inspection we find that the action for perturbations in the bouncing model, written in terms of physical variables, is obtained from the action for perturbations in the genesis model by sending δ → 0 [this includes also the first term in the right-hand side of (A.4)]. Thus, the sufficient condition for the absence of the strong coupling problem is obtained from the result of Ref. [47] in the limit δ → 0 and reads µ < 1 .
This result is quoted in Sec. 3.
Let us illustrate this constraint by considering the quadratic and cubic action for tensor modes. The complete expression is [63] S (2) hh + S We recall that F T ∝ (−t) −2µ , see (28a). To figure out the associated strong coupling scale, we introduce the canonically normalized field (omitting indices) and find that the interaction term is (modulo numerical factor) Thus, on dimensional grounds, the strong coupling scale is This scale is much higher than the classical scale E class =t −1 for µ < 1, as promised.

B CONTRACTING GENESIS: SUBTLETY OF NUMERICAL SOLU-TION
As we pointed out in Sec. 4.2, corrections to the leading asymptotics of classical solutions at early times behave as (−t) −δ . This makes straightforward numerical treatment problematic for small δ. Here we sketch our way of dealing with this problem. We consider for definiteness the model of Sec. 4.2 with δ = 0.1.
We study the time interval −∞ < t < t 1 , where t 1 is negative and |t 1 | is large enough, so that it is a very good approximation to use the asymptotics Then coefficient u 2/δ factors out in equations of motion (54), and coefficients in these equations become linear polynomials in u. This form of equations of motion enables one to solve them in a straightforward way.
Thus, corrections are small for small enough u = u 0 , where initial conditions are imposed. In practice we choose u 0 = 10 −7 , which corresponds to t 0 = −6 · 10 71 . This huge number is the reason why we have invented our procedure.
We solve the equations of motion written in terms of u until u becomes roughly of order 1; in practice we choose u 1 = (−ct 1 ) −δ = 0.67, so that t 1 = 3000. At that time Eq. (B.1) is still a good approximation. Then we continue solving equations of motion using time t, with obvious matching at t = t 1 .