{'fraction_non_alphanumeric': 0.03288699234844026, 'fraction_numerical': 0.012580929958799295, 'mean_word_length': 4.474892588614393, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stable bounce and inflation in non-local higher derivative cosmology

One of the greatest problems of primordial inflation is that the inflationary space-time is past-incomplete. This is mainly because Einstein's GR suffers from a space-like Big Bang singularity. It has recently been shown that ghost-free, non-local higher-derivative ultra-violet modifications of Einstein's gravity may be able to resolve the cosmological Big Bang singularity via a non-singular bounce. Within the framework of such non-local cosmological models, we are going to study both sub- and super-Hubble perturbations around an inflationary trajectory which is preceded by the Big Bounce in the past, and demonstrate that the inflationary trajectory has an ultra-violet completion and that perturbations do not suffer from any pathologies.


Introduction
Primordial inflation is one of the most compelling theoretical explanations for creating the seed perturbations for the large scale structure of the universe and the temperature anisotropy in cosmic microwave background radiation [1]. Since inflation dilutes everything, the end of inflation is also responsible for creating the relevant matter required for the Big Bang Nucleosynthesis besides matching the perturbations observed in the CMB, see for instance [2,3], where inflation has been embedded completely within a supersymmetric Standard Model gauge theory which satisfies all the observed criteria, for a recent review on various models of inflation, see Ref. [4].
In spite of the great successes of inflation, it still requires an ultra-violet (UV) completion. As it stands, within Einstein's gravity inflation cannot alleviate the Big Bang singularity problem, rather it pushes the singularity backwards in time [5,6]. In a finite time in the past the particle trajectories in an inflationary space-time abruptly ends, and therefore inflationary trajectory cannot be made past eternal [7]. Within the framework of Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology, this can be addressed by ensuring the inflationary phase to be preceded either by an "emergent" phase where the space-time asymptotes to a Minkowski space-time in the infinite past [8], or by a non-singular bounce where a prior contraction phase gives way to expansion. Unfortunately, within the context of General Relativity (GR), as long as the average expansion rate in the past is greater than zero, i.e. H av > 0, the standard singularity theorems due to Hawking and Penrose hold, and the fluctuations grow as the universe approaches the singularity in the past, which inevitably leads to a collapse of the space-time [9] commonly referred to as the Big Bang (see also [7]). Thus, in order to provide a geodesic completion to inflation one needs to modify GR. 1 There have been some recent progress in ameliorating the UV properties of gravity in the context of a class of "mildly" non-local 2 higher derivative theory of gravity [12][13][14], involving terms that are quadratic in the Reimann curvature but contains all orders in higher derivatives. It was found that such an extension of Einstein's gravity renders gravity ghost-free and asymptotically-free in the 4 dimensional Minkowski space-time. 3 As a result the ordinary massless graviton sees a modified interaction and propagation at higher energies, but at low energies in the infra-red (IR) it exactly yields the properties of Einstein's gravity. It has been demonstrated that such a modification of Einstein's gravity ameliorates the UV behaviour of the Newtonian potential, and also permits symmetric non-singular oscillating solutions around the Minkowski space-time [14] signalling that the action can resolve the cosmological singularities, see also [12,13,19]. Similar non-local higher derivative theories of gravity have also been considered in connection with black holes [20], inflationary cosmology [21], string-gas cosmology [22] and in efforts to understand the quantum nature of gravity [23]. 4 There was also a proposal to solve the cosmological constant problem by a non-local modification of gravity [26] (see also [27]).
What is rather intriguing is that non-local actions with an infinite series of higher derivative terms as a series expansion in α , the string tension, have also appeared in stringy literature. One mostly finds them in the open string sector, notably in toy models such as p-adic strings [28], zeta-strings [29], strings on random lattice [30,31], and in the open string field theory action for the tachyon [32]. For more general aspects of string field theory see for instance [33]. There have also been progress in understanding quantum properties of such theories, i.e. how to obtain loop expansion in the string coupling constant g s , which has lead to some surprising insights into stringy phenomena such as thermal duality [34], phase transitions [35] and Regge behavior [31]. It is indeed a wishful thinking to derive such an action in the closed string sector involving gravity. 5 In Ref. [12], with the means of an exact solution, it was shown how a subset of the actions that have been proposed more recently in Ref. [14], can ameliorate the Big Bang cosmological singularity of the FLRW metric. The action discussed in [12] only contained the Ricci scalar and it's derivatives up to arbitrary orders and yielded non-singular bouncing solution of hyperbolic cosine type. However, it was not clear how generic these solutions were, whether these solutions were stable under small perturbations, and whether these theories contain other singular solutions or not.
Some of these issues, in particular concerning only time dependent fluctuations around a bouncing solution were addressed in Ref. [13], where the authors have analysed the stability of the background trajectory for very long wavelength (super-Hubble) perturbations, either when the cosmological constant is positive, i.e. Λ > 0, or negative, i.e. Λ < 0. 6 The latter provided a way to construct cyclic cosmologies, and in particular cyclic inflationary scenarios that have been studied in [40], and in [41] where cosmological imprints were also discussed.
In this paper we aim to focus on the case when Λ > 0. This particular case is very interesting from the point of view of inflation, as it paves the way for a geodesically complete paradigm of inflation in past and future. It was found that there exists a de Sitter solution as an attractor for large enough times before and after the bounce. The weakening of gravity in the UV regime facilitates the bounce for the FLRW metric. Our main goal will be to study what happens to the non-singular background solution under arbitrary classical and quantum fluctuations (super and sub-Hubble type) during the bounce period. We have a two-fold aim. First, we want to verify that at late times there are no fluctuations that are growing during inflation. This would clearly signal an instability in the system which would certainly jeopardize the inflationary mechanism of creating superHubble fluctuations which becomes a constant at late times. Second, we want to check that the bounce mechanism itself is robust, i.e. the perturbations do not grow unboundedly near the bounce. In the process of our investigations, we will develop techniques to deal with non-local perturbation equations that one encounters in these theories, and this we believe is a major achievement of the paper. We hope that this will help us in analyzing other cosmological issues in future. Finally, let us comment on a rather interesting feature we noted. It turns out that the requirement of the bounce automatically introduces an extra scalar degree of freedom, and if the mass of this scalar is light, it may be able to play the role of a curvaton [4]. We however leave a more detailed phenomenological investigation of such a scenario for future.
The paper is organized as follows: In section 2, we review some of the results for non-local gravity models and generalize some of the techniques previously used to find exact solutions. In section 3, we derive the non-local perturbation equations for the Bardeen potentials around a large class of cosmological background solutions. Next, in section 4 we proceed to solve these equations around the "inflationary" de Sitter late-time attractor solutions that exist in these models. In particular, we obtain the criteria for which the perturbations freeze out or decay after crossing the Hubble radius which allows one to maintain the inflationary mechanism of producing scale-invariant CMB fluctuations. Thereafter, we look at the stability of the fluctuations in section 5 around the bounce. Finally in section 6, we conclude by summarizing our results and discussing their implications for inflationary cosmology.

String-inspired Non-local Gravity
In this section we begin by introducing the String theory motivated non-local gravity models. We will then generalize the methods previously employed in [12,13] to obtain exact solutions in these theories. We will explore the various non-singular cosmologies that are present in these theories and in particular focus on solutions which can geodesically complete the inflationary space-time. Although the inflaton potential energy does not remain a strict constant and changes slowly, for technical reasons in our analysis we will treat the potential energy as a "cosmological" constant. This should not however affect the issues we plan to address in this paper.

Action and Equations of Motion
In Refs. [12,13] the cosmology of a special class of ghost free actions was investigated. The action has the form where M P is the Planck mass: M 2 P = 1/(8πG N ), G N is the Newtonian gravitational constant, Λ is the cosmological constant, M * is the mass scale at which the higher derivative terms in the action become important, and λ is essentially a book keeping device to help us keep track of the higher derivative corrections. L M is the matter Lagrangian. We use the convention where the metric g has the signature (−, +, +, +).
The analytic function F( /M 2 * ) = n 0 f n n is an ingredient inspired by String Theory as an expansion in the string tension (α ). We recall that the action of the covariant d'Alembertian on a scalar is given by where ∇ µ is the covariant derivative. Variation of action Eq. (2.1) yields the following equations: where T M µ ν is the energy-momentum tensor of matter and is the Einstein tensor. We assume that the matter stress energy tensor obeys the usual conservation equation, ∇ µ T M µ ν = 0, to ensure consistency of the modified Einstein's equations, which we will refer to in future as Einstein-Schmidt (ES) equations 7 . The ES equations can be written in a reasonably compact form:

4)
7 The generalization of the Einstein's equations to include terms such as R n R was analyzed by Schmidt [42], having an infinite series of such terms give rise to Eq. (2.1).
where we have introduced two additional quantities: Making any progress with the ES equations may seem daunting at first, but previous studies have shown that the trace equation can be tractable since all we have to deal with is the Ricci scalar: (2.6) or, in the more compact notation

General Ansatz for finding Exact Solutions
It has been shown in Refs. [12,13,43] that the following ansatz is useful in finding exact solutions. From Eq. (2.8) one gets the following relationships: If the scalar curvature R satisfies (2.8), then equations (2.2) have the following form: where F (1) is the first derivative of F with respect to the argument, and the trace equation (2.6) becomes especially simple: where We proceed to consider a traceless radiation along with a cosmological constant. This will simplify many of our calculations. Then the simplest way to obtain a solution to Eq. (2.11) with T M µ µ = 0 is to put A 1 = A 2 = 0, and impose This implies Thus we see that within the context of the ansatz, the ES equations and the corresponding solutions are characterized by three independent parameters, r 1 , λf 0 , and λF 1 . Although λ is only a book-keeping device, and one can set it to one if desired, it is interesting to look at the analytic properties of the different parameters with respect to λ since the λ → 0 limit can mimic the large time limit when one expects to recover GR and all the higher derivatives become small. Now, r 1 is determined by (2.12) and is therefore independent of λ, but r 2 depends non-analytically on λ. 8 This could seem to be a potential problem in recovering General Relativistic solutions. As we will see later in explicit examples, although r 2 may have non-analytic dependence, it is really just an intermediate variable of convenience, and the solutions for the scale-factor need not have any singular dependence in λ (see section 2.4). This is why one is able to recover the GR solutions in appropriate limits.
It is important to note that by virtue of the Bianchi identities, solving the trace of the ES equations essentially guarantee that we have a solution for the complete ES equations, which now simplifies to (2.14) In general one is required to include radiative sources to get exact solutions [13], and we need to make sure that the radiative energy density if present is positive and not ghost-like. We will discuss this point in the specific context of cosmological solutions in the following subsection. Let us emphasize that equations up to this point are general and do not take into account the properties of the metric. Thus the above procedure can be used to obtain metrics which are cosmological (time dependent solutions) in nature or which have spatial dependence. This is an important generalization on the analysis discussed in [13], and a major step forward in investigating these non-local models.

Cosmological Solutions, Bounce and Stability
We are now going to focus on spatially flat FLRW type cosmological solutions that was obtained in [12,13] in the presence of a cosmological constant and radiation energy density, with a traceless energy momentum tensor. In the above expression a is the scale factor and a 0 refers to its value at the bounce. In the cosmological context the ansatz Eq. (2.8) becomes a third order differential equation for the Hubble parameter, H =ȧ/a, where dot defines derivative w.r.t physical time. The solutions to the above equation, their stability, and asymptotic behaviours have been studied in details in [13]. It was found that several non-singular solutions exist where the universe bounces from a phase of contraction to a phase of expansion, thus providing a resolution to the Big Bang singularity problem. By looking at the 00 component of the ES equations, one can calculate the radiation energy density at the bounce point, ρ 0 , which must be positive: where h 1 =Ḧ/M 3 * characterizes the acceleration of the universe at the bounce point and plays the role of an 'initial condition'. There are three cosmological scenarios one can envisage from Eq. (2.16) : • Λ < 0, r 1 > 0 ⇒ r 2 > 0: One obtains a cyclic universe scenario where the universe undergoes successive bounces and turnarounds. Qualitatively speaking the bounce is caused by the higher derivative modifications, while the turnarounds occur when the radiation energy density cancels the negative cosmological constant. Such scenarios can give rise to "cyclic inflationary" models if one is able to incorporate entropy production [40,41], but we will not consider these scenarios here any further.
• Λ > 0, r 1 < 0 ⇒ r 2 > 0: One obtains a bouncing universe scenario where a phase of contraction gives way to a phase of super-inflating expansion phase. An exact example is the Ruzmaikin-type solution [44,45], which incidentally contains no radiation [43]. Unfortunately, in these scenarios the universe is stuck in a quantum gravitational phase accelerating faster and faster, and we never recover a GR phase. Clearly such solutions are not phenomenologically viable and hence we disregard these scenarios.
• Λ > 0, r 1 > 0 ⇒ r 2 < 0: This represents a geodesic completion of an inflationary spacetime via a non-singular bounce. It is easy to see that the action Eq. (2.1) without matter admits constant curvature vacuum solutions, the de Sitter solution or the Minkowski solution, depending upon the value of the cosmological constant. These background solutions are in fact exactly the same as that of the general relativistic de Sitter phase, as for constant R, all the higher derivative terms in the field equations vanish, and one just obtains from Eq. (2.6) This is precisely the de Sitter solution the space-time asymptotes to in the more general case, as the radiation dilutes away very quickly during the cosmological constant dominated phase.
It is this last case that is relevant for inflationary cosmology, and is the main focus of our paper. The higher derivative extension of gravity provides us with a plausible answer to the question, "what happened before inflation?". The general picture that emerges is that in the presence of a positive cosmological constant, if a small causal patch is approximately homogeneous and isotropic at Planckian densities, the universe would undergo a non-singular bounce, following which the universe will start to inflate quickly diluting any matter (radiation) present in the system, and we essentially will have a pure de Sitter. Of course, to be consistent with the inflationary phenomenology, Λ cannot be a strict constant but has to evolve, possibly with a scalar field, the inflaton [4]. For technical reasons we haven't been able to incorporate this dynamics yet, but our main aim here is to assess the robustness of the bounce mechanism and the inflationary scenario in terms of perturbations.
What we expect is that with the dilution of radiation energy density during the vacuum energy dominated phase, any fluctuations present in the fluid would also get diluted away, so that the fluctuations can become dominated by the quantum fluctuations of the inflaton field which we know can produce the observed approximately scale-invariant spectrum in CMB. Since we do not have any inflaton in our scenario, what we want to check is whether there are any growing or unstable modes of fluctuations around the geodesically complete de Sitter space-time. If present, this would signal either the fallibility of the bounce mechanism, or the breakdown of the inflationary mechanism of generating fluctuations.

The Hyperbolic Cosine Bounce
It is rather fortunate that one can actually find an exact solution for the higher derivative extension of gravity under consideration, which serves as a good illustration of a geodesic complete space-time. It is given by This solution satisfies the ansatz with the specific parameter combination Substituting r 2 = −6r 2 1 into Eqs. (2.13) and (2.17) (2.20) Note that the radiation is non-ghost like, provided F(r 1 ) < 0. (2.21) This turns out to be a general constraint which must be satisfied in order to realize geodesically complete de Sitter space-time, and we are going to assume that this constraint is satisfied. The above specific example (2.19) serves as a nice illustration of some of our general results and we will fall back on it as appropriate in the remainder of the paper. Finally, one may wonder about the special case when F 1 = 0, which means there is no radiation. Although this is under current investigations, indications are that such theories will suffer from the presence of ghosts. Effectively there is a double pole in the propagator at p 2 = r 1 .

Perturbations
In the previous section we have provided a general algorithm for finding cosmological background solutions in non-local gravity models which in particular can realize inflationary space time at late times. We have also provided concrete examples of such solutions. What we now want to address is the robustness of the paradigm by studying the fluctuations around the background. This will let us test whether the bounce mechanism is stable, as well as whether the inflationary mechanism of generating scale-invariant perturbations can survive in the presence of higher derivative modifications. In this section we will provide a general framework for dealing with the scalar modes and obtain the coupled (non-local) differential equations describing the evolution of the Bardeen potentials. Since this is a somewhat technical section, we have provided a summary subsection 3.4, where we enumerate the results in a self-contained manner.

Perturbations on higher derivative action of gravity
To compute perturbation equations we use the fact that the background configuration obeys the ansatz (2.8) and conditions (2.13). Let us introduce the following notations to separate the background and the perturbations: where the subscript B stands for background. Also, let us denote the variation of the ansatz (2.8) as follows: While computing the variation of the ES equations, what we encounter most often is the variation δ( n R), which thanks to the ansatz (2.8) sums up to and this will turn out to be the key simplification 9 . This can be used to obtain where we have introduced Note that, expanding F( ) in the Taylor series, one can see that δ(F( )R) has no pole. Also thanks to F (1) (r 1 ) = 0, Z( ) is an analytic function. Similarly, we obtain where we have used conditions Eq. (2.13), and in particular F (r 1 ) = 0, to obtain Eq. (3.7).
9 Similar relationships are useful for the perturbation analysis in models with non-local scalar field as shown in Ref. [39].

Perturbing the trace equation
Assuming δT M µ µ = 0, which is true for radiation, we get the trace equation (2.7) for perturbations linearized about a background solution: Substituting Eqs. (3.5) and (3.7), and using condition A 1 = 0 to eliminate M 2 P , we obtain This is a very important equation because it gives us information about ζ. It can be rewritten as Pζ = 0, (3.10) where As we see, operator P depends not only on d'Alembertian but also on partial derivatives and has non-constant coefficients. Nevertheless, we have been able to find a non-local differential equation for one of the perturbative variable ζ, which one can use different techniques to solve. The form of the perturbed trace equation demonstrates the importance of the function ζ.

Scalar Perturbations
The metric perturbations can be divided into 4 scalar, 4 vector and 2 tensor degrees of freedom, according to their transformation properties with respect to the three spatial coordinates on the constant-time hypersurface. Different types of perturbations do not mix at the first order [46]- [50]. In this paper, we are going to focus on the scalar perturbations leaving the detailed analysis of vector and tensor modes for separate projects. Scalar metric perturbations are given by 4 arbitrary scalar functions φ(τ, x a ), β(τ, x a ), ψ(τ, x a ), γ(τ, x a ) in the following way where τ is the conformal time related to the cosmic one as a(τ )dτ = dt. To avoid gauge fixing issues we use gauge-invariant variables. There exist two independent gauge-invariant variables (the Bardeen potentials), which fully determine the scalar perturbations of the metric tensor [47]- [50] where χ = aβ + a 2γ , ϑ = β + γ , H(τ ) = a /a, and differentiation with respect to the conformal time τ is denoted by a prime. Explicit calculations show that ζ is a gauge-invariant function, and it is therefore not possible to choose such gauge that ζ = 0. Now, the function ζ can be written as 14) where δR GI is the gauge invariant part of the curvature variation: (3.15) and we have used Hereafter k is the comoving wavenumber. The formula (3.10) thus provides us with one equation for Φ and Ψ. It is essential to get another relationship among the Bardeen potentials and one can accomplish this task by considering the (i = j) components of the perturbed ES equations, Eq. (2.4). Trivially, all the terms with δ i j go away, δK i j = 0 and δΓ 0 mn = 0, and only few non-trivial pieces survive. After some algebra we get the following equation ( 3.17) This provides us with the second equation which we need in order to solve for the Bardeen potentials.
It is curious to note that for radiative sources we were able to decouple the evolution equations for the the metric perturbations from the fluctuations in fluid. For the sake of completeness we have provided the equations that determine the fluid fluctuations from the Bardeen potentials in the appendix.

Perturbation summary
In this section we are able to summarize two differential equations for the two Bardeen potentials, Φ and Ψ. They read as follows: Clearly, equations (3.18) and (3.19) provide us with a coupled system of equations, albeit non-local, for the two Bardeen potentials, Φ, Ψ. In our analysis, we never really need the quantities describing the matter fluctuations, but they can indeed be obtained by looking at the other components of the ES equations. We provide some details of that in the appendix.

Perturbations at late-times in the de Sitter limit
Evidently it is a rather difficult task to solve the evolution equations for the Bardeen potentials for the most general case, even if ζ vanish 10 . We can gain however significant insight into how the perturbations evolve within certain approximations. In particular, we are interested in addressing two important questions: (i) do we recover the GR limit during the inflationary phase, this would ensure that the inflationary mechanism of generating near scale-invariant perturbations can be trusted despite the higher derivative modifications to gravity, and (ii) we also want to see whether any pathologies develop around the bounce, this would lead to robustness of the bouncing mechanism, which several other non-singular alternatives seem to lack. In this section we will address the first question.

Localizing equations
At later times, t → +∞, solution (2.19) tends to a de Sitter space-time, with R → 4Λ. This simplifies the perturbation equations (3.10) and (3.17) remarkably. At large times R → −r 2 /r 1 , and therefore we get where the subscript 'dS' indicates that we substitute R B = −r 2 /r 1 . Resulting perturbation equations are as follows: In this limit W is an analytic function with constant coefficients. We can therefore solve the first equation assuming that we know the roots of the algebraic (or transcendental) equation W(ω 2 )| dS (ω 2 − r 1 ) = 0 [37]. In the case of simple roots ω i this reduces to 11 There is always one root ω 2 1 = r 1 but there could be other roots of W(ω 2 ), ω 2 i with i = 2..N . Once δR GI is found to be, say δR i corresponding to the root ω 2 i , we can express δR in terms of Φ and Ψ and solve the equation along with (4.2). The most general solution is obviously an arbitrary superposition of all the solutions. Before proceeding further a couple of remarks are in order. Firstly, the fact that solutions to non-local equations can be approximated as a linear superposition of solutions to a collection of local equations is a rather interesting feature of these models, but it is not surprising and has been discussed in the literature. Secondly, the fact that there always exists a root ω 2 1 = r 1 is also not surprising. A straight forward application of the results obtained in [14] tells us that for actions of the form (2.1), there is always a scalar degree of freedom, the famous Brans-Dicke theory belongs to this class as well, r 1 is simply the mass squared of this scalar. Finally, we note that the roots, their number and values depend on F, and in particular one can choose F's such that there are no roots other than r 1 . 10 There have been steady progress in trying to solve and understand some of the properties of such infinite differential equations, see for instance [51]. 11 If some roots ωi are multiple roots they produce equations of the form ( − ω 2 i ) m i δRGI = 0, where mi is the multiplicity, and these must be treated differently. For the case of double roots see for instance [38]. For simplicity we assume that only simple roots are present.

Trivial case: δR GI = 0
Trivial possibility of Eq. (4.3) implies Φ − Ψ = 0 (4.5) as it is in GR. What we then have to solve is Substituting the appropriate background functions for de Sitter, and R B ≡ 12H 2 , we can solve the above equation in terms of elementary functions where and we remind the readers that the inflationary scale factor in conformal time is given by At large times clearly η → 0 as the scale factor is increasing exponentially, and all the terms die out. This is precisely the mode that one would obtain in ordinary GR in the presence of a cosmological constant and radiative sources. Within a few e-folds, the radiation density and it's fluctuations become irrelevant and the quantum vacuum fluctuations of the inflaton becomes the dominant source of metric perturbations.

Non-trivial case, δR GI = 0
We now consider solutions to the equation Exact solutions to the above equation are known in terms of the conformal time: where J ν , Y ν are Bessel functions of the first and the second kind respectively, and Notice that when all d 1k , d 2k 's vanish we recover the well known GR limit Φ − Ψ = 0. We are primarily interested in the late time super-Hubble limit, k|τ | 1, of these modes. The asymptotic forms of the Bessel functions only depend on the real part of ν, lets call it ν R . Then we have We are now in a position to identify the different cases: • If ω 2 i < 0, then ν 2 is positive, and ν > 3/2. This means the Y ν function grows faster that (−τ ) −3/2 as τ → 0, which in turn implies that Φ − Ψ grows at late times making the inflationary trajectory unstable, and the perturbations unviable for inflationary cosmology.
• If ω 2 i ≈ 0. This is a rather interesting and special case corresponding to the presence of a massless degree of freedom. As it is clear from the asymptotics, in this case Φ − Ψ tends to a constant at late times. In fact, if one assumes the usual Bunch-Davis type subHubble vacuum initial conditions for this field, then one will end up with a scale invariant spectrum for these modes. This is very similar to the curvaton scenario [4] that is routinely considered in inflationary cosmology and we plan to study the phenomenological implications of this mode in more details in the future.
• If ω 2 i is complex, then ν 2 is also complex and whether these modes are dangerous or not depends on the value of ν R . According to the discussion above, for phenomenological consistency of the model we must have ν R ≤ 3/2.
Thus to summarize, our criteria for preserving the inflationary mechanism is that all the roots, ω i , should be such that the corresponding For real roots, this simply means that they must also be non-negative. In particular since ω 2 = r 1 is always a root, we must have r 1 > 0. As discussed in section 2.3, this is consistent with the observations in [13], where it was found that when r 1 < 0, the de Sitter solution is no longer a stable solution, but it is when r 1 > 0. Therefore, there exists a large class of models which are compatible with inflationary predictions.

Perturbations across the bounce
Having determined the late time behavior of the different perturbation modes around our inflationary bouncing solution, our next task is to see how they behave near the bounce.
Since we have already provided the criteria to ensure that none of the fluctuations grow during inflation, all we now have to check is that there are no modes which become singular at any finite time. In particular we will focus near the bounce, as within a few e-folds of expansion the scale-factor is expected to resemble the de Sitter phase which we have already analyzed, at least this is what happens in the hyperbolic cosine bounce (2.19). In fact, for most part we will restrict to this special case to keep the analysis tractable.

Bounce Limit
Near the bounce, even though the equations can again be simplified considerably it is a more complicated limit. The main difference with the previous section is that we can approximate R B by a constant but cannot assume that all its derivatives vanish. We are however, going to look at symmetric bouncing solutions, so that near the bounce we can make the following assumptions about the background quantities: where the subscript b stands for the bounce and numbers indicate the number of derivatives. Under these assumptions ζ ≈ −2R 2 Φ + ( B − r 1 )δR GI (5.2) and equations (3.10) and (3.17) simplify to Note that even if R 2 = 0 equations do not become identical to the ones in the previous section because non-local operators become slightly different due to the new constant R b instead of R dS . It is more difficult do disentangle the above equations but on the other hand the d'Alembertian operator is simpler, thanks to the above mentioned assumptions: 5.5) and δR GI becomes a bit shorter As the first step we solve the homogeneous equation for ζ (5.3). Then we express Φ from (5.4) since the only differential operator acts on Ψ there, and as the last step substitute everything in (5.2). The result will be an equation for Ψ of the fourth order. The main question here is whether there are modes blowing up exactly at the bounce.
In what follows in this Section we are going mainly to specialize to the exact solution (2.19) so that we can use some nice relations among the background quantities. These kinds of relations are not unique to this solution and the results can be easily generalized. Such general considerations however do not shed more light on the problem and the most important features are already captured in this particular solution.
In order to consider the general bouncing solutions, and not just the hyperbolic cosine bounce, we must avoid the specific background relations for particular solutions. Although tedious, it is certainly possible to do so. Indeed, we just need to keep values R b and R 2 general, although they are always related to each other, thanks to our general ansatz, as We stress that the parameters r 1 and r 2 are not initially constrained. The only reason why we are not elaborating on this general case is that this would significantly complicate all the formulae and would not provide any new insights into the problem. There is nothing special about the hyperbolic cosine bounce except the fact that it has a closed analytic form and also obeys some very elegant relations among its parameters.

Trivial case, ζ = 0
The simplest solution to Eq. (5.3) is ζ = 0. This simplifies Eq. (5.4) to and we can deduce Φ from here as It is important to note that from Section 2.3, see Eq. (2.19), R b = 3r 1 , and therefore equations can be separated in the small k 2 limit. Specializing to the known exact solution we have for k 2 = 0Ψ + r 1 Ψ = 0 (5.10) Having r 1 > 0 we get bounded oscillations across the bounce. Then Eq. (5.3) reads where we have used R 2 = −3r 2 1 . These are again two bounded oscillating modes. If k 2 > 0 things are more messy. Again discussing only the exact solution from section 2.3, we have 6Ψ + 4k 2 Ψ + 6r 1 Ψ = 2k 2 Φ (5.12) Using ζ given by Eq. (5.2) we yield for ζ = 0: This equation has four solutions Ψ = c ±± e ±ν ± t , ν ± = 1 √ 6 −(9r 1 + 4k 2 ) ± 9r 2 1 + 4k 4 (5.14) where signs are independent. We see that all the modes pass the bounce point smoothly and in the limit k → 0, modes for Φ and Ψ would mix together giving the above obtained solutions in k = 0 case. What is reassuring is to note that the short wavelength fluctuations (k → ∞) makes the solutions more and more oscillatory rather than contributing to any exponential growth of fluctuations. This is an encouraging sign in terms of whether these kinds of theories can lead to a truly UV complete theory of gravity.
One may worry that the right hand side of eq. (5.9) may become zero for a finite value of k for other more general backgrounds. This however, does not change the analysis of the equations and just shifts the point when perturbations decouple to some finite value of the wavenumber k. Finally, we note that these four modes should merge with the GR mode and the ω 2 = r 1 mode that was discussed in the previous section during the de Sitter phase.

Non-trivial case, ζ = 0
As in the previous section we first solve a simple homogeneous equation of the form 5.16) and the exact behaviour depends on the value of ω i . It is clear, nevertheless, that it passes the bounce point smoothly. Such a non-zero ζ serves as an inhomogeneous source term in right hand side of the two remaining equations. Using Eq. (5.4), we determine Φ as

This equation has a solution
where as before we have used particular relations for the exact solution Eq. (2.19). Substituting the results in Eq. (5.2), we get Extra modes related to the inhomogeneous part will be either oscillatory functions of the frequency determined by ζ i , or resonance modes if frequencies of ζ i and Ψ coincide. We however stress that even the modes which look dangerous as they may grow exponentially do not blow up while passing the bounce point, and therefore must merge eventually with the late time de Sitter modes discussed in the previous section. They may continue to grow, in which case they must be eliminated according to our criteria (4.16), or they must be modes which eventually die out.

Conclusions
In this paper we studied the perturbations of a special class of string-motivated higher derivative extension of Einstein's gravity, which is well-known to ameliorate the Big Bang cosmological singularity of the FLRW metric. More specifically, our action contains the Ricci scalar but is non-local in the sense that it contains higher derivative terms up to all orders. It can yield non-singular bouncing solutions similar to that of the hyperbolic cosine bounce which was first found in [12]. The higher derivatives lead to weakening of the gravity in the UV regime which facilitates the bounce for the FLRW metric in presence of non-ghost like radiation and the cosmological constant. For the first time we analyzed equations for the scalar perturbations for such non-local higher derivative theories of gravity in complete generality. We provided the general nonlocal equations that the two Bardeen potentials satisfy, and studied the evolution of the generic fluctuations in the relativistic fluid and the scalar potentials during the bounce as well as at late times corresponding to the de Sitter attractor solution. Although we have not included the dynamics and perturbations of the inflaton field, our analysis assessed the robustness of the inflationary scenario in terms of perturbations. We were able to arrive at a simple criteria such that the modes (small and large k ) pass the bounce phase smoothly, and consequently decay during the de Sitter phase. This means, at least at the classical level there are no pathologies present during the evolution. We also noted the possibility of realizing a curvaton scenario given the extra scalar degree of freedom present in these models are light.
One of the issues we haven't discussed in our paper is anisotropic deviations from the bouncing background. This requires looking at homogeneous (time dependent) but anisotropic Bianchi I type of models. (The isotropic homogeneous mode was already considered in [13] and issues of stability were discussed.) This is an important issue to determine the robustness of the bouncing mechanism as bouncing/cyclic models are notoriously plagued with Mixmaster type chaotic behavior related to the growth of anisotropies in a contracting universe. The generalization of the application of the ansatz (2.8) that has been performed in this paper to arbitrary metrics (potentially inhomogeneous and anisotropic) should help in addressing this issue, but we leave this for future considerations. For stability analysis related to anisotropies in other non-local gravity and string-inspired models one can see [25,52] and references therein. It is probably worth pointing out though that anisotropies are expected to be the largest at the bounce point, and thus if there is a "smooth-enough" bouncing patch, anisotropies are only going to decrease at large times in the past or future.
For simplicity, some of our analysis was done for the hyperbolic cosine bounce, but we discussed how the results should generalize to other bouncing solutions. This strongly suggests that higher derivative actions of gravity considered in this paper do indeed yield a non-singular bounce without a pathology, and in doing so can provide a model of inflation, both in the past and the future.
The resulting expression for ε then is One interesting special case is R B → −3r 1 , but this is not a limit for the known background solutions.