Non-Gaussianities and tensor-to-scalar ratio in non-local R^2-like inflation

In this paper we will study $R^2$-like inflation in a non-local modification of gravity which contains quadratic in Ricci scalar and Weyl tensor terms with analytic infinite derivative form-factors in the action. It is known that the inflationary solution of the local $R+R^2$ gravity remains a particular exact solution in this model. It was shown earlier that the power spectrum of scalar perturbations generated during inflation in the non-local setup remains the same as in the local $R+R^2$ inflation, whereas the power spectrum of tensor perturbations gets modified due to the non-local Weyl tensor squared term. In the present paper we go beyond 2-point correlators and compute the non-Gaussian parameter $f_{NL}$ related to 3-point correlations generated during inflation, which we found to be different from those in the original local inflationary model and scenarios alike based on a local gravity. We evaluate non-local corrections to the scalar bi-spectrum which give non-zero contributions to squeezed, equilateral and orthogonal onfigurations. We show that $f_{NL}\sim O(1)$ with an arbitrary sign is achievable in this model based on the choice of form-factors and the scale of non-locality. We present the predictions for the tensor-to-scalar ratio, $r$, and the tensor tilt, $n_t$. In contrast to standard inflation in a local gravity, here the possibility $n_t$>0 is not excluded. Thus, future CMB data can probe non-local behaviour of gravity at high space-time curvatures.


Introduction
After the Planck mission, inflationary cosmology took a turn towards the Starobinsky inflationary model based on the modified R + R 2 gravity with one-loop corrections from matter quantum fields [1] (shortly dubbed as the (local) R 2 inflation afterwards), or to models with strongly non-minimally coupled scalar fields like the Higgs inflationary model [2], or to models emerged within string and supergravity (SUGRA) theories which have similar (very flat plateau-like) behaviour of their effective potentials in the Einstein frame [3]. These highlighted models comprise a set of most successful ones from the point of view of the observational data matching but neither of them features a manifestly well established UV completion setting.
In an attempt to tackle the issue of the UV completion for models of inflation on par with proposals to explain local higher order curvature terms like R 2 or non-minimal couplings of scalars by considering supergravity models [4,5] there are recent developments on studying R 2 inflation solution within an analytic infinite derivative (AID) gravity framework. This latter approach was shown to be absolutely consistent with existing cosmological observational data [6][7][8][9][10] while providing a viable track towards complete theory of quantum gravity [11][12][13][14][15][16][17][18][19]. We thus stick to this model in our effort to explore and analyze non-Gaussian perturbations and the tensor-to-scalar ratio during inflation in the present paper.
Analytic infinite derivative gravity essentially generalizes Einstein's General Relativity (GR) with quadratic in Ricci scalar, R and Weyl tensor, W terms in the action with analytic at zero functions of the d'Alembertian operator, also named form-factors, inserted in between like RF R ( )R where is the covariant d'Alembertian operator [14,20,21]. Analyticity of form-factors at zero implies a smooth low-energy limit of the model. This type of gravity models is also often abbreviated as IDG for infinite derivative gravity. We however would mainly stick to the AID abbreviation to stress the analytic behavior of the form-factors for low momenta. This is an essential mathematical point and this distinguishes these type of gravity modifications from those having non-analytic dependence on derivatives, like it is already encountered in one-loop effective quantum gravity action [22][23][24], for example, non-analytic forms like logarithms of the d'Alembertian or inverse of the d'Alembertian appear in the effective action as a result of integrating out matter fields [25,26].
Further it is the most general extension of GR around maximally symmetric spacetimes as it captures behavior of linear perturbations around such space-times in any analytic gravity generalization (see [27] for details). Inspiration for such a modification of gravity naturally comes from string field theory (SFT) where infinite derivative operators like e /M 2 s , with M s being a string scale M s M p enter the vertex terms of the string interaction [28][29][30][31][32][33][34]. M p here is the Planck mass as usual. Discussing gravity we name M s the scale of non-locality. Such a gravity is known to be ghost-free (unitary) improving thereby the proposal by Stelle [35,36] which suffered from the presence of the spin-2 ghost. Absence of ghosts requires the presence of an infinite tower of derivatives and thus forces one to consider essentially non-local form-factors [21]. This guaranties the powercounting renormalizability and raises the hope that Black-Hole and Big Bang singularities are being resolved in gravity theories of this kind. A number of crucial questions has been investigated recently regarding the resolution of singularities [37][38][39][40][41][42][43][44][45][46][47][48][49] in this approach. The arguments above explain why such an infinite derivative gravity tends to become a UV complete gravity theory.
Obviously, the higher derivative form-factors provide the crux of the gravity modifications which we study in this paper. We have freedom to choose form-factors such that the excitations spectrum is exactly like in a local R 2 gravity modification. Namely, there is a propagating massive scalar and a massless graviton in the theory. This scalar degree of freedom is identified with the 'scalaron' that drives inflationary expansion and is responsible for nearly scale invariant fluctuations which in turn seed the large scale structure formation. Although the structure of form-factors is constrained by the conditions on the theory to be ghost-free around maximally symmetric backgrounds [21,50], we lack for the time being a more fundamental approach of deriving them within the scope of string theory or SFT. Worth mentioning that recently some progress has been made to fix form-factors within asymptotic safety approach to quantum gravity [49,51]. However, still cosmological observations is the only major way to get a guidance for any gravity modification and thus inflationary cosmology and CMB observations play a vital role.
It was shown in [10] that an analytic non-local extension of the R 2 inflation produces up to the leading order in the slow-roll expansion exactly the same value of the spectral index, n s , of the Fourier power spectrum of primordial scalar perturbations generated during inflation, but the tensor to scalar power spectra ratio, r, as well as the tensor power spectrum consistency relation, can get modified. Moreover, the modification of the tensor power spectrum is solely due to the form-factor in the quadratic in Weyl tensor term in the action.
In the present paper, we aim to understand and constrain the analytic infinite derivative form-factors in the considered gravity model using the current constraints on inflationary parameters [52] such as (n s , r) and the non-Gaussianity parameter (f N L ) in the squeezed, equilateral and orthogonal configurations which read as n s = 0.9649 ± 0.0042 at 68%CL , r < 0.064 at 95% CL , (1.1) and with respect to Planck, BICEP2/Keck Array and BK15 with TT,TE, EE+lowE+lensing data [52,53]. Non-Gaussianity is an important tool to understand primordial Universe and especially the nature of scalaron or inflaton 1 , especially whether it is a canonical field, or a noncanonical one, or if multiple fields are responsible for the period of inflation, or if the inflation is driven effectively by a single field in the presence of relatively heavy fields, or if inflation happens in a higher order scalar tensor theories [54][55][56][57][58][59]. Every model produces a distinct signal for f N L for different triangular templates of 3-momenta. In the case of single canonical field models of inflation non-Gaussianities are very small as it was shown in [54]. Local R 2 model can be seen as a single canonical field model upon conformal transformation of the R + R 2 gravity action to the Einstein frame where a minimally coupled scalar field with a very flat plateau shape potential appears [60]. Any inflationary framework beyond a canonical scalar field may lead to different observable signatures of f N L [61,62]. In our case there is only one propagating scalar, namely the curvature perturbation, which is approximately time-independent on super-Hubble scales. However, the presence of higher derivative and even non-local operators can in principle lead to non-trivial signatures in the bi-spectrum (see for instance earlier observation of large non-Gaussianities in the context of SFT inspired p-adic inflation [63,64] where a minimally coupled non-local scalar field drives inflation).
The paper is organized as follows. In Section 2 we present the modified gravity model to be studied and discuss the generic structure of form-factors for the theory to be ghost free. We summarize the general conditions on form-factor F R ( ) for the theory to admit the local R 2 inflationary solution which satisfies the equation R = M 2 R where M is the scalaron mass. In Section 3 we discuss second order inflationary perturbations and present quantitative results for scalar and tensor power spectra and their tilts for sample formfactors. In Appendix A we provide equations of motion, slow-roll conditions and further technicalities on second order scalar perturbations. In Section 4 we compute the 3rd order variation of the studied gravity action around the local R 2 inflationary solution up to the leading order in the slow-roll approximation. Using the mode functions computed in Section 2 we calculate the 3-point correlation functions and the non-linearity parameter f N L as functions of momenta. Namely, we calculate the non-Gaussianity parameter f N L for local, equilateral and orthogonal configurations. To do quantitative analysis, we select sample form-factors and discuss the bounds on the scale of non-locality, M s , within the current CMB constraints. This section is supplemented by detailed computations in Appendices B and C. The Conclusion Section summarizes the performed analysis.

Analytic infinite derivative gravity and inflationary solution
In this section, we will review the equations of motion and inflationary solutions of AID gravity. The notations commonly used throughout the paper are as follows. The reduced Planck mass in our notation is where G is the Newtonian constant. The metric signature we work with is (−, +, +, +). The 4-dimensional indices are labeled by small Greek letters and 3-dimensional quantities are labeled by i, j = 1, 2, 3. Throughout the paper a bar over quantities is used to indicate their background values for Friedmann-Lemaître-Robertson-Walker (FLRW) metric g µν = (−1, a 2 , a 2 , a 2 ), where a = a(t) is the scale factor and t is the cosmic time. We use notation and a dot to represent the differentiation with respect to conformal time τ and cosmic time t respectively. ( †) and ( ‡) denote first and second derivatives with respect to the argument for various functions. (1,2,3) denotes in most cases the order of the performed variation. R, R µν , W µνρσ denote the Ricci scalar, Ricci tensor and Weyl tensor respectively. Throughout the paper, the sign " ≈ " denotes the leading order in the slow-roll approximation. Also we use * to indicate the quantities evaluated at the Hubble radius exit k ≈ aH, where k is the co-moving momentum and H =ȧ/a is the Hubble parameter.
The action we study contains non-local quadratic in curvature terms which was shown to be most general action around maximally symmetric space-times [9,10] Here s = /M 2 s with M s being the scale of non-locality and λ is a dimensionless parameter useful to control the effect of higher curvature contributions as the whole. The form-factors are analytic at zero and can be Taylor expanded as follows where f Rn , f W n are dimensionless constants. The vacuum equations of motion for the above action are given by (A.1) in Appendix A.
The form-factors should not be arbitrary as they are restricted by the ghost-free conditions. For instance, the generic structure of form-factors that leads to an extra scalar field of mass M and a massless tensor degree of freedom around maximally symmetric space-times backgrounds was studied in [14,21].
In the inflationary context we have to consider the structure of form-factors around the de Sitter background studied in [21] which we discuss in the next Sections. To study inflation, we start with FLRW metric which is conformally flat and as such the corresponding Weyl tensor vanishes. Therefore the trace of equations (A.1) for conformally flat backgrounds becomes It was shown in [10] that the local R 2 inflationary solution which satisfies the following equation¯ becomes the only (inflationary) solution 2 of AID gravity as long as form-factor F R ( s ) obeys the following conditions 9) where N ∼ 50 − 60 is number of e-foldings before the end of inflation. Note that the Hubble-slow-roll parameters , η used here in the original Jordan frame differ from the potential slow-roll parameters V , η V usually known in the context of scalar 2 The local R 2 inflationary solution satisfies (2.4) with the scale factor a ≈ a0(ts − t) − 1 6 e − M 2 (ts−t) 2

12
, H =ȧ a ≈ M 2 6 (ts − t) + 1 6(ts − t) , where ts 1 M denotes the end of inflation. 3 This can be seen heuristically by expanding the quadratic Ricci scalar part of the action (2.1) as In the above action we see that at high curvature regime R M 2 quadratic curvature terms are naturally dominant. However, it is easy to see that R 2 term is scale invariant (i.e., the term is invariant under local scale transformations gµν → e λ gµν ) and all the higher derivative terms are not scale invariant. So it is obvious to see that the hierarchy M 2 M 2 s is essential for the model to be compatible with cosmological data. field inflation (in the Einstein frame representation of the R + R 2 model) [60]. The second Hubble-slow-roll parameter in our case is η = + d ln 2dN ≈ 1 24N 2 (note that the calculation of its actual value requires using the second, next to leading, order of the slowroll approximation presented in Eq. (2.5)). On the other hand, in the Einstein frame V |η V |. Since all our study is in Jordan frame, in all the computations we apply the slow-roll approximation up to the leading order in .
From the considerations (2.8) and (2.9) we assume the hierarchy of mass scales in the theory as shown on the scheme Fig. 1. Further we refer to the condition (2.9) as the slow-roll approximation.

Second order perturbations and power spectra
Inflationary observables of R 2 -like inflation in AID gravity related to two point correlations were studied in detail in [8][9][10]. In this section, we will briefly review the second order perturbations of scalar and tensor parts and highlight the points essential for the rest of the paper without additional referencing. We start with the following perturbed line element in terms of gauge invariant Bardeen potentials (Φ, Ψ) and transverse and traceless tensor fluctuation h ij This is equivalent to fixing the Newtonian gauge. Study of perturbed linear equations of motion in the slow-roll approximation, i.e. in the quasi-de Sitter regime gives among other important results The solution of the above equation leads to Φ + Ψ ≈ 0 in the quasi-de Sitter limit in full analogy with the local R 2 inflationary model [23]. Proceeding with the construction of an action for a canonical scalar variable one must restrict the form-factor in order to avoid ghosts and this leads to the following form of the operator function F R ( s ): where γ S is an entire function of the d'Alembertian operator. Imposing the conditions (2.6) on the form-factor F R ( s ) in (3.3), we can deduce that We can also see that condition (3.4) being used in (3.3) does not give us any relation between scalaron mass and the non-locality scale. The second order action for the canonically normalized scalar has the form where the subscript (S) stands for the scalar part and u ≈ 2F 1R Ψ. A solution for Fourier modes ofũ = au upon the spatial Fourier transform yields in terms of Hankel functions [60] Adiabatic vacuum initial conditions for the observable range of Fourier modes of perturbations following from their quantization in the deep sub-Hubble limit |kτ | 1 are c 1k = 1, c 2k = 0. The Bunch-Davies initial conditions taken literally would correspond to imposing these conditions for all Fourier modes k including the limit k → 0. However, this is not possible for the most realistic inflationary models including the local R 2 model. For our purposes, it is sufficient to impose the adiabatic vacuum initial conditions for all Fourier modes which are deep inside the Hubble radius at the beginning of the inflationary stage. In the leading order in slow roll, the sub-Hubble limit of the mode function can be approximated asũ The curvature perturbation is defined as and is (approximately) time-independent on super-Hubble scales k aH. The primordial power spectrum and its tilt can be computed as where N is the number of e-folds and = −Ḣ H 2 ≈ 1 2N . We can notice that the scalar spectral index remains exactly same as for the local R 2 model. From the Planck data [66] the power spectrum at the pivot scale k * = 0.05Mpc −1 is P * R ∼ 2.2 × 10 −9 . Using this we get the mass of scalaron, Hubble parameter and Ricci scalar during slow-roll as The second order action for the canonically normalized tensor perturbations in the de Sitter approximation (applyingR M 2 ) is where the subscript (T ) stands for tensor part and γ T here is an entire function 4 which is related to the form-factor as Following the standard procedure for getting the tensor power spectrum and its tilt one yields The crucial difference here in comparison with local R 2 model is that the tensor power spectrum is scaled by an exponential factor of γ T evaluated at the pole of the tensor modē s =R 6M 2 s and accordingly tensor tilt also gets modified with a term proportional to the derivative of γ T at¯ s =R 6M 2 s . The ratio of tensor-to-scalar power spectra is given by (3.14) In the local R 2 model r = 12 N 2 = 3(1 − n s ) 2 as it follows from the original computation of scalar and tensor power spectra generated during inflation provided in [65]. We can clearly notice that if γ T −R 2M 2 s = 0 one exactly recovers the same prediction as in the local R 2 model. A deviation from the local R 2 model happens if From (3.14) and (3.13) we can see that the values of (r, n t ) explicitly depend on the ratio of scalar curvature to the square of M 2 s and the choice of entire function γ T ( s ). The scalar curvatureR during inflation can be read from (3.10). The two observables (r, n t ) can be fixed by the parameter space From the latest Planck 2018 data [52] we can deduce the following constraint (for N = 55) We have no constraint on the tensor tilt however. Imposing a heuristic limit on tensor tilt −4 n t 4 we get a constraint on parameter space (3.16) as To illustrate novel configurations which become accessible thanks to the presence of the AID operators let us consider the following form of the entire function where the approximation sign means that contributions of higher orders of the d'Alembertian operator are negligible when one evaluates this entire function at s =R 6M 2 s . From the point of view of the UV completion the entire function here should provide a suppression of the propagator at large momenta and as such one should worry about the sign of the coefficient in front of the highest degree of the d'Alembertian operator [67]. On the other hand this implies that given there are terms O( s ) in (3.19) one can freely choose signs of parameters α 1,2 .
In this particular example we see that r depends on α 2 and M s because the coefficient of α 1 vanishes when evaluated at s =R 6M 2 s . But the tensor tilt n t depends on (α 1 , α 2 , M s ). In figures Figs. 2 and 3 we present the predictions for (r, n t ) for some values of parameter space (α 1 , α 2 ) and we take M s H following the hierarchy in the scheme Fig. 1. Figure 2. In the above plots, we depict (r, n t ) versus the scale of non-locality M s for the formfactor (3.19) and N = 55. We have taken α 2 = 0.52, 0.28, 0.13, 0.065 (corresponding to brown full and dashed, black full and dashed lines respectively). In the right panel the lines with n t > 0 correspond to α 1 = 2.5 and the lines with n t < 0 correspond to α 1 = −2.5. In both the plots, we can notice that in the limit M s → M p we recover the predictions of the local R 2 model.
In the above plots, we depict (r, n t ) versus the scale of non-locality M s for the formfactor (3.19) and N = 55. We have taken α 2 = −0.06, −0.04, −0.025, −0.01 (corresponding to solid red, dashed red, solid blue and dashed blue lines respectively). In the right panel the lines with n t > 0 correspond to α 1 = 2.5 and the lines with n t < 0 correspond to α 1 = −2.5. In both the plots, we can notice that in the limit M s → M p we recover the predictions of the local R 2 model.
Given we have free parameter space related to the form-factor between the Weyl tensors in the action, the predictions of the non-local R 2 -like model for (r, n s , n t ) lies within the future of scope of CMB experiments. Detection of primordial gravitational waves in the view of future CMB experiments such as BICEP2/Keck, CMB S-4, SO, LiteBIRD and PICO can fix the form-factor and the scale of non-locality [68][69][70][71][72][73][74][75]. In figure Fig. 4 we compare the predictions of non-local R 2 -like model with the local R 2 model in (n s , r) plane. Even though r and n t are not detected so far, the latest Planck data presents a likelihood analysis for the values of n t which is an important parameter for inflationary setup as shown in figure Fig. 5. In the single field models of inflation we get r = −8n t (which holds for the local R 2 model as well) and this is in general proposed to be a crucial test. In multifield models or non-canonical models of inflation this consistency relation gets violated due to isocurvature perturbations and/or non-zero sound speed of curvature perturbation [76,77]. In many of the inflationary setups (except some special cases [78]) tensor tilt is found to be negative and it is regarded as a unique test for inflationary framework given that some alternative frameworks to the inflationary paradigm predict a blue tilt of tensor fluctuations [79].
As it becomes obvious in our case of a non-local R 2 -like inflation we can achieve the blue tilt. Therefore, we stress that a detection of a blue tensor tilt cannot rule out inflation in general.

Third order perturbations and inflationary bi-spectrum
In this section, we compute the 3rd variation of the action (2.1) around (2.4) within the slowroll approximation of the local R 2 inflationary solution. Then we compute the inflationary bi-spectrum of AID gravity and we especially quantify non-local corrections to the local R 2 gravity up to the leading order in the slow-roll approximation. Using the standard methods, the expectation value of 3-point correlations can be computed as [54,56,[80][81][82] where k i are wave vectors and H int is the interaction Hamiltonian. It is related to a perturbation of the Lagrangian (2.1) expanded up to the 3rd order in curvature perturbations (L 3 ) as H int ≈ −L 3 that is a valid approximation during slow-roll inflationary regime [54,56].
The integral (4.1) describes non-linear evolution of inflationary perturbations produced in the adiabatic vacuum initial state which we compute in the limit when cosmological scales exit the Hubble radius. Since during quasi-de Sitter expansion τ ≈ − 1 aH , it is a good approximation to calculate the integral in the limit τ e → 0 [56]. In the Fourier space, the three-point correlation function of curvature perturbations is defined as is called the bi-spectrum and |k i | = k i . Due to the translational invariance, the total momentum K = k 1 + k 2 + k 3 is conserved. Scale invariance and rotational invariance imply that the 3-point function B R (k 1 , k 2 , k 3 ) has to be homogeneous function of degree −6 that means and the number of independent variables of the 3-point function reduces to 2 which are the ratios k 2 k 1 , k 3 k 1 [83]. In the template of CMB observations, non-linear corrections in the curvature perturbation expressed as [84,85] where R g is the Gaussian random field and the f N L is the non-linearity parameter 5 . Then is called the amplitude of bi-spectrum. In the momentum space, f N L is a function of three wave numbers and it is related to the amplitude of bi-spectrum as The f N L defined above is often called as the reduced bi-spectrum. We calculate below the action (2.1) in the cubic order in curvature perturbations and the corresponding 3-point correlations (4.1) up to the leading order in the slow-roll approximation using the following procedure 6 . We presented details of computations and results in the Appendix C. 5 The factor of 3 5 comes from the relation between curvature perturbation R and the Bardeen variable Φ (Newtonian potential) during matter domination which reads as R = − 5 3 Φ, 6 Usually, in scalar field models of inflation, the 3rd order action is computed using the Arnowitte-Deser-Misner (ADM) formalism where spatial slicing (reparametrizing spatial coordinates) is chosen such that perturbations of a scalar field always vanish. In this gauge choice the spatial part of the metric becomes e 2R δijdx i dx j [80]. We however do computations in the Jordan frame and thus do not require such a gauge choice related to slicing of a 3-dimensional metric. We stick to a more convenient in our case gauge invariant approach which leads to the gauge invariant metric (3.1). This corresponds to choosing the Newtonian gauge.
• We use the fact that Φ + Ψ ≈ 0 during the inflationary period as it was shown in [9,10] and noted in the previous Section.
• From the second order action for canonical variable Υ, we can deduce that within the quasi-de Sitter approximation. In other words, we treat Ψ as an eigenmode of the d'Alembertian in the computation of the 3rd order action up to the leading order slow-roll approximation. To perform this step, first we must eliminate the terms in the 3rd order action which are proportional to linearized equations of motion via a suitable field redefinition [54].
• We carefully keep terms up to the leading order in and neglect higher order contributions treating ≈ const. In this regard, to convert the 3rd order action in Ψ into curvature perturbation R, we make substitution Ψ ≈ − R which follows from (3.8).
• Using the following approximations, we can bring the 3rd order action into a convenient form to compute 3-point correlations (4.1) where O ( s ) can be any analytic non-local operator which can be Taylor expanded in terms of d'Alembert operators. The second relation above is the result of the following general commutation relation [86] where φ is a scalar. In the slow-roll approximation, R µν ≈ R 4 g µν , so we can approximate it as (4.10) • We must be careful when applying the quasi-de Sitter approximation, especially when infinite derivatives are acting on a quantity. For example, the background Ricci scalar should not be treated as a constant until all the infinite tower of d'Alembertians are applied, otherwise we might end up with a zero contribution. As a consequence we would miss some terms when collecting next order non-zero contributions. For example, let us consider the following term Looking carefully at the above derivation, we can see that passing from the first line to the second one, our assumption of approximatingR ≈ const can be perfectly justified for every action of the d'Alembert operator, but using the same approximation in the next step would give us a null result due to the fact that 0 that follows from (2.6). Thus, we do not treatR = const in the second line, rather we apply the background solution¯ R = M 2R andṘ ≈ −2RH , and then we capture next order terms in the slow-roll approximation that provides us with new interactions of curvature perturbations and leads to a non-zero contribution to the 3-point function (c.f. (C.13)). The crucial lesson here is if any result of a computation gives zero, we must go to next to leading order in slow-roll, as far as we want to find a non-zero answer.
• In our approach, perturbed modes which began in the adiabatic vacuum state when being deep inside the Hubble radius during inflation are evaluated when they left this region (rigorously speaking, after some amount of e-folds, large but still much less than H 2 /|Ḣ|, when they had reached their constant value in the super-Hubble regime). Given the fact the leading mode of curvature perturbations remains constant after that as far as k aH, the three point correlations of R do not evolve there, too. This means that they keep information of primordial interactions of scalar modes during the period −∞ < τ < τ e where τ e ∼ − 1 K is a (conformal) time after few efoldings of the Hubble radius exit [87].
Using the above steps, we write below the result of a long and tedious computation presented in the Appendix C. Our obtained cubic order action in R of AID gravity in the leading order slow-roll approximation is where T i 's are dimensionless constants which can be read from (C.46). Here we present the final result for the amplitude of bi-spectrum after numerous simplifications and neglecting terms that are higher order in slow-roll: where the terms A i indicate the contributions from 8 types of interactions of curvature perturbations which can be read from (C.37) to (C.44) where K = k 1 + k 2 + k 3 and The non-Gaussianity parameter f N L can be obtained by substituting (4.13) in (4.6).
In the AID gravity, we can qualitatively see that we depart from the local theory due to the additional non-local interactions of the curvature perturbation listed from (C.40) to (C.45). This can for example imply that the long wavelength mode that is exited the Hubble radius can strongly interact with the short wavelength mode that is crossing the Hubble radius. Although this effect is present in the case of standard single field inflation [54,88], it is significantly suppressed by slow-roll. Enhanced non-Gaussianities beyond the standard single field only known to occur in the context of non-canonical and/or multifield models of inflation [88,89]. Our case significantly differs from this so far well known physics of inflation. We have obtained enhanced non-Gaussianities due to non-localities present in AID gravity leading to non-local interactions of the only propagating scalar mode of the theory 7 . A trivial observation we can make is that the contributions (C.40) to (C.45) do not appear in local theory. From the observational point of view, three configurations of reduced bi-spectrum f N L for squeezed (k 1 → 0, k 2 = k 3 = K 2 ), equilateral (k 1 = k 2 = k 3 = K) and orthogonal (k 1 = k 2 = K/4, k 3 = K/2) are the most relevant.
We can easily verify that our computation of bi-spectrum reduces to the result of the local R 2 model in the limit T * NL → 0 and F ( †) R R 4M 2 s → 0. Especially, we can verify here the result of squeezed shape (k 1 which is exactly the Maldacena consistency relation found in the single scalar field inflation [54]. The shape of reduced bi-spectrum in non-local R 2 can be distinguishable from that in the local R 2 theory [83] due to the non-local corrections. Considering the general structure of the form-factor in (3.3), we obtain three shapes of reduced bi-spectrum f N L in the leading order approximation as (4.17) In the above result of the reduced bi-spectrum, we have presented only the leading order terms in from both the local and non-local contributions (c.f., (4.13)). The parameter space that determines the different shapes of reduced bipsectrum (4.17) is we recover the result of local R 2 theory in the leading order [80].
To see quantitatively the measure of non-Gaussianities, let us consider the following couple of choices of entire function γ S which are compatible with (3.4). Moreover we would require an entire function to be compatible with (A.13) which reflects an assumption that the higher derivative form-factors do not break the slow-roll approximation scheme (see Appendix A.2 for a more detailed explanation).
Our first simplest choice of entire function is as below The second choice we consider below is higher order in d'Alembertian with a property of γ S R 4M 2 s = 0. This implies the second terms in (4.17) vanish but the 3rd terms can be present and dominate the local contribution depending on the scale of non-locality (4.20) As in the previous Section the approximation sign means that there are higher degree terms in the d'Alembertian operator which would guarantee the supression of the propagator at large momenta but on the other hand provide a freedom choosing signs of arameters β 1,2 .
For the above form-factors (4.19) and (4.20) we plot below in figure Fig. 6 f N L as a function of M s for squeezed, equilateral and orthogonal configurations. Plots are made for N = 55 and β 1 = 1, β 2 = −2. We can read that for M s 5.6 × 10 −5 M p , we are well within the bounds of (1.2). If the scale is 5.6 × 10 −5 M p < M s < 6.5 × 10 −5 M p , we get f N L ∼ ±O (1). We can also notice that all f N L values for different configurations are negative for the case of (4.19) and positive in the case of (4.20). The crucial and significant result here is that the R 2 -like inflation in AID gravity can give f N L ∼ O(1). Moreover, the non-Gaussianity parameter can be of either sign on contrary to the local R 2 model which yields only positive values. The large-scale structure observations provide a test-bed for large values of non-gaussianities [90,91]. Detecting non-Gaussianities will fix up to some extent the scale of non-locality and the form-factor. More importantly, we obtain f sq N L ∼ O(1) in our case which has been only known to occur in multifield models of inflation where isocurvature perturbations appear. In our case, although we have only a curvature perturbation here which gets frozen on super-Hubble scales, corrections due to the presence of non-locality can result in f sq N L ∼ O(1). Similarly, detectable levels of f eq N L , f ortho N L have been so far known to be possible in non-canonical models of inflation where the sound speed of curvature perturbations is much less unity. In our case, even though the sound speed of curvature perturbations is unity during inflation, non-local interactions give us f eq N L , f ortho N L ∼ O(1). This makes the R 2 -like inflation in AID gravity an interesting target with respect to future CMB data.

Conclusions
In this paper we have studied an embedding of the local R 2 inflationary solution in the framework of a quadratic in curvatures non-local gravity with analytic infinite derivative form-factors. The essence of the present paper is the computation of the scalar non-Gaussinaities (3-point scalar correlations) and the subsequent attempt to constrain the higher-derivative form-factors based on the existing observational data.
A strong point of R 2 -like inflation in the presented non-local gravity setting is that the scalar spectral index n s remains the same as in the local R 2 gravity and thus best fits the current CMB data [52]. However, the tensor-to-scalar ratio and the tensor tilt get modified from the standard consistency relation due to the addition of the Weyl tensor squared term with an analytic non-local form-factor in the action [10]. The analysis reveals that the tensor-to-scalar ratio can acquire any value from the present upper bound all the way down to 10 −11 . Furthermore, the tensor tilt is not sign constrained in this model a priori and can take any value in the range −O(1) < n t < O(1). In particular this disregards possible claims that a sign of the tensor tilt may prove or disprove the inflation as such and moreover a possible blue tensor tilt would highly favor the model outlined in the present paper. In general the performed analysis makes the presented scenario to be a very interesting candidate to be tested within the scope of several future experiments such as BICEP2/Keck, CMB S-4, SO, LiteBIRD and PICO [68][69][70][71][72][73][74][75].
As the main advance in this paper we have computed following the standard methods [54,82] in full scalar non-Gaussianities by means of deriving the 3rd order variation of the action around the R 2 -like inflationary solution up to the leading order in the slow-roll approximation. We demonstrate that the scalar bi-spectrum contribution comes only from the quadratic in Ricci scalar term. Assuming suitable form-factors which are compatible with ghost-free conditions of the theory, we have shown that this model is consistent with the current CMB data and leads to predictions that can be tested with future CMB experiments. Namely, our results show that it is possible to obtain f sq N L ∼ O (1) and moreover a sign of the non-Gaussianity parameter is not fixed. Usually any detection of such a large and positive non-Gaussianity is attributed to the presence of multiple scalar fields or to an effect of heavy fields giving rise to non-trivial evolution of curvature perturbations [62,89], or to non-slow-roll inflationary scenarios [57]. Our findings present a counterexample to this common point of view since we obtain just a single adiabatic mode of scalar perturbations which can lead to detectable values of non-Gaussianity in the squeezed (local) limit due to the presence of non-local interactions already at the tree level. Furthermore, building a chance to gain large non-Gaussianity parameter of either sign seems to be a unique feature of the present model.
Moreover, R 2 -like inflation in our setup would also lead to detectable non-Gaussianities in equilateral and orthogonal shapes which is again due to non-local interactions of curvature perturbations when its Fourier modes exit the Hubble radius during inflation. In turn, this presents an alternative to the occurrence of similar large non-Gaussianities in non-canonical models of inflation [55] in which the sound speed of curvature perturbations is much less than unity. We emphasize that the sound speed of curvature perturbations in our case is unity during inflation, and different non-Gaussian shapes are solely due to the effect of non-local interactions.
Therefore, any detection of f N L ∼ O (1) or any departure from a standard single field consistency relation can be easily attributed to higher derivative effects in the early Universe. A possible negative values of would be detected non-gaussianity parameters will raise even more the chance for our analytic infinite derivative gravity modification setup to become the favorite one. In particular such effects can be put to test in the large-scale structure observations [90][91][92]. Also our results suggest that the shape of the form-factor F R ( s ) can be probed by observations.
In all the instances detectable deviations from so called standard expectations for parameters are subject to the ratio of M s /M p which is engaged in the hierarchy with the scalaron mass as expressed in (2.8) and (2.9) For M s approaching the Planck scale, our results reduce to the ones known from the local R 2 model providing no novel detectable signatures of the model. However, the smaller is the scale of the higher derivative gravity modification, the greater is a possible departure towards large non-Gaussianities in particular. To get a more detailed picture, one needs to constrain the shapes of form-factors which is a well more complicated task requiring more theoretical thinking rather than comparison with observational data. Nevertheless, computing other types of 3-point correlations with the aim to get even more constraints from existing datasets is a very interesting question for forthcoming publications.

A Equations of motion, slow-roll approximation and perturbations
The equations of motion for the action (2.1) are [15] (A.2) In the above expressions ∇ and ; denote covariant derivatives.

A.1 Slow-roll approximation
During the quasi-de Sitter expansion,R is nearly constant that means In terms of the conformal time, the de Sitter space-time is defined by Quasi-de Sitter space is a slight departure from the exact de Sitter which can be defined by the slow-roll parameter where H =ȧ a is the Hubble parameter, M is the scalaron mass and N = ln a f − ln a is the number of e-folds during inflation counted from its end (a = a f ) backwards in time. The slow-roll parameter above is for the local R 2 inflation (2.5). The following relations are useful in many of our computations where H =ȧ a , H = a a are the Hubble parameter in the cosmic and conformal time respec-tivelyR We note that slow-roll would also mean M 2 = 6H 2 a 2 = 1 2R → 0 which can be easily deduced using (A.8) and (2.4) in the limit A.2 On Φ + Ψ ≈ 0 during quasi-de Sitter expansion One of the prime result of [10] is that second order action of scalar perturbations in AID gravity remains the same as the one in the local R 2 gravity, around the inflationary solution (2.5) that satisfies (2.4), up to the leading order in slow-roll approximation. Here we discuss this point briefly and heuristically show that non-local corrections at the second order level are indeed negligible. The exact from of second order action of AID gravity around FLRW space times that satisfy the background equation (2.4) reads as [10] where δ local is the second order variation of an action Moreover it was shown by solving perturbed trace equation of AID gravity that ζ = 0 in the leading order in the slow-roll approximation, and this is exactly equivalent to Φ + Ψ ≈ 0.
Since the first order variations of the Weyl tensor are proportional to Φ + Ψ ≈ 0, there will be no scalar contributions from the Weyl tensor term. As a result, we get the scalar power spectra like in the local R 2 model up to the leading order corrections as we discussed in Section 3 and we can also easily get it from (A.9). Now, we can heuristically translate this result into a condition on the form-factor F R ¯ s by computing the contributions of the second term in (A.9) in the next to leading order in the slow-roll approximation and showing them to be negligible compared to the contribution of the local term (A.10). To show this, we consider Φ+Ψ ≈ 0 and¯ Ψ ≈ M 2 Ψ. Now calculating slow-roll corrections to the second term in (A.9), we obtain where we applied the approximationR ≈ const, and F is the second derivative of the form-factor evaluated at the value of scalaron mass square. To declare that the contribution (A.12) is negligible compared to the local term from the quadratic part in the local action (A.10), we require Considering the hierarchy of scales as (2.8), we can easily satisfy the above condition if the second derivative of the form-factor evaluated at M 2 satisfy F

B Useful relations for the computation of bi-spectrum
In this section, we provide computations of perturbed curvature quantities and action of infinite derivatives on them within the slow-roll approximations. Many of the following simplifications were used in computing the 3-point correlations in the previous section. Perturbation functions are not space-homogeneous. It is convenient to perform a spatial Fourier transform which for some function ϕ(τ, x) in the flat FLRW metric reads ϕ(τ, x) = ϕ(τ, k)e i k x d k . (B.1) Using the above Fourier representation, we can writē 2) The first and second variations of the d'Alembertian operator are Below we provide the useful list of perturbations of scalar curvature up to the 3rd order Below we derive some recurrent relations when d'Alembertian operators act on perturbed quantities (C.1) Here we performed integration by parts and used the background solution (2.4). The contribution from the Einstein-Hilbert term is negligible because the R 2 term dominates during inflation.
Since Φ + Ψ ≈ 0 during inflation, we compute (C.1) up to the cubic order in Ψ ≈ R in the leading order slow-roll approximation. From [9], we know that the first variation of Weyl term only contributes to tensor perturbations and does not lead to scalar perturbations because Φ + Ψ ≈ 0 during inflation. Notice that all the terms in (C.1) contains at least one first order variation of the Weyl tensor. Therefore, these contributions can be approximated by zero for the scalar 3-point correlation. As a result, the cubic order scalar contributions can only arise from the quadratic Ricci scalar term in the action.
Moreover, for the computation of 3-point correlations (4.1) we use the mode functions evaluated from the second order action of curvature perturbations evaluated is the de Sitter approximation (3.5) which we recall here again expressing in terms of curvature perturbation R as The above action in the quasi-de Sitter approximation gives two-point correlations possessing a nearly scale invariant power spectrum which is Gaussian. Non-Gaussianities arise from 3-point correlations which are computed from the 3rd order action in the next to leading order in the slow-roll approximation. In the computation of the 3-order action it is important to eliminate terms proportional to the equation of motion for R (it can be obtained by varying the action (C.2)) via a suitable field redefinition of R. This is step is crucial to extract the leading non-linear contributions to curvature perturbations [54,80].
In the local R 2 gravity, it is known that non-Gaussianities are small [54]. However, in the non-local gravity we have non-local contributions to the bi-spectrum which can be read from (C.1). Note that in the calculation of bi-spectrum, we can perform the computation using the leading order in the slow-roll approximation.
In ( Below, we perform the term by term computation of the 3rd order action (C.1) in terms of the curvature perturbation R and the evaluation of 3-point correlation functions Note that in all the following calculations, we only write the interaction terms which are leading order in slow-roll. Also in the calculations we encounter terms proportional to the quantity which is found to be negligible (for the form-factors of the form (3.3)) compared to remaining terms which contribute significantly to the bi-spectrum. Therefore, we drop all these terms for economic reasons. Also we do not write the interaction terms which only lead to imaginary values of 3-point correlations, as by definition the 3-point correlation function is a real quantity (4.1). First we compute the local contribution of the 3rd order action (C.1) which is the 3rd variation of the local action below where ∂L 2 ∂R is the term proportional to the equation of motion of R following from (C.2). This term can be eliminated by a field redefinition of R → R + 4 R 2 which leads to a modification of the Gaussian term as This term leads to local contributions to the bi-spectrum which are (C.37) to (C.39) Now we calculate 3-point correlations from the following non-local term which does not involve variation of the form-factor This term contributes to (C.40) Let us consider the following term of (C.1) which involve first variation of the formfactor F R ( s ) (C.11) Passing from the first line to the second one, we have carried integration by parts and used the background solution¯ R = M 2R . Passing from the second line to the third one, we have substituted δ 0 from (B.4) and the result derived in (B.10). This term gives the following contribution to the bi-spectrum (C.40)-(C.42) Let us now consider the following non-local contribution of (C.1) (C.13) The contributions form this term are (C.40), (C. 44) and (C.45) (C.14) Let us now consider the following two non-local terms which require second variation of the form-factor F R ( s ). They can be written as two kinds of variations in the following way δR .

(C.15)
Calculating the first part in (C.15), we obtain (C.16) The contributions from this term are (C.40), (C. 44) and (C.45) Let us compute the second term The leading contributions from this term are (C. 41) and (C.42) Similarly, the following non-local contribution can be written as two parts (C.21) The first term in (C.21) can be simplified as The bi-spectrum contribution from this term is (C.38) Now the second term in (C.21) can be worked out as The terms involving 3rd variation of the form-factor can be split into 4 different type of terms as below We neglect the first interaction term RRR in (C.28) in comparison with the first term in (C.24) within the slow-roll approximation.
Let us now consider the third term in (C.26) and evaluate it as (C.30) The leading contribution from this term is (C.38) Now we collect all the terms above and calculate the bi-spectrum contribution from each type of interaction substituting the curvature perturbation evaluated in the adiabatic vacuum initial state (3.8) in (4.1). To illustrate, in the rest of the computations here, we calculate the 3-point correlation contribution from the interaction R∇R·∇R of (C.6) using (3.8) and (4.1) −i2 5 π 7 δ 3 (k 1 + k 2 + k 3 ) 1 1. Calculating the bi-spectrum contribution due to the interaction R∇R · ∇R, we get (C.37) 2. Calculating the bi-spectrum contribution due to the interaction RR 2 , we get 3. Calculating the bi-spectrum contribution due to interaction R 3 , we obtain 2 5 π 7 δ 3 (k 1 + k 2 + k 3 ) 1 i k 3 i P * 2 R T * 3 − K 3 3 + 2Kk 1 k 2 + k 1 k 2 k 3 3 + perms . (C. 39) 4. Calculating the bi-spectrum contribution due to the interaction RRR , we obtain 5. Calculation of the bi-spectrum contribution due to the interaction ∇R · ∇RR gives 2 5 π 7 δ 3 (k 1 + k 2 + k 3 ) 1 6. Calculating the bi-spectrum contribution due to the interaction R 3 , we obtain 7. Calculating the bi-spectrum contribution due to the interaction R ∇R · ∇R , we get 8. Calculating the bi-spectrum contribution due to the interaction RR R , we get 9. Calculating the bi-spectrum contribution due to the interaction R∇R · ∇R , we get 2 5 π 7 δ 3 (k 1 + k 2 + k 3 ) 1 i k 3 i P * 2 R T * 9 (k 1 · k 2 ) k 2 Here Ei(z) is the integral exponential function. The integral has to be evaluated in the limit τ e → 0 while convergence at large τ is made possible by the oscillatory behaviour at τ → −∞. Some of the integrals do diverge in the limit τ e → 0 where we must follow the guidance well explained in [87] which prescribes us to evaluate integrals up to the conformal times corresponding to the Hubble radius crossing scale, such that K ∼ aH ∼ − O(1) τe . For the purposes of practical computations one fixes Kτ e = −1.