$R^2$ inflation to probe non-perturbative quantum gravity

It is natural to expect a consistent inflationary model of the very early Universe to be an effective theory of quantum gravity, at least at energies much less than the Planck one. For the moment, $R+R^2$, or shortly $R^2$, inflation is the most successful in accounting for the latest CMB data from the PLANCK satellite and other experiments. Moreover, recently it was shown to be ultra-violet (UV) complete via an embedding into an analytic infinite derivative (AID) non-local gravity. In this paper, we derive a most general theory of gravity that contributes to perturbed linear equations of motion around maximally symmetric space-times. We show that such a theory is quadratic in the Ricci scalar and the Weyl tensor with AID operators along with the Einstein-Hilbert term and possibly a cosmological constant. We explicitly demonstrate that introduction of the Ricci tensor squared term is redundant. Working in this quadratic AID gravity framework without a cosmological term we prove that for a specified class of space homogeneous space-times, a space of solutions to the equations of motion is identical to the space of backgrounds in a local $R^2$ model. We further compute the full second order perturbed action around any background belonging to that class. We proceed by extracting the key inflationary parameters of our model such as a spectral index ($n_s$), a tensor-to-scalar ratio ($r$) and a tensor tilt ($n_t$). It appears that $n_s$ remains the same as in the local $R^2$ inflation in the leading slow-roll approximation, while $r$ and $n_t$ get modified due to modification of the tensor power spectrum. This class of models allows for any value of $r<0.07$ with a modified consistency relation which can be fixed by future observations of primordial $B$-modes of the CMB polarization. This makes the UV complete $R^2$ gravity a natural target for future CMB probes.


Introduction
Finding a gravity theory consistent with the concepts of quantum field theory is a longstanding problem. General Relativity (GR) [1] was known not to be ultra-violet (UV) complete from the very beginning. Hence one is forced to modify GR in order to construct any self-consistent model of quantum gravity. Moreover, generalizing GR one has not give it up altogether, as it is heavily supported by absolutely all measurements in the low energy or infra-red (IR) regime including the recent direct discovery of gravitational waves [2].
One of the most obvious and at the same time very promising generalization to consider is M 2 P 2 R + R 2 /(6M 2 ) Lagrangian instead of just M 2 P 2 R as in GR where as usually R is the Ricci scalar, M −2 P = 8πG, G is the Newtonian constant, c = 1 and M becomes the mass of what is the propagating scalar in this model, dubbed scalaron. We refer it hereafter as the local R 2 model or just R 2 model (or Lagrangian, etc.). Besides being the simplest one and having only one free parameter M which value is fixed to M ≈ 1.3 × 10 −5 M P by the observed Fourier power spectrum P R of primordial scalar (matter density) perturbations in the Universe, this generalization has two major advantages. First it was proven in [3,4] that it is renormalizable, i.e. UV complete in the scalar sector. Second, a dramatically successful model of inflation [5][6][7] known as "Starobinsky inflation" is provided by the R 2 Lagrangian. With the latest Cosmic Microwave Background (CMB) measurements by the PLANCK mission [8,9] and the more recent BICEP2/Keck Array experiments [10], R 2 inflation produces an excellent fit for the key inflationary predictions.
The advantage of renormalizability of the local R 2 gravity has an unfortunate fate to be spoiled by the non-unitarity as a spin-2 ghost appears in the physical spectrum. This ghost is a manifestation of the Ostrogradski instability [11] due to higher derivatives and it appears as long as the required for the full renormalization W 2 term, with W being the Weyl tensor, is included. It is not hopeless to try to remove such a ghost and few ways are known in principle. Ghosts may become unphysical in constrained systems [12,13]. Also one can try to consider special constructions like Horndeski theories [14] in which higher derivatives in the action still result in a second order equations of motion (EOM). Another way is to promote the Lagrangian to a non-local model such that infinitely many derivatives form some operator which does not create new poles in the propagator and consequently does not generate new physical degrees of freedom. On this way possible operators which we may encounter are: analytic in derivatives like exp( ) with being the covariant d'Alembertian operator, non-analytic in derivatives like 1/ , having logarithms like log( ), etc.
It was shown already in [15] that a systematic accounting of one-loop corrections from quantum matter fields to the R 2 gravity leads to infinite derivative logarithmic functions of the d'Alembertian in the action. Theories with analytic infinite derivative (AID) operators in the action naturally appear in string theory when the string field theory (SFT) [16,17] is considered. Also p-adic string theory [18] is an example of a model featuring AID Lagrangians. Both of these stringy models are unitary and UV complete theories. A study of gravity theories having similar AID operators was initiated in [19]. This led recently to an intensive study of AID gravity theories [20,21] which were shown moreover to be easily made ghost free by adjusting the AID operators. This study was focused on a quadratic in R Lagrangian. Note that absence of ghosts in such a setup can be achieved actually only by introducing an infinite number of derivatives, i.e. non-local operators. Further questions of renormalizability [22], presence of a ghost-free bounce [23] and an amelioration of singularities [21,[24][25][26][27] were addressed emphasizing in all instances a possibility to resolve successfully and consistently the problems in the framework of a quadratic in curvatures AID gravity modification (AID quadratic gravities/models/theories/etc. in short). On top of this ad-hoc AID scalar field models with a minimal coupling to GR were proven to be interesting in tackling the Dark Energy problem [28][29][30].
There are several reasons to stick exactly with analytic differential operators. Primar-ily these are the presence of a well defined low-energy limit and their native appearance in a more fundamental approach which is SFT. However, as have been already mentioned the presence of exactly infinite number of derivatives (i.e. non-local operators) is a requirement to avoid ghosts. Interest in exactly quadratic in curvatures Lagrangians is stemmed from the fact that in many applications it is enough to study highly symmetric backgrounds, in particular Maximally Symmetric Space-times (MSS, which are in fact (Anti)-de Sitter ((A)dS) and Minkowski in the space-time dimension greater than 2) and linear perturbations around them. Linear perturbations in turn are described by the quadratic variation of the action. It was proven by explicit construction in [31], that starting with a very generic action for the metric field such that the Lagrangian is analytic in curvatures and covariant derivatives, and focusing on the task of studying linear perturbations around MSS, one ends up with a quadratic in curvature action with analytic functions of the covariant d'Alembertian operator. In the most general case these analytic functions of derivatives become AID operators. No other combinations of derivatives apart from d'Alembertian and its AID functions appear. This is exactly the AID quadratic gravity and this is the most general and the only relevant Lagrangian we need to use in studying fluctuations around MSS.
The full gravity theory does not have to be just quadratic in curvatures. The point is that only the quadratic in curvatures part of some more general theory is responsible for the structure of propagators. This structure in turn in vast amount of situations determines whether the theory is unitary or not. However, another crucial property of a theory, its renormalizability, may require higher curvature terms present in the action [32][33][34]. Nevertheless the already observed properties of AID quadratic models make it inevitable to ask whether these theories are capable to eventually grow up to a full non-perturbative quantum theory of gravity. A significant step in this direction with a positive outcome was made recently in [35] where super-renormalizable or finite quantum gravity candidates around MSS are constructed.
Given the success of AID quadratic gravity, it is natural to study whether it can admit inflationary solutions for some range of curvature. This is because cosmic inflation is not only a very successful theory of the early Universe [5,36,37] but also at the same time for the moment is the best test-bed to challenge modified gravity theories. Viable models of inflation which can be parametrized by a number of free dimensionless parameters which values have to be fixed from observational data produce definite predictions about postinflationary space-time metric perturbations given that an inflationary stage lasts long enough. The simplest models like the Starobinsky one have only one such parameter, so their predictive power is high. Further it explains the emergence of the Standard Model physics through the reheating mechanism [38,39]. However, in spite of a very large number of other inflationary models already excluded by observations, still there remains a sufficient amount of them which remain viable, too, see e.g. [40,41]. In many of these models having a larger number of free parameters, additional scalar fields are introduced and gravity remains Einsteinian up to the Planck curvature (unity in our notations). Even if we restrict ourselves, as we do in this paper, to purely geometrical models of inflation and modified gravity, still even in the R 2 model observations can probe a large curvature regime only up to R = 4N M 2 ∼ 240M 2 where N is the number of e-folds from the end of inflation back in time, that is much less than unity. In this regard AID gravity theories become natural candidates accounting the fact that they can be made ghost-free and tend to be renormalizable.
To have inflation is equivalent to guarantee a presence of solutions with a long enough nearly dS expansion and a subsequent graceful exit from this regime. To compute the key inflationary parameters one has to study linear perturbations around this nearly dS expansion phase. As we have explained above, AID quadratic gravity action is the maximal possible generalization one should ever consider for this purpose. This, as it was proven, covers considerations of inflation in arbitrary general original gravity theory as long as its action is analytic in curvatures and derivatives, and an appropriate inflationary solution exist. One equally can maintain a structural connection with other theories, like SFT, while this is not obviously necessary.
A first and successful try of embedding R 2 inflation into quadratic in R AID gravity was performed in [42]. In a more recent paper [43] it was argued that a particular quadratic in R and in W Lagrangian with AID operators is a renormalizable at least by power-counting and ghost-free gravity theory. The local R 2 inflation can be seamlessly embedded in this AID quadratic Lagrangian as well. Parameters of the new model allow to maintain a good agreement with the observational data at easy.
The main purpose of this paper is to deepen from the side of inflation the study of the AID quadratic gravity model undertaken in [43]. In what follows we will provide more support for the particular action used in that paper. The advance of the current analysis is the proof that the AID quadratic action which was derived in [31] as the least non-redundant action for studying linear perturbations around MSS is in fact redundant thanks to Bianchi identities. Similar ideas were used in [44] doing computations around the Minkowski background. Here we provide the full treatment around MSS and this is the purpose of Section 2. Further we prove a very important statement that under certain assumptions the space of background space-homogeneous solutions in our AID model is identical to the space of backgrounds of a local R 2 gravity. This is a very important step since it allows to claim that the classical inflation remains an attractor behavior in the case of AID quadratic gravity. This is done in Section 3. As the main accent of the present paper we systematically derive the inflationary parameters following from our model keeping the leading order in the slow-roll approximation throughout the whole computation. In particular we compute spectral tilts and tensor to scalar ratio. Note, that the previous studies assumed an exact dS background in the course of computation and applied the slow-roll approximation only starting from the action for canonical perturbation variable. The technique developed in this paper opens ways to restrict tighter the parameters of the new theory and to meet more and more toughly squeezed observational constraints. All the inflationary computations related to our model are accumulated in Section 4. In Section 5 we discuss the main results obtained in the paper and outline open questions. At last, extensive Appendices contain all the notations used in the paper as well as most technical pieces of the derivations.

Most general AID quadratic gravity action around MSS
One of the most crucial results of [31] provides a most generic action for studying linear perturbations around MSS. Consider the following action (all notations without an explicit explanation in the main text hereafter are accumulated in the Appendixes): where λ is a dimensionless constant which is convenient to control the magnitude of the R 2 modification and Λ is an in principle possible cosmological constant term. Briefly the reduction is done by carefully accounting all possible terms which may contribute nontrivially to the second variation around MSS (and dropping all other terms). The fact which is heavily used on this way is that all curvature tensors on MSS are annihilated by covariant derivatives. An important assumption essential for the actual computations and which was discussed in the Introduction is that all functions F are analytic. To be precise we need at the moment to have these function analytic around zero. This is indeed required from the physical point of view. We want functions F( ) reduce to constants or vanish in a low-energy limit because we have to restore GR at very low energies. There is also other way to understand this. In writing F( ) we always assume that there is an energy scale of the gravity modification M, which we name the scale of non-locality as in principle we may have infinite derivative operators (M should not be mixed with the much lower energy scale M at which the R 2 /6M 2 term in the local R 2 inflationary model becomes comparable to the GR term R). This scale enters as F( /M 2 ). Even though for most of our technical steps we can put M = 1 we still want to have a local or trivial limit once M → ∞ in order to be able to eventually restore GR. Hence, we come to the conclusion that functions F must be analytic at least in the origin.
Proposition: Action (2.2) is redundant in describing linear fluctuations around MSS. This proposition can be proven to be true because the previous analysis did not make use of Bianchi identities which is the cornerstone of the succeeding further reduction. To start with, action (2.2) can be rewritten as 3) The purpose of using the Weyl tensor W and L-tensor is their simplicity. Both are identically zero on MSS. Moreover, W is zero on any conformally flat background (including spatially flat FLRW). 1 The term to be attacked by Bianchi identities is the L-piece. A good reason to tackle this term somehow is that being simple on the background it produces tremendous complications while trying to compute perturbations.
To make the long story short we put all the technical details of the reduction into Appendix B. Upon a lengthy but a straightforward procedure the full resulting action relevant for study of linear perturbations around MSS vacua of (2.1) becomes We consider this action as a significant simplification of (2.3) for several reasons: 2 (i) it contains only Ricci scalar and Weyl tensor and no Ricci tensor or its linear combination with the metric. Weyl tensor enters only quadratically and being identically zero on any conformally flat manifold does not contribute to conformally flat background solutions. Importantly, spatially flat FLRW metric is conformally flat.
(ii) as such, any solution already found in the literature with only RF R ( )R piece in the action is a solution to equations of motion which one can derive from our new action.
(iii) linear perturbations of Weyl tensor are very simple using (1 + 3) decomposition of the ADM formalism. These were computed in [43] and one can track computations relevant to our AID models in application to inflation to the end. Actually, perturbations of a possible term with any of the second rank tensors (Ricci, Schouten, Einstein or L-tensor) turn out to be very much complicated and seem to be intractable. It is worth stressing that actions (2.2) and (2.3) are not fully equivalent. They are equivalent as long as at most linear perturbations around MSS are considered. As a consequence, non-MSS may be solutions to EOM derived from one action and not from another. For example, local R 2 inflationary background is a solution to EOM derived from action (2.3) and is not a solution as long as the quadratic term with a second-rank tensor is We still consider using L-tensor is preferred as it is identically zero on MSS. 3 Classical dynamics of AID quadratic gravity without Λ

EOM and solution construction
The major focus of the present paper is on the consideration of inflation in AID quadratic gravity and for that purpose the cosmological term is not needed. It is shown already in [5] that in fact a cosmological term spoils the inflation. To proceed with actual computation we recite action (2.4) dropping the cosmological term Λ: This action was studied and many technical details were elaborated in [23,42,43]. We are going to use them without extensive further referencing. EOM which one can derive from action (3.1) (see [46]) read: The trace equation reads Terms linear in Weyl tensor are not present in the trace equation because their trace vanishes by construction on any space-time. Terms O(W 2 ) can be found in [47]. We are interested in cosmological solutions of the spatially flat FLRW type. First this implies that the Weyl tensor vanishes and as such it does not manifest itself in the trace equation neither in the background nor in linear perturbations. Second, such solutions for the metric are space-homogeneous and isotropic. This means that system of equation (3.2) has essentially two distinct equations. The standard choice is the trace equation and the ( 0 0 )-equation. However, presence of Bianchi identities guarantees that given we have a solution to the trace equation with zero RHS then it will be a solution to the whole system of equations modulo a possible radiation source (which is conserved and is traceless). We are thus focused on solving the trace equation (3.3) which is a non-linear differential (nonlocal) equation on R and all the differential operators are of the form of d'Alembertian.
We start solving the trace equation by reminding that originally it was proposed in [20] to use an ansatz to construct solutions. First we note that the original ansatz also had a free constant term r 2 in the right hand side but it is not compatible with the absence of the cosmological term. Also we note that this ansatz was indeed helpful to construct several exact solutions to equation of motion.
It is instructive to show sketchy how the technique works. Substituting (3.4) into (3.3) we restore the result obtained in [20] The way to solve the latter equation is to assume the algebraic conditions Since we constrain here only parameters we get the most what we can using (3.4). If we do not impose the above conditions then we must satisfy additional conditions on R which can be shown to trivialize possible solutions to just one R = 0. We accumulate the details supporting this claim in Appendix C. This will become useful in the coming Subsection.

Proof that (3.4) is general solution to (3.3)
Now we formulate the main claim of this Section and in fact a very important statement for the development of AID quadratic gravity theories in general. Let us start with noting that having a physical attitude to the problem we formulate here sufficient and not obligatory necessary conditions.
The first condition (i) just serves for the setting of the present paper to discuss a space-homogeneous inflationary space-time and is simple to account. Technically during the proof the only property to be exploited will be the space-homogeneity of the metric in the synchronous frame. As such the proof itself can be applied to more general spacehomogeneous metrics. For example to anisotropic backgrounds like Bianchi I or other. However, we need the metric to be conformally flat to eliminate the Weyl tensor squared terms in the trace equation (3.3). This would allow us to claim that the restrictions imposed by this proposition on the space of solutions apply to the full trace equation. Only for this purpose we stick to spatially flat FLRW metrics only. This implies that given the Weyl tensor dependent term is not included in the action (3.1) one can relax condition (i) to: (i) the metric is space-homogeneous in the synchronous frame.
The second condition (ii) needs more explanations though. We name Fourier harmonics the eigenfunctions of the d'Alembertian operator such that where w i are constants. Generically we expect the spectrum of the d'Alembertian is continuous even though this is not crucial. We name ϕ i the Fourier harmonics in analogy with the flat space-time where they reduce to the plane waves which in turn are used do define the Fourier transform. A crucial property of the Fourier transform in the flat space-time is that the corresponding harmonics form a basis in the domain of square integrable functions L 2 . Or in other words any square integrable function can be presented as a linear superposition of plane waves. In our model the situation seems to be more involved as a priori these nice properties of the Fourier transform in the flat space-time cannot be elevated to a curved background. It is known from the spectral analysis of the Beltrami-Laplace (BI) operator on Riemannian manifolds that indeed the eigenmodes of the BI operator form a basis in L 2 as long as the manifold is compact or has a boundary [48]. In most cosmological applications the space-time manifold is however pseudo-Riemannian (i.e. the metric is not positively defined and d'Alembertian operator replaces the BI operator), non-compact and without a boundary. In this situation general theorems do not help and presently one has to consider systems case-by-case. Paper [49] provides an explicit proof that in two notable cases of dS and (A)dS space-times indeed the eigenmodes of the d'Alembertian operator form a basis for square integrable functions. This remains valid in a special situation when there is no spatial dependence is present. It is an important situation though since in the vast majority cosmologically viable backgrounds are space-homogeneous. Naturally, it is the case for the present paper as well. Technically, this implies that the d'Alembertian operator lacks of spatial derivatives and eigenmodes ϕ i depend on time only.
Coming to physical grounds we stress that the regime of the space-time evolution of interest in the present paper is the nearly dS expansion. This in combination with the results in [49] provides some hint that our condition (ii) in the proposition above is sensible. However, there is one more physically important argument why a physically viable space-time must have such a structure that the d'Alembertian operator eigenmodes form a basis. Namely, we expect that our model can be quantized. To have this happen we have to have a vacuum and creation and annihilation operators which in the canonical quantization scheme appear as operator coefficient in front of Fourier modes in which a given classical solution is decomposed. Given a situation that Fourier modes do not form a basis (i.e. the set of modes is not enough to represent any function) we will hit a problem that certain classical configurations cannot be quantized in a canonical way. This simple consideration shows that the fact that eigenmodes of the d'Alembertian operator form a basis is necessary to implement the canonical quantization scheme. This gives us even a stronger hint that we indeed want the condition (ii) in the proposition to be satisfied.
Finally, we do not specify explicitly the domain of functions on which the completeness of the Fourier decomposition is true. We presume that in most cases we need to have it either for functions from L 2 or functions with a compact support which is a more plausible case as long as time evolution of a classical system is considered. This will be noted just below as well.
Therefore, in proving the proposition we assume that the scalar curvature R as any function can be represented as and w i are constants. Possible constants in front of ϕ i in the decomposition of R are absorbed inside of ϕ i for simplicity. Few comments are in order here. First, one should not be confused with the fact that R itself depends on the metric as for the time being it is just some function of time. Second, one should not worry about possible non-trivial asymptotics of R in past or future infinities (which may render it non-square integrable) and consider only a given time interval during which our model describes the evolution of the Universe. This will waive doubts of the square integrability since the function gains the compact support by construction. In other words it is equivalent to saying that we work in a given coordinate patch. Also, here we explicitly come to the special situation mentioned above that all functions depend on time only since a space-homogeneous background is considered. The corresponding simplification will become crucial to fulfill the proof. Using (3.8) one readily computes and further Notice that for i = j we have using the Taylor series expansion ω ii = F R (w i ) where the superscript (1) denotes the derivative with respect to an argument. Substituting all of that into (3.3) and accounting that the Weyl tensor vanishes one yields To prove the proposition we have to show that no (non-trivial) solutions to (3.11) exist as long as R is a superposition of more than a single Fourier eigenmode. First we note that the technique of equating coefficient to zero does not work in this general case. Indeed, the quadratic in ϕ i term in (3.11) can be eliminated by requiring Since however different ϕ k are eigenfunctions of d'Alembertian with different eigenvalues they are linearly independent. This means that in order to satisfy the latter equation we must require M 2 P − 6λF R (w 1 )w k = 0 for each k and as such all w k are equal. We thus effectively come back to the situation R = ϕ 1 like it is served by (3.4).
Thus we must keep the quadratic terms in (3.11) and solve it as a differential equation on ϕ i . Satisfying (3.11) will necessarily produce stringent constraints since the resulting solution for R must be identical to the Ricci scalar constructed from the metric. Note, that in the beginning of the proof we have mentioned that R is just some function of time.
Here we explicitly make reference to its relation to the metric. This, however, in no way complicates the use of desired spectral properties of the d'Alembertian.
Going further one can pass to modified quantitiesφ i = ϕ i + c i where we have done shifts by constants defined as Hence we rewrite (3.11) as where we have used the fact that all ϕ i are space-homogeneous and depend only on time.
The sign change in front of ∼φ 2 i is due to the signature of the metric and also we assume that g 00 = −1. Also a common factor λ has been cancelled. Interestingly, we recognize in the latter formula the conserved integral of energy originating from a sigma-model-type dynamical system.
To make the succeeding analysis more transparent we rewrite the last formula using matrix notations as follows˙ where R is a vector made ofφ i and w = diag(w 1 , w 2 , . . .) and ω is a matrix formed by ω ij .
We use a simple transposition as all quantities are real and matrices are symmetric from the physical origin of the problem. We diagonalize matrix ω by choosing an appropriate matrix D. We can always do this because if ω cannot be diagonalized then some values w i are identical and we must just drop equivalent terms from decomposition (3.8). Denoting d 2 = D T ωD and using further redefined functions Q = dD T R we get a canonically normalized diagonal term with derivatives. The whole expression transforms aṡ We can simplify the things even more by diagonalizing the matrix ν choosing an appropriate matrix M . Denoting m 2 = M T νM and Here the most crucial achievement that matrix m is diagonal. Differentiating the latter equation with respect to the time t one getṡ As noticed above, all ϕ i are linearly independent and the same are P i . Indeed, since the matrices which define the quadratic form are non-degenerate this guarantees that P i are linearly independent. We thus can consider only the second order linear equations in the latter expression as all of them must be satisfied independently. Since moreover matrix m is diagonal we readily find each P i as where we have assumed m = diag(m 1 , m 2 , . . . ). Returning to (3.8) we can rewrite the corresponding expression in matrix notations as well. That isR Note that the latter equation is valid for any space-homogeneous metric as long as g 00 = −1 and as such H is the Hubble function only in the case of a spatially flat FLRW metric. Passing to variables P i we getP To prove the proposition we must show that solutions (3.18) are incompatible with (3.20) for more than one component vector P . Being lucky that we could construct solutions for P i explicitly we just substitute them into (3.20). The resulting expression is where P ± = diag(P i± e ±m i t ). Each component P i is a different exponent and in order to satisfy the latter equation we must put to zero coefficients in front of each of them. Moreover, we must have the constant term on the right hand side vanish. If H = 0 we essentially must require the matrix m to be zero and this is equivalent to having all w i = 0 and as such we come back to the situation R = 0 which is just a sub-case of (3.4) and in no way requires any more general form of R then a single Fourier mode. This completes the proof of the proposition during which we actually have never used that the metric must be exactly of a spatially flat FLRW type.
In a slightly exotic situation such that there is a space-homogeneous metric which generates vanishing factor H in (3.19) in a combination with a non-constant R one need to have a separate consideration regarding the space of solutions to EOM.

Discussion on classical dynamics
We just have proven a very important fact related to theories of type (3.1): all spacehomogeneous conformally flat background solutions are subject to equation (3.4) in combination with conditions (3.6).
To understand what happens we must examine conditions (3.6) which tell that a nontrivial solution (i.e. a solution more involved than a constant R) exists only if there is a point r 1 such that function F R (r 1 ) being considered as a function of r 1 as its parameter obeys two independent algebraic conditions. Moreover, a would be solution must obey an equation which can be derived from a local R 2 gravity. It was elaborated in [43] what a Lagrangian for a local model must be written such that its equation of motion yields R = r 1 R. So essentially we should worry whether function F R (r 1 ) provides a chance to have at least some solution.
The other point of view is to consider presence of a solution as a criterion for function F R (r 1 ) such that it provides a choice of points r 1 at which conditions (3.6) are true. Since functions F X ( ) are not constrained so far apart from being analytic at the origin one may wonder about the space of solutions. Indeed, it is possible that a generic function F R (r 1 ) has many or may be even infinitely many points in which r 1 F 1 is the same value and plus to this F (1) R (r 1 ) = 0 at these points. This is some sense would mean that our model include multiple copies of local R 2 gravity.
Even though mathematically possible we are going to remind that the operator functions F X ( ) get severely constrained by demand that no new excitations must appear in the spectrum. Indeed, as it was derived in all the details in [45] the quadratic Lagrangian for the spin-0 mode of the metric around the Minkowski space-time (where we must fix the operator functions) is In order to contain the spectrum of excitations and to have inflation we must require that the expression in curly brackets has exactly one zero which would correspond to the scalaron. This can be achieved by demanding Here σ( ) must be an entire function. For the definiteness we assume that σ(0) = 0 and in order not to lose generality we introduce σ 0 . Recall that F(0) = f 0 and evaluating left and right sides of (3.23) at¯ = 0 we get 1 = −σ 0 M 2 yielding σ 0 = −1/M 2 . Next, evaluating (3.23) at¯ = r 1 and accounting (3.6) we get This implies that r 1 = M 2 . Differentiating (3.23) with respect to the d'Alembertian, evaluating the result at¯ = r 1 and accounting (3.6) one gets This implies σ(r 1 ) = 0. The above results confirm the derivation done in [43]. However the more important observation is that the condition r 1 = M 2 together with the demand that only one excitation can exist guarantees that from the point of view of the physics of our model only a single unique point r 1 is allowed. This is a very powerful statement because it implies that as long as space-homogeneous and conformally flat metrics are considered our quadratic AID gravity has exactly the same space of solutions as a local R 2 gravity. In particular this means that the inflationary background will remain an attractor behaviour as it was found originally for the local R 2 model in [5].
To conclude this Section we notice that there are already mentioned limitations of our analysis: we are talking about space-homogeneous and conformally flat solutions (while allowing anisotropic metrics like of the Bianchi I types in the absence of the Weyl tensor term from the very beginning), we do not have matter sources apart from perhaps radiation. Also, we do not have the cosmological term in the action. It was mentioned above that the presence of a cosmological term requires (3.4) to be modified as follows In this case, for instance, the whole proof of the proposition from the previous Subsection must be reconsidered if at all possible. We keep the generalization of our analysis to the models with the cosmological term as well as other interesting directions of further developments for future projects.

General considerations
As it was understood in Section 2 any arbitrary analytic in curvatures and derivatives gravity action (2.1) is governed around MSS in D = 4 by action (2.4) or (3.1) if no cosmological term is present. We note that the cosmological term is not needed at all for the inflation to happen. Inflation is the nearly dS expansion phase and thus action (3.1) is the correct starting point to compute inflationary parameters. Therefore having computed and confronted with measurements parameters of inflation from action (3.1) one can further try to approach the full theory of gravity.
As it was further crucially proven in Sections 3, given one seeks for space-homogeneous and conformally flat metrics, our theory features only solutions which can be found in a local R 2 gravity. Technically this means that any solution must obey equation (3.4) and due to the ghost-free conditions only one non-zero r 1 parameter is allowed. Hence, the local inflationary background is seamlessly embedded in our most general consideration. We stress that this particular background is the minimal choice for the physically viable cosmological inflation which must be a long enough nearly dS phase of expansion with a graceful exit. Thanks to the proof of absence of other solutions in our model we guarantee that the local R 2 inflation remains an attractor solution which is an important property from the physical point of view.
Surely, perturbations are expected to be different from a local gravity. Since the main results in consideration of any inflationary model are related to perturbations we focus on analyzing them as detailed as possible and this is the main focus of this Section. Note that the best done so far computations of inflationary parameters [43,50] treated the background as exact dS until the action for canonical perturbation variables. In the present analysis we are going to maintain consistently the leading slow-roll approximation throughout all the computations. To control the slow-roll approximation we use the pretty much standard slow-roll parameter ǫ defined as Its use is elaborated in Appendix D.
Thus, technically we should work with the solution to equation (3.4) which was solved already in [5]. The inflationary phase of this solution is known approximately [5,43] and one can check that the slow-roll approximation is valid for it for appropriate choice of parameters. Here t s is the time of the end of inflation. Also, in principle, the solution to the whole system of Einstein equations may require radiation ρ r . Its amount is given by (see [43]): and in the leading slow-roll approximation the radiation source vanishes (see Appendix D). The latter consideration of energy density of radiation emphasizes that given we do not consider its perturbations we are limited by the linear order in the slow-roll parameter in our analysis as any further expansion would require to include perturbations of the radiation fluid into consideration.

Action for perturbations
In order to analyze perturbations and their properties one can either analyze the linear variation of EOM or an action for perturbations which is the second order variation of the background action. While the results must be the same irrespectively of the approach some steps may be more simple in either of them. Linearization of EOM worked well in previous papers but a construction of the second order variation of the action was possible only around MSS. Below we present for the first time the second order variation of (3.1) around any solution satisfying (3.4) and conditions (3.6). To do so we introduce an auxiliary local action The answer for the second variation of (3.1) turns out to be unexpectedly simple and it reads after some laborious steps outlined in Appendix E where δ local = λF 1 2 2(R + 3r 1 )δ 0 −R 2 δ g + (δR) 2 is exactly the second order variation coming from a local action (4.4). Quantities δ 0 = (δ 2 √ −gR)/ √ −g and δ g = (δ 2 √ −g)/ √ −g are computed explicitly in Appendix E. Further, This is essentially a variation of (3.4) and it would be identically zero in a local case but may not be assumed as such as (3.4) is not an EOM in our AID action. ζ actually becomes an essential quantity in our model. Finally, is an operator analytic in¯ as can be shown using the Taylor series expansion in combination with (3.6). As usual quantities with bars over them are the background ones. First we note that the latter action for perturbations is derived without any assumptions on the background apart from the fact that it satisfies (3.4) and conditions (3.6) are satisfied. It is a general action valid in all regimes and not only around the actual dS expansion of the Universe. This is a complete expression for all modes: scalars, vectors and tensors. The action differs from a known answer for a local R 2 inflationary model by just one term square in ζ. This term is actually non-local due to the presence of the operator Z 2 ( ). The last term containing the Weyl tensor variations is that simple because the Weyl tensor is zero on the backgrounds of interest and as such nothing else can survive upon the second variation.
We proceed by considering scalar and tensor perturbations meaning that the classification is with respect to the 3-dimensional symmetry group.

Scalar perturbations
The perturbed line element for scalar perturbations in terms of Bardeen potentials (Φ, Ψ) reads The gauge invariant perturbation of the scalar curvature is given by This is used instead of δR in the definition of ζ in (4.6) as long as we pass to gauge invariant variables. As it is shown above ζ essentially measures the difference of our model from the local R 2 gravity (see action (4.5)). Also we recall that being a variation of (3.4) it is zero in a local R 2 gravity but is not obligatory trivial in our case. However, ζ is governed by a linear and homogeneous though non-local equation (F.1) which we recite here for the completeness where We are going to explore the solutions to the above equation in the leading slow-roll approximation. This generalizes our consideration in [43] because there we have taken the pure dS background for inflation which is the zero level approximation. Strictly speaking however, having the background (4.2) up to the leading slow-roll correction we have to follow the same approximation in computing perturbations.
Using the details accumulated in Appendix F we state the essence of this consideration that in the leading slow-roll approximation the only possibility for ζ is a trivial solution and δW does not contain scalar perturbations at all. This has a major consequence. Just looking at action (4.5) we see that the scalar perturbations during inflation inflation at the leading order are always governed by a local R 2 action irrespectively of what the full gravity theory is. Technically this is seen through the simplification of all the system of linearized scalar perturbation equations presented in Appendix F. Effectively it behaves the same way as in the pure dS leaving the analysis and the results of [43] intact. As a particular important consequence we gain Φ + Ψ = 0 . (4.12) One can see from the explicit expressions in Appendix F that δW depends solely on Φ + Ψ and as such the contribution which would come from the Weyl tensor piece vanishes.
Having said this we can straightforwardly utilize the results of [42,43] to write down the action for a canonical variable which is Here Υ is the canonical variable in question related to Bardeen potentials as 2F 1R Ψ = Υ. The operator W(¯ ) is W(¯ ) = 3F R ¯ + R + 3r 1 Z 1 ¯ . (4.14) Comparing operator W with the expression for the spin-0 propagator around dS background found in [45] we see that for a consistent theory (around the dS background which is the case during inflation) we should demand for some entire function γ(¯ ).
One may wonder about the denominator F R (¯ ). Naively one would expect that the whole fraction W(¯ )/F R (¯ ) must be an exponent of some entire function to avoid extra poles in the propagator. It is however not always necessary as the denominator would contribute to poles of the propagator only if it has poles on its own. In a particular case that the denominator itself is an entire function or in a situation that would be propagator poles are beyond the domain of validity of our effective theory (in simple words the excitations are heavier than the effective theory scale, in our case heavier than the scale of inflation) one should not worry about the presence of the denominator in the operator function. The detailed analysis of such a situation is presented in [35].
Going further we should proceed with the quantization of perturbations and an evaluation of the two-point function in order to deduce the power spectrum and the corresponding scalar spectral tilt. What is intriguing however, accounting the fact that we have to compute the final quantities at the position of the pole for the canonical variable, i.e. at¯ = r 1 and using that Z 1 (r 1 ) = 0 we are going to get answers identical to those in a local R 2 theory as W(r 1 )/F R (r 1 ) = 3. On the other hand this is not a surprise as in the situation when ζ = δW = 0 we see that action (4.5) is nothing but a second order variation of an action for a local R 2 gravity.
Actual results for the scalar power spectrum P R for R = Ψ + HδR GI /Ṙ and the corresponding scalar tile n s can be found in [42,43] and are as follows where N is the number of e-folds. We thus double argued: first, by using the explicit action for the scalar perturbations and, second, by the rederivation of the action for a canonical variable in the scalar sector that irrespectively of what is the general full gravity theory, an inflation would lead always to the same universal predictions in the scalar sector up to the leading order in the slow-roll correction. One however would get absolutely new corrections coming truly from the nonlocal operators as long as next to leading orders in the slow-roll approximation or higher, i.e. three-or more, -point correlation functions are considered.

Tensor perturbations
Computation of the tensor perturbations was done in [43] and already accounts the leading slow-roll approximation. The action for the canonical variable is where h µν is transverse and traceless and the factor 4 instead of 2 in the denominator is useful as one has to multiply further by 2 to account for two polarizations. The extra operator P( ) appears because of the original AID operators and reads as Noticec that a constant F W (¯ ) results in the second pole in the spin-2 Lagrangian and this is exactly the Weyl ghost observed by Stelle in [3,4]. Demanding that no new (and necessarily ghost) spin-2 excitations appear we must have either P(¯ ) = const > 0 or P(¯ ) = e −2ω(¯ ) where ω(¯ ) is an entire function in full analogy with γ(¯ ) in (4.15). The first choice results in a non-analytic F W (¯ ). The second choice on the one hand evades a ghost but on the other hand proves that only a truly non-local operator can generate a ghost-free spectrum. As the result presence of a non-constant ω(¯ ) is inevitable. The result for the power spectrum of tensor modes without slow-roll corrections as it was got in [43] is (4.19) Note that function ω(¯ ) must be evaluated at the position of the pole which isR/6 and also we restore M, the scale of non-locality in our model. Here we advance our study in comparison with [43] by computing also the tensor tilt and finding it up to the leading order in the slow-roll approximation. However, to achieve this, the next order slow-roll correction may need to be added to the tensor power spectrum (4.19), since in the case of the local R 2 model the leading term in P T does not depend on N = − ln k + const [7] (c.f. the recent paper [51] in this connection, too). A careful computation gives (4.20) We use this to obtain where prime is the derivative with respect to the argument and we have used that ǫ = 1/(2N ) with N the number of e-folds. Note that if ω ′ R 6M 2 = 0 we recover the tensor tilt of the Starobinsky model, i.e. n t = − 3 2N 2 which is a red tilt n t < 0.

r and modified consistency relation
Using a standard (local) results for the scalar power spectrum as advocated above (4.16) and modified tensor power spectrum (4.19) the tensor to scalar ratio is given by Therefore the presence of the Weyl tensor squared term in the action modifies the single field consistency relation (r = −8n t ) as follows (4.23) Using the fact that our computations do not depend technically on whether operators are of finite or infinite order in derivatives one can readily compare our results with an analogous derivation done in [52] where a pure Weyl tensor squared term was considered. In that case our answers can be shown to match. We however stress again that only a truly non-local operator, i.e. AID operator is necessary to get rid of the Weyl ghost.

Conclusions
R 2 always stood as the most successful theory of inflation and it is now the best fit for the most recent Planck data. In a recent study [43] this model realized in the context of non-local gravity that was shown to be UV complete in the sense of having no ghosts (unitarity) and being super-renormalizable (or finite). This was the significant theoretical development which embeds R 2 inflationary paradigm in a finite theory of quantum gravity. The present paper further extends this previous study with more rigorous analysis of the action, generalized solutions for EOM and derivation of inflationary parameters that can be tested in the future CMB data.
We have started by deriving explicitly that a most general theory of gravity that contributes to the linearly perturbed EOM around MSS contains the Einstein-Hilbert term, R 2 and Weyl tensor squared terms with AID operators in between, and the cosmological constant, see (2.4) and (2.5). Using the power of Bianchi identities we were able to reduce the final action presented in [45] effectively eliminating the Ricci tensor squared term. Our proof applies to any theory of gravity in any dimension. It is worth mentioning that SFT [17] provides a natural motivation for such kind of AID actions.
We have proceeded by presenting a rigorous mathematical proof that the trace equation of a local R 2 gravity without the cosmological term i.e., R = r 1 R is the only solution to EOM of our AID quadratic gravity also without the cosmological term as long as spatially flat FLRW metrics are considered. This means that even though we have complicated the local higher derivative theory of gravity [3] with AID operators, the space of background solutions remains the same satisfying R = r 1 R. In the situation when the Weyl tensor term is not included from the very beginning the claim remains true for any metric which is space-homogeneous in the synchronous reference frame, for example, anisotropic metrics, in particular Bianchi I type configurations, etc.
Further, we have derived the full perturbed second order action of AID quadratic gravity without the cosmological term around general backgrounds satisfying R = r 1 R and conditions (3.6), significantly boosting the previous studies where perturbations were only computed around MSS [23,42,43]. This full second order action is undoubtedly useful for further studies in the framework of AID gravities.
One of the crucial immediate task to be done in a forthcoming study is to extend these achievements to the models involving the cosmological term and to other types of metric, for example space-inhomogeneous or those describing spherically symmetric solution. This will allow one to attack in full, for example, the study of non-singular and ghost-free bounce configurations [23,46] which require the cosmological term to be present in the action. Also one can use the AID quadratic gravity framework to reconsider the problem of the curvature singularity which was proven to be generic in the case of Bianchi I metric in a local R 2 model [53].
Using the above described tools which are constructed in the present paper for the first time, we have come with the inflationary predictions of AID quadratic gravity model such as scalar spectral index, tensor to scalar ratio and tensor tilt consistently computing them in the leading order of the slow-roll approximation.
On this way we have proven that in the leading slow-roll approximation scalar perturbations in our model are equivalent to those in a local R 2 gravity. Our analysis thus makes it transparent that the scalar power spectrum remains the same as in a local R 2 model proving therefore that our model is to be the best fit with the present constraint n s = 0.968 ± 0.006.
The tensor power spectrum however gets modified exactly due to a differential operator in the Weyl tensor squared term introducing thereby a new parameter associated solely with this operator. As a result, the tensor to scalar ratio r gets a correction by a parameter that can give any value of r < 0.07 following (4.22). As an interesting but not a surprising consequence the computed tensor tilt deviates from a local R 2 model and thus the consistency relation gets modified as in (4.23). This resembles the results obtained in [52] with a huge difference that our model can avoid ghosts by promoting an operator in the Weyl tensor squared term to an AID operator. In particular, Eq. (4.18) and the discussion thereafter gives a very simple and clear explanation that the only way to defeat the Weyl ghost is to introduce a truly infinite derivative operator as long as one allows only analytic dependence on derivatives.
Our current AID quadratic gravity model modifies the tensor power spectrum and consequently r by a new parameter which is associated in this model with a quantum gravity prescription in the UV regime. This is in contrast to many other "Starobinsky"like models in the market [54][55][56][57] which modify only the scalar power spectrum. The tensor tilt in our model gets a new parameter related to the scale of non-locality M. The value of this new parameter can be fixed by future observations of primordial B-modes. Therefore, inflation in AID quadratic gravity meets all the current CMB constraints by PLANCK and is undoubtedly a very interesting and natural target for future CMB probes. We emphasize also that despite the fact that we can have any value of r < 0.07 in this model the energy scale of inflation remains the same as the latter is determined by scalar perturbations. This is a noteworthy feature of our model which is absent in the scalar field attractor models, i.e. so called α-attractor models [58].
Our results for AID quadratic gravity theory provide a foundation for studying not only inflation but also bounce, black holes, late time acceleration etc. in this framework. Given the theoretical progress we have achieved in the present paper future studying of reheating, non-gaussianities and other crucial questions are very important and timely. Intensifying the study of more inflationary parameters in combination with constraints from the observational camp would allow to narrow, for instance, the scale of non-locality and to start shaping the non-perturbative quantum gravity.
The last property in the latter line is sometimes called algebraic Bianchi identity.
The (differential) Bianchi identity is: We note down the following important second rank tensors Rg µν , (A.9) Traceless Ricci: All these tensors are symmetric. An important third rank tensor is the Cotton tensor: The fourth rank Weyl tensor is: The Weyl tensor has all the symmetry properties of the Riemann tensor (A.4) and it is absolutely traceless, i.e. W µ αµβ = 0 . Moreover it is invariant under the conformal rescaling, i.e.Ŵ µ αβγ = W µ αβγ forĝ µν = Ω 2 (x)g µν . (A. 13) This implies that the Weyl tensor is zero on conformally flat manifolds (i.e. when the metric can have the form ds 2 = a(x) 2 η µν dx µ dx ν where η µν is the Minkowski metric with the same signature). In fact, one should keep in mind that in D ≥ 4 vanishing Weyl tensor is a necessary and a sufficient condition for the space-time to be conformally flat.
Applying the Bianchi identity to the Weyl tensor one can find Squaring the last equality we get Next one can compute Going further one finds We see that in 4 dimensions the combination is transverse. We have used here the normalization of the Schouten tensor and the traceless property of the Weyl tensor. The latter combination is named Bach tensor. It is symmetric, traceless, transverse (for D = 4) and has the conformal weight −2 (i.e. scales as Ω(x) −2 upon scaling of the metric by Ω(x) 2 ). The spatially flat FLRW Universe metric is t is the cosmic time and a(t) is the scale factor. We intrinsically assume zero spatial curvature. The Hubble function is H =ȧ/a with dot denoting the derivative with respect to t. Equivalently we write the FLRW metric (A.19) as τ is the conformal time such that adτ = dt. Spatially flat FLRW Universe is conformally flat and the Weyl tensor in it is identically zero. The background quantities in the latter metric are We use the index "0" for the τ -component of any tensor in order not to confuse with just a small Greek letter (the cosmic time is used less often and wherever needed we designate it with the index t). Latin small letters from the middle of the alphabet are spatial indices, prime is the derivative with respect to the conformal time τ . Further: where s is the d'Alembertian operator acting on scalars.
In using non-local operators F X ( ) where X is some notational index we always assume that these operators are analytic functions of their arguments. This allows to write them in a Taylor series representation The metric perturbations are introduced as For any other quantity ϕ apart from the metric we useφ for its background value and δϕ for linear corrections. A spatial Fourier transform used to study perturbations is k is the spatial comoving momentum.

B Reduction of (2.3) to (2.4)
As the first step of this procedure we notice that one can use Schouten tensor S µν instead of L µν in (2.3). This is a linear transformation of L-tensor and it leads to some redefinition of functionF R ( ). Using that F L ( ) is an analytic function we can write it as a Taylor series and the corresponding expression is The zero term in this series can be dropped thanks to the presence of the Gauss-Bonnet (GB) invariant in 4 dimensions. With -s in between there is no such an obvious possibility. However, we can write the series without the zero term as follows We have moved one d'Alembertian to the left using an integration by parts since we are doing our computation under the integral.
The second step is to express S µν using Bianchi identity (A.16) and the right most S µν through L µν and g µν R. Schematically without explicit coefficients this can be written as Here x i , y i are numeric coefficients. We observe the following by primary inspection: terms proportional to ∼ W LL or ∼ LLL would not survive under the second variation around MSS as they contain 3 multipliers which are trivial on the background; since Weyl tensor and L-tensor are fully traceless any contraction of them with δ-symbol vanishes. Terms which still can contribute are written below: A careful examination reveals that the last term in the first line of the latter expression can survive the second variation only if each of two L-tensors is varied. As such this term with R taking its background valueR resembles the L-piece in action (2.3). We thus can tackle this term recursively looping back to (B.1) and considering a new expression Obviously, this is essentially the same expression as in (B.1) with new series coefficients. Iterating more and more (infinitely many) times we eliminate this term entirely. At each iteration the zero term with n = 0 can be dropped thanks to the use of the GB invariant. Further, first term on the second line in (B.3) can be effectively absorbed in RF R ( )R piece. The second term of the last line only contributes to a second variation when both L-tensor multipliers are varied and R takes its background value. Thus, it is again a local L 2 term and we drop it thanks to the use of the GB invariant.
So for the moment what we are left with is just first two terms from the first line of (B.3). We rewrite these two terms for convenience: The second term here is somewhat simpler to play with and we will show using its example that one can commute covariant derivatives without paying attention to extra contribution exactly like it happens in Minkowski background. For the zero term in the series we simply write (modulo the constant coefficient where we have used the Bianchi identity for the Einstein tensor (A.8). So this is just another contribution to RF R ( )R piece. Given we have one d'Alembertian operator we can write the following (again we omit the constant coefficient (B.7) The first term on the right can be treated in full analogy with (B.6) with one extra d'Alembertian inside. All other pieces can be in fact can be treated again as (B.6). Indeed, a non-zero second variation would only be present if R in the first factor ∂ ν R and L-tensor in the second factor are varied. As such, Riemann and Ricci tensors take their background values and all the terms but first in the second parenthesis can be absorbed in the zero series term. With the same strategy in mind, for any n the n-th term in the series can be redistributed into (n − 1)-th term and RF R ( )R. Finally this recursion absorbs everything of interest into RF R ( )R. The first term in (B.5) in fact produces indeed non-vanishing remaining contributions. For the first term in series expansion without d'Alembertian we have (we again omit constant coefficient x 1 f L1 y 1 below) Given we have at least one d'Alembertian, we will first replace L µν → S µν which is possible since the difference is proportional to the metric and the Weyl tensor is traceless, and second, we will use (A.16) to express S µν like in (B.2). This yields (modulo overall constant multiplier x 1 y 1 ) Terms ∼ W LL or ∼ LLL cannot contribute to the quadratic variation as they are cubic in quantities which are identically zero on the background. The term x 6 L µ ν R can only survive having R taking its background valueR. Otherwise, the second variation of this term in the action will be zero. As such, this term resembles the LHS of the equality with one d'Alembertian less. One can recursevely apply (A.16) as long as no d'Alembertians are left at all. Without d'Alembertian oeprators in between the structure like in (B.8) is reproduced.
The term x 2 ∇ ν ∂ µ R can be transformed as follows (B.10) The first transformation is possible thanks to the same logic as explained after (B.7). According to it we move derivatives absorbing all newly emerging Riemann tensor terms inside of terms with one d'Alembertian operator less. The second transform is an equality which uses the definition of the Bach tensor B µν (see (A.18)). This tensor is transverse and as such this contribution vanish. This is exactly like that in D = 4. In higher dimensions the divergence of the analog of Bach tensor gains non-vanishing contributions of the form ∼ W L and as such similar in their properties to the second term which is already present. This second term (and possible its complements in D > 4) produces ∼ W L∂R contribution which does not survive the second variation as all three factors must be varied to be nontrivial. The last unaccounted term in (B.9) which actually survives is Thus to the moment we have the following result. We started with (B.1) and ended up with two non-vanishing contributions (B.8) and (B.11) while also functionF R ( ) got redefined due to absorption of similar terms into it.
In order to simplify (B.8) we move one derivative by integration by parts to L-tensor and use (A.15) to write an equivalent expression (modulo a constant factor): Here we followed the logic explained after (B.7) according to which we can move derivatives like in Minkowski space. We also have replaced one Weyl tensor with the Riemann tensor as their difference vanishes upon contraction with another (traceless) Weyl tensor. We can now apply the Bianchi identity (A.7) and transform the last expression to Here we moved derivatives like in Minkowski space again and returned to the Weyl tensor as this transformation is identical. We see that the last term in the last expression is just the LHS in (B.12) with an opposite sign. Therefore, expression in (B.8) is equivalent modulo a constant factor to and can be absorbed inside ofF W ( ). In order to simplify (B.11) we have to repeat the same sequence of moves like for (B.8) just four times longer. The bottom line is that under the condition that R satisfies the equation (3.4) with r 1 > 0 and a spatially flat FLRW space-time is assumed. The simplest way to solve this problem is to reduce the differential equation (C.1), which is the second order with respect to the Hubble function H ≡ȧ/a (the scale factor a(t) itself does not enter due to invariance under rescaling of all spatial coordinates), to an expression containing R(t) only which should be an identity if R is not a constant. Note first that one such solution is R ≡ 0, and then no conditions on A, B and r 1 arise. Let us assume further that R is not identically equal to 0. Let B = 0. Theṅ where we take the minus sign forṘ corresponding to decrease of R with time for definiteness. If r 1 = 0 and R = 0,Ṙ = 0 too. By differentiating Eq. (C.1) with respect to time and dividing both its sides byṘ, one gets  4)), too, we notice that this case having zero measure in the space of initial conditions for Eq. (3.3) does not require special consideration further.
We finally note that one can extend the above proof to any metric which is spacehomogeneous in the synchronous frame.
D Slow-roll approximation using (4.1) The slow-roll parameter ǫ defined in (4.1) and we note it down again here for references: During inflation which is a nearly dS expansion we have ǫ ≪ 1 and use it as the small parameter for our approximations. It is useful to compute We require ǫ ′ ≪ Hǫ which is the condition to have inflation last long enough. The background scalar curvature and its conformal time derivatives becomē R ′′ ≈ −2HR(ǫ ′ + Hǫ) + 6H 2 ǫ 2R ≈ −2H 2 ǫR ≈ HR ′ .

(D.3)
Substituting these expressions into (3.4) one finds an important relation among our parameters As such, radiation contributes at least at the order ǫ 2 .

E Second order variation of action (3.1) around (3.4) and (3.6)
We will actively use where both quantities are analytic in thanks to condition F ′ R (r 1 ) = 0 (see (3.6)). We also remind few sums: This allows to compute (δF R ( ))R = Z 1 ( )(δ )R , The most tedious piece in varying (3.1) is the one with F R . Its variation contains 10 terms which explicitly can be written as where we use as well as the above derived relations.
To advance further we compute some quantities involving h µν δg µν = −h µν , δΓ ρ µν = γ ρ µν = and for the d'Alembertian acting on scalars We do not use the subscript "s" because no variations of other (acting on tensors for example) d'Alembertians are used in the paper. Using this latter expression we can show utilizing the integration by parts This in combination with the definition (4.6) simplifies (E.2) to is in fact the second order variation of the Einstein-Hilbert action. Accounting the Einstein-Hilbert and Weyl tensor terms in (3.1) and using conditions (3.6) we arrive to (4.5).
F Linearized EOM and proof that ζ = 0 in slow-roll approximation The variation of the trace equation (3.3) reads δE = ∂ µR ∂ µ + 2r 1R Z 2 ¯ k + 3F R ¯ k + R + 3r 1 Z 1 ¯ k ζ = 0 , (F.1) where ζ is defined in (4.6) and Z 1 , Z 2 are defined in Appendix E. The variation of the ( i j )-equation with i = j in the system (3.2) yields δE i j = −2λ k i k j a 2 F 1 (R + 3r 1 )(Φ − Ψ) + Υ + 2λc i j = 0 . (F.2) The variation of the ( 0 i )-equation in the system (3.2) yields The variation of the ( 0 0 )-equation in (3.2) yields Here we have used the notations Υ = F(¯ k ) − F 1 k − r 1 ζ + F 1 δR GI and Ξ = F(¯ k ) − F 1 (¯ k − r 1 ) 2 ζ , We see from the last equation that the single derivative term with ζ ν i in the trace equation can be recast in other ζ-s but also in the same ζ with an extra time-dependent factor. The orthogonality of Bessel and Neumann functions together with the appearance of this extra τ -factor proves ζ = 0 is the only solution to (F.1) in the leading slow-roll approximation. Going further one can work out the same technique using the quasi dS scale factor and Hubble parameter in terms of the conformal time (F.10) One can show upon construction of the ǫ-corrected d'Alembertian and corresponding eigenmodes that that an appearance of new linearly independent functions upon computing the derivatives is the blocking issue for any solution apart from ζ = 0.