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

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, 11] while providing a viable track towards completing theory of quantum gravity [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. We therefore 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 $R\mathcal{F}_{R}(\square)R$ where $\square$ is the covariant d’Alembertian operator, see: [6, 23, 15, 24, 16, 25, 26, 27]. 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 it distinguishes this type of gravity modifications from those having non-analytic dependence on derivatives, like it is already encountered in one-loop effective quantum gravity action [28, 29, 30, 31], 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 [32, 33].

Furthermore, it is the most general extension of GR around maximally symmetric space-times as it captures behavior of linear perturbations around such space-times in any analytic gravity generalization (see [26] for details). Inspiration for such a modification of gravity naturally comes from string field theory (SFT) where infinite derivative operators like $e^{{\square}/{\mathcal{M}_{s}^{2}}}$, with $\mathcal{M}_{s}$ being a string scale $\mathcal{M}_{s}\lesssim M_{p}$ enter the vertex terms of the string interaction [34, 35, 36, 37, 38, 39, 40]. $M_{p}$ here is the Planck mass as usual. Discussing gravity we name $\mathcal{M}_{s}$ as the scale of non-locality. Such a gravity is known to be ghost-free (unitary) improving thereby the proposal by Stelle [41, 42] 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 [27]. This guaranties the power-counting 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 [23, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55] 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 [27, 56], 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 [54, 57]. However, still observational cosmology 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 [11] 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 [58] such as $\left(n_{s},\,r\right)$ and the non-Gaussianity parameter ($f_{NL}$) in the squeezed, equilateral and orthogonal configurations which read as

$\displaystyle n_{s}=0.9649\pm 0.0042\,\,\textrm{at}\,\,68\%\textrm{CL}\,,\quad r<0.064\,\,\textrm{at}\,\,95\%\,\textrm{CL}\,,$

(1.1)

and

$$f_{NL}^{sq}=-0.9\pm 5.1\,\quad f_{NL}^{equiv}=-26\pm 47\,\quad f_{NL}^{ortho}=-38\pm 24\,,\,\textrm{at}\,\,\,\,68\%\,\,\textrm{CL}.$$

(1.2)

with respect to Planck, BICEP2/Keck Array and BK15 with TT,TE, EE+lowE+lensing data [58, 59].

Non-Gaussianity is an important tool to understand primordial Universe and especially the nature of scalaron or inflaton¹¹1Which is a scalar degree of freedom coming from addition of hypothetical matter to GR., especially whether it is a canonical field, or a non-canonical 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 [60, 61, 62, 63, 64, 65, 66, 67]. Every model produces a distinct signal for $f_{NL}$ 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 [60]. 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 [68]. Any inflationary framework beyond a canonical scalar field may lead to different observable signatures of $f_{NL}$ [69, 70]. In our case there is only one canonical 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 [71, 72] 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 $\mathcal{F}_{R}\left(\square\right)$ for the theory to admit the local $R^{2}$ inflationary solution which satisfies the equation $\square 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 form-factors. 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_{NL}$ as functions of momenta. Namely, we calculate the non-Gaussianity parameter $f_{NL}$ 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, $\mathcal{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.

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 $\frac{M_{P}^{2}}{2}=\frac{1}{16\pi G}$ 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_{\mu\nu}=(-1,a^{2},a^{2},a^{2})$, where $a=a(t)$ is the scale factor and $t$ is the cosmic time, like $\bar{R}$ for the background value of the Ricci scalar, $\bar{\square}$ for the background value of the d’Alembertian operator, etc. We use notation ^′ and a dot to represent the differentiation with respect to conformal time $\tau$ 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_{\mu\nu},\,W_{\mu\nu\rho\sigma}$ denote the Ricci scalar, Ricci tensor and Weyl tensor respectively. Throughout the paper, the sign “$\,\approx\,$” denotes the leading order in the slow-roll approximation.

The action we study contains non-local quadratic in curvature terms which was shown to be most general action around maximally symmetric space-times [15, 27, 26, 10, 11]

$$S=\int d^{4}x\sqrt{-g}\left(\frac{M_{p}^{2}}{2}R+\frac{\lambda}{2}\bigg{[}R\mathcal{F}_{R}\left(\square_{s}\right)R+W_{\mu\nu\rho\sigma}\mathcal{F}_{W}\left(\square_{s}\right)W^{\mu\nu\rho\sigma}\bigg{]}\right)\,.$$

(2.1)

Here $\square_{s}={\square}/{\mathcal{M}_{s}^{2}}$ with $\mathcal{M}_{s}$ being the scale of non-locality and $\lambda$ 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

$$\mathcal{F}_{R}\left(\square_{s}\right)=\sum_{n=0}^{\infty}f_{Rn}\square_{s}^{n},\quad\mathcal{F}_{W}\left(\square_{s}\right)=\sum_{n=0}^{\infty}f_{Wn}\square_{s}^{n}\,.$$

(2.2)

where $f_{Rn},f_{Wn}$ 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 [15, 26, 27].

In the inflationary context we have to consider the structure of form-factors around the de Sitter background studied in [27, 26] 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

$$E=(M_{P}^{2}-6\lambda\square\mathcal{F}_{R}(\square_{s}))R-\lambda(\mathcal{K}^{\mu}_{\mu}+2\tilde{\mathcal{K}})=0\,,$$

(2.3)

where $\mathcal{K}^{\mu}_{\mu},\,\tilde{\mathcal{K}}$ are infinite derivative square in scalar curvature terms defined in (A.2). In the sequel we put parameter $\lambda$ equal 1 for simplicity. It was shown in [11] that the local $R^{2}$ inflationary solution which satisfies the following equation

$$\bar{\square}\bar{R}=M^{2}\bar{R}\,,\,$$

(2.4)

becomes the only (inflationary) solution²²2The local $R^{2}$ inflationary solution satisfies (2.4) with the scale factor $\displaystyle a$ $\displaystyle\approx a_{0}(t_{s}-t)^{-\frac{1}{6}}e^{-\frac{M^{2}(t_{s}-t)^{2}}{12}}\,,\quad H=\frac{\dot{a}}{a}\approx\frac{M^{2}}{6}(t_{s}-t)+\frac{1}{6(t_{s}-t)}\,,$ (2.5) $\displaystyle\bar{R}=12H^{2}+6\dot{H}$ $\displaystyle\approx\frac{M^{4}(t_{s}-t)^{2}}{3}-\frac{M^{2}}{3}+{O}\left((t_{s}-t)^{-2}\right)\,,$ where $t_{s}\gg\frac{1}{M}$ denotes the end of inflation. of AID gravity as long as form-factor $\mathcal{F}_{R}(\square_{s})$ obeys the following conditions

$$\mathcal{F}_{R}^{\left(\dagger\right)}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)=0\quad,\quad\frac{M_{p}^{2}}{2}=3\lambda M^{2}\mathcal{F}_{1}\,,$$

(2.6)

where $\mathcal{F}_{R}^{\left(\dagger\right)}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)=\frac{\partial\mathcal{F}_{R}}{\partial\square_{s}}\Big{|}_{\square_{s}=\frac{M^{2}}{\mathcal{M}_{s}^{2}}}$ and $\mathcal{F}_{1}=\mathcal{F}_{R}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)$. Since cosmological observations indicate the nearly scale invariant power spectrum of scalar fluctuations, we consider the following hierarchy of scales³³3This can be seen heuristically by expanding the quadratic Ricci scalar part of the action (2.1) as $$S=\int d^{4}x\sqrt{-g}\left[\frac{M_{p}^{2}}{2}R+\frac{M_{p}^{2}}{12M^{2}}R^{2}+O\left(\frac{M_{p}^{2}R\square R}{M^{2}\mathcal{M}_{s}^{2}}\right)\right]\,.$$ (2.7) In the above action we see that at high curvature regime $R\gg 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_{\mu\nu}\to e^{\lambda}g_{\mu\nu}$ ) and all the higher derivative terms are not scale invariant. So it is obvious to see that the hierarchy $M^{2}\ll\mathcal{M}_{s}^{2}$ is essential for the model to be compatible with cosmological data. [10, 11]:

$$M^{2}\ll\mathcal{M}_{s}^{2}\lesssim M_{p}^{2}\,,$$

(2.8)

During inflation the following slow-roll condition for the background solution (2.5) is satisfied [73, 68]:

$$\epsilon=-\frac{\dot{H}}{H^{2}}\,\approx\frac{M^{2}}{6H^{2}}\approx\frac{1}{2N}\ll 1\implies M^{2}\ll H^{2}\,,$$

(2.9)

where $N\sim 50-60$ is number of e-foldings before the end of inflation.

Note that the Hubble-slow-roll parameters $\epsilon,\,\eta$ used here in the original Jordan frame differ from the potential slow-roll parameters $\epsilon_{V},\,\eta_{V}$ usually known in the context of scalar field inflation (in the Einstein frame representation of the $R+R^{2}$ model) [68]. The second Hubble-slow-roll parameter in our case is $\eta=\epsilon+\frac{d\ln\epsilon}{2dN}\approx\frac{1}{24N^{2}}\ll\epsilon$ (note that the calculation of its actual value requires using the second, next to leading, order of the slow-roll approximation presented in Eq. (2.5)). On the other hand, in the Einstein frame $\epsilon_{V}\ll|\eta_{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 $\epsilon$.

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.

Inflationary observables of $R^{2}$-like inflation in AID gravity related to two point correlations were studied in detail in [9, 10, 11]. 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 $\left(\Phi,\,\Psi\right)$ and transverse and traceless tensor fluctuation $h_{ij}$

$$ds^{2}=a^{2}\left(\tau\right)\left[-\left(1+2\Phi\right)d\tau^{2}+\left(\left(1-2\Psi\right)\delta_{ij}+2h_{ij}\right)dx^{i}dx^{j}\right]\,.$$

(3.1)

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

$$\left[\mathcal{F}_{1}\left(\bar{R}+3r_{1}\right)+\left(\bar{\square}-2H^{2}\right)\mathcal{F}_{W}\left(\bar{\square}+6H^{2}\right)\right]\frac{\Phi+\Psi}{a^{2}}=0\,.$$

(3.2)

The solution of the above equation leads to $\Phi+\Psi\approx 0$ in the quasi-de Sitter limit in full analogy with the local $R^{2}$ inflationary model [29]. 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 $\mathcal{F}_{R}(\square_{s})$:

$$\mathcal{F}_{R}\left(\square_{s}\right)=\mathcal{F}_{1}\frac{3e^{\gamma_{S}\left(\square_{s}\right)}\left(\square_{s}-\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)+\left(\frac{\bar{R}}{\mathcal{M}_{s}^{2}}+3\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)}{3\square_{s}+\frac{\bar{R}}{\mathcal{M}_{s}^{2}}}\,,$$

(3.3)

where $\gamma_{S}$ is an entire function of the d’Alembertian operator. Imposing the conditions (2.6) on the form-factor $\mathcal{F}_{R}\left(\square_{s}\right)$ in (3.3), we can deduce that

$$\gamma_{S}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)=0\,.$$

(3.4)

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

$$\delta^{(2)}S_{(S)}=\frac{1}{2}\int d\tau d^{3}x\sqrt{-\bar{g}}u\left(\bar{\square}-M^{2}\right)u\,,$$

(3.5)

where the subscript $(S)$ stands for the scalar part and $u\approx 2\mathcal{F}_{1}\bar{R}\Psi$. A solution for Fourier modes of $\tilde{u}=au$ upon the spatial Fourier transform yields in terms of Hankel functions [68]

$$\tilde{u}_{k}=\frac{\sqrt{\pi}}{2}(-\tau)^{1/2}\left[c_{1k}H_{\nu_{s}}^{(1)}\left(-k\tau\right)+c_{2k}H_{\nu_{s}}^{(2)}\left(-k\tau\right)\right]\,.$$

(3.6)

Adiabatic vacuum initial conditions for the observable range of Fourier modes of perturbations following from their quantization in the deep sub-Hubble limit $|k\tau|\gg 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 $\bf{k}$ including the limit $k\to 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

$$\tilde{u}_{k}=\frac{\mathcal{H}}{\sqrt{2k^{3}}}e^{ik\tau}(1-ik\tau)\,.$$

(3.7)

The curvature perturbation is defined as

$$\mathcal{R}=\Psi+H\frac{\delta R}{\dot{\bar{R}}}=\Psi-\frac{24H^{3}}{24H\dot{H}}\Psi\approx-\frac{1}{\epsilon}\Psi=-\frac{1}{2\epsilon}\frac{1}{\sqrt{3\mathcal{F}_{1}\bar{R}}}\frac{H}{\sqrt{2k^{3}}}e^{ik\tau}\left(1-ik\tau\right)\,.$$

(3.8)

and is (approximately) time-independent on super-Hubble scales $k\ll aH$. The primordial power spectrum and its tilt can be computed as

$$\mathcal{P}_{\mathcal{R}}=\left.\frac{H^{2}}{16\pi^{2}\epsilon^{2}}\frac{1}{3\mathcal{F}_{1}\bar{R}}\right|_{k=aH}\,,\quad n_{s}\equiv\left.\frac{d\ln\mathcal{P}_{\mathcal{R}}}{d\ln k}\right|_{k=aH}\approx 1-\frac{2}{N}\,,$$

(3.9)

where $N$ is the number of $e$-folds and $\epsilon=-\frac{\dot{H}}{H^{2}}\approx\frac{1}{2N}$. We can notice that the scalar spectral index remains exactly same as for the local $R^{2}$ model. From the Planck data [74] the power spectrum at the pivot scale $k_{\ast}=0.05\text{Mpc}^{-1}$ is $\mathcal{P}_{\mathcal{R}}^{\ast}\sim 2.2\times 10^{-9}$. Using this we get the mass of scalaron, Hubble parameter and Ricci scalar at $k_{\ast}=aH$ as

$$M\approx 1.3\times 10^{-5}\times\frac{55}{N_{\ast}}M_{p},\quad H\approx\sqrt{\frac{M^{2}}{6\epsilon}}\approx 5.6\times 10^{-5}\times\frac{55}{N_{\ast}}M_{p},\quad\bar{R}_{\ast}\approx 220M^{2}\times\frac{55}{N_{\ast}}\,.$$

(3.10)

Here $N_{\ast}$ is the number of $e$-foldings before the end of inflation for pivot scale $k_{\ast}$.

The second order action for the canonically normalized tensor perturbations in the de Sitter approximation (applying $\bar{R}\gg M^{2}$) is

$\displaystyle\delta^{2}S_{(T)}=$

$\displaystyle\frac{1}{2}\int d^{4}x\sqrt{-\bar{g}}h_{ij}^{\perp}e^{2\gamma_{T}\left(\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}\right)}\left(\bar{\square}-\frac{\bar{R}}{6}\right)h^{\perp{ij}}\,,$

(3.11)

where the subscript $(T)$ stands for tensor part and $\gamma_{T}$ here is an entire function⁴⁴4Note that in the notation of [11], $\omega\left(\square_{s}\right)=\gamma_{T}\left(\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}\right)$. which is related to the form-factor as

$$\mathcal{F}_{W}\left(\square_{s}\right)=\frac{\mathcal{F}_{1}\bar{R}}{\mathcal{M}_{s}^{2}}\frac{e^{\gamma_{T}\left(\square_{s}-\frac{2}{3}\frac{\bar{R}}{\mathcal{M}_{s}^{2}}\right)}-1}{\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}}\,.$$

(3.12)

Following the standard procedure for getting the tensor power spectrum and its tilt one yields [11]

	$\displaystyle\mathcal{P}_{\mathcal{T}}$	$\displaystyle=\left.\frac{1}{12\pi^{2}\mathcal{F}_{1}}\left(1-3\epsilon\right)e^{-2\gamma_{T}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right)}\right|_{k=aH}\,,$		(3.13)
	$\displaystyle n_{t}\equiv\left.\frac{d\ln\mathcal{P}_{\mathcal{T}}}{d\ln k}\right|_{k=aH}$	$\displaystyle\approx\left.-\frac{d\ln\mathcal{P}_{\mathcal{T}}}{dN}\left(1+\frac{1}{2N}\right)\right|_{k=aH}$
		$\displaystyle\approx\left.-\frac{3}{2N^{2}}-\left(\frac{2}{N}+\frac{3}{2N^{2}}\right)\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\gamma_{T}^{(\dagger)}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right)\right|_{k=aH}\,,$

The crucial difference here in comparison with local $R^{2}$ model is that the tensor power spectrum is scaled by an exponential factor of $\gamma_{T}$ evaluated at the pole of the tensor mode $\bar{\square}_{s}=\frac{\bar{R}}{6\mathcal{M}_{s}^{2}}$ and accordingly tensor tilt also gets modified with a term proportional to the derivative of $\gamma_{T}$ at $\bar{\square}_{s}=\frac{\bar{R}}{6\mathcal{M}_{s}^{2}}$.

The ratio of tensor-to-scalar power spectra is given by [11]

$$r=\left.\frac{12}{N^{2}}e^{-2\gamma_{T}\left(\frac{-\bar{R}}{2\mathcal{M}_{s}^{2}}\right)}\right|_{k=aH}\,.$$

(3.14)

In the local $R^{2}$ model $r=\frac{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 [73]. We can clearly notice that if $\gamma_{T}\left(\frac{-\bar{R}}{2\mathcal{M}_{s}^{2}}\right)=0$ one exactly recovers the same prediction as in the local $R^{2}$ model. A deviation from the local $R^{2}$ model happens if

$$\gamma_{T}\left(\frac{-\bar{R}}{2\mathcal{M}_{s}^{2}}\right)\Bigg{|}_{k=aH}\sim O\left(1\right)\,.$$

(3.15)

From (3.14) and (3.13) we can see that the values of $\left(r,\,n_{t}\right)$ explicitly depend on the ratio of scalar curvature to the square of $\mathcal{M}_{s}^{2}$ and the choice of entire function $\gamma_{T}\left(\square_{s}\right)$. The scalar curvature $\bar{R}$ during inflation can be read from (3.10). The two observables $\left(r,\,n_{t}\right)$ can be fixed by the parameter space

$$\Bigg{\{}-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}},\,\gamma_{T}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right),\,\gamma_{T}^{(\dagger)}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right)\Bigg{\}}\Bigg{|}_{k=aH}\,.$$

(3.16)

From the latest Planck 2018 data [58] we can deduce the following constraint (for $N_{\ast}=55$)

$$r<0.064\implies\gamma_{T}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right)\Bigg{|}_{k=aH}>-1.32\,.$$

(3.17)

We have no constraint on the tensor tilt however. Imposing a heuristic limit on tensor tilt $-4\lesssim n_{t}\lesssim 4$ we get a constraint on parameter space (3.16) as

$$-2N\lesssim-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\gamma_{T}^{(\dagger)}\left(-\frac{\bar{R}}{2\mathcal{M}_{s}^{2}}\right)\Bigg{|}_{k=aH}\lesssim 2N\,.$$

(3.18)

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

$$\gamma_{T}\left(\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}\right)\approx\left(\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}\right)\left[\frac{\alpha_{1}}{3}\left(\square_{s}-\frac{\bar{R}_{\ast}}{6\mathcal{M}_{s}^{2}}\right)+\alpha_{2}\left(\square_{s}-\frac{2\bar{R}}{3\mathcal{M}_{s}^{2}}\right)\right]\,.$$

(3.19)

where the approximation sign means that contributions of higher orders of the d’Alembertian operator are negligible when one evaluates this entire function and its derivative at $\square_{s}=\frac{\bar{R}}{6\mathcal{M}_{s}^{2}}$. 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 [75]. On the other hand this implies that given there are terms $O(\square_{s})$ in (3.19) one can freely choose signs of parameters $\alpha_{1,2}$.

In this particular example we see that $r$ depends on $\alpha_{2}$ and $\mathcal{M}_{s}$ because the coefficient of $\alpha_{1}$ vanishes when evaluated at $\square_{s}=\frac{\bar{R}}{6\mathcal{M}_{s}^{2}}$. But the tensor tilt $n_{t}$ depends on $\left(\alpha_{1},\,\alpha_{2},\,\mathcal{M}_{s}\right)$. In figures Figs. 2 and 3 we present the predictions for $\left(r,\,n_{t}\right)$ for some values of parameter space $\left(\alpha_{1},\,\alpha_{2}\right)$ and we take $\mathcal{M}_{s}\gtrsim H$ following the hierarchy in the scheme Fig. 1.

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 $\left(r,\,n_{s},\,n_{t}\right)$ 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 [76, 77, 78, 79, 80, 81, 82, 83]. In figure Fig. 4 we compare the predictions of non-local $R^{2}$-like model with the local $R^{2}$ model in $\left(n_{s},\,r\right)$ 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 [84, 85]. In many of the inflationary setups (except some special cases [86]) 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 [87].

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.

In this section, we compute the 3rd variation of the action (2.1) around (2.4) within the slow-roll 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 [60, 88, 64, 89, 90]

$$\langle\mathcal{R}\left(\mathbf{k_{1}}\right)\mathcal{R}\left(\mathbf{k_{2}}\right)\mathcal{R}\left(\mathbf{k_{3}}\right)\rangle=-i\int_{-\infty}^{\tau_{e}}ad\tau\langle 0|[\mathcal{R}(\tau_{e},\,\mathbf{k_{1}})\mathcal{R}(\tau_{e},\,\mathbf{k_{2}})\mathcal{R}(\tau_{e},\,\mathbf{k_{3}}),\,H_{int}]|0\rangle\,,$$

(4.1)

where $\boldsymbol{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 ($\mathcal{L}_{3}$) as $H_{int}\approx-\mathcal{L}_{3}$ that is a valid approximation during slow-roll inflationary regime [60, 64]. 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 $\tau\approx-\frac{1}{aH}$, it is a good approximation to calculate the integral in the limit $\tau_{e}\to 0$ [64]. In the Fourier space, the three-point correlation function of curvature perturbations is defined as

$\displaystyle\langle\mathcal{R}\left(\mathbf{k_{1}}\right)\mathcal{R}\left(\mathbf{k_{2}}\right)\mathcal{R}\left(\mathbf{k_{3}}\right)\rangle=\left(2\pi\right)^{3}\delta^{3}\left(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}}\right)\mathcal{B}_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)\,,$

(4.2)

where $B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)$ is called the bi-spectrum and $|\boldsymbol{k_{i}}|=k_{i}$. Due to the translational invariance, the total momentum $\mathbf{K}=\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}}$ is conserved. Scale invariance and rotational invariance imply that the 3-point function $B_{\mathcal{R}}\left(k_{1},\,k_{2},\,k_{3}\right)$ has to be homogeneous function of degree $-6$ that means

$$B_{\mathcal{R}}\left(\lambda k_{1},\,\lambda k_{2},\,\lambda k_{3}\right)=\lambda^{-6}B_{\mathcal{R}}\left(k_{1},\,k_{2},\,k3\right)\,.$$

(4.3)

and the number of independent variables of the 3-point function reduces to 2 which are the ratios $\frac{k_{2}}{k_{1}},\,\frac{k_{3}}{k_{1}}$ [91].

In the template of CMB observations, non-linear corrections in the curvature perturbation expressed as [92, 93]

$$\mathcal{R}=\mathcal{R}_{g}-\frac{3}{5}f_{NL}\left(\mathcal{R}_{g}^{2}-\langle\mathcal{R}_{g}\rangle^{2}\right)\,,$$

(4.4)

where $\mathcal{R}_{g}$ is the Gaussian random field and the $f_{NL}$ is the non-linearity parameter⁵⁵5The factor of $\frac{3}{5}$ comes from the relation between curvature perturbation $\mathcal{R}$ and the Bardeen variable $\Phi$ (Newtonian potential) during matter domination which reads as $\mathcal{R}=-\frac{5}{3}\Phi$,. Then

$$B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)=4\pi^{4}\frac{1}{\prod_{i}k_{i}^{3}}\mathcal{P}_{\mathcal{R}}^{2}A_{\mathcal{R}}\left(k_{1},\,k_{2},\,k_{3}\right)\,,$$

(4.5)

where $A_{\mathcal{R}}\left(k_{1},\,k_{2},\,k_{3}\right)$ is called the amplitude of bi-spectrum. In the momentum space, $f_{NL}$ is a function of three wave numbers and it is related to the amplitude of bi-spectrum as

$$f_{NL}=-\frac{5}{6}\frac{A_{\mathcal{R}}\left(k_{1},\,k_{2},\,k_{3}\right)}{\sum_{i}k_{i}^{3}}\,.$$

(4.6)

The $f_{NL}$ 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⁶⁶6Usually, 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^{2\mathcal{R}}\delta_{ij}dx^{i}dx^{j}$ [88]. 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 presented details of computations and results in the Appendix C.

• •

  We use the fact that $\Phi+\Psi\approx 0$ during the inflationary period as it was shown in [10, 11] and noted in the previous Section.

• •

  From the second order action for canonical variable $\Upsilon$, we can deduce that

  $$\square_{s}\Psi\approx\frac{M^{2}}{\mathcal{M}_{s}^{2}}\Psi\implies\square_{s}^{n}\Psi\approx\left(\frac{M}{\mathcal{M}_{s}}\right)^{2n}\Psi\,.$$

  (4.7)

  within the quasi-de Sitter approximation. In other words, we treat $\Psi$ 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 [60].

• •

  We carefully keep terms up to the leading order in $\epsilon$ and neglect higher order contributions treating $\epsilon\approx\mathrm{const}$. In this regard, to convert the 3rd order action in $\Psi$ into curvature perturbation $\mathcal{R}$, we make substitution $\Psi\approx-\epsilon\mathcal{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)

  	$\displaystyle\bar{\square}_{s}\mathcal{R}\approx$	$\displaystyle\,\frac{M^{2}}{\mathcal{M}_{s}^{2}}\mathcal{R}\implies$	$\displaystyle\mathcal{O}\left(\bar{\square}_{s}\right)\mathcal{R}\approx$	$\displaystyle\,\mathcal{O}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)\mathcal{R}\,$		(4.8)
  	$\displaystyle\bar{\square}_{s}\mathcal{R}^{\prime}\approx$	$\displaystyle\,\left(\frac{\bar{\square}_{s}}{\mathcal{M}_{s}^{2}}+\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\mathcal{R}^{\prime}\implies$	$\displaystyle\mathcal{O}\left(\bar{\square}_{s}\right)\mathcal{R}^{\prime}\approx$	$\displaystyle\,\mathcal{O}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}+\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\mathcal{R}^{\prime}\,,$

  where $\mathcal{O}\left(\square_{s}\right)$ 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 [94]

  $$\nabla_{\mu}\square_{s}\phi=\square_{s}\nabla_{\mu}\phi-\frac{R_{\mu\nu}}{\mathcal{M}_{s}^{2}}\nabla^{\nu}\phi\,.$$

  (4.9)

  where $\phi$ is a scalar. In the slow-roll approximation, $\bar{R}_{\mu\nu}\approx\frac{\bar{R}}{4}g_{\mu\nu}$, so we can approximate it as

  $\displaystyle\nabla_{\mu}F(\square_{s})\phi\approx$

  $\displaystyle\,\mathcal{F}\left(\square_{s}-\frac{R}{4\mathcal{M}_{s}^{2}}\right)\nabla_{\mu}\phi\,.$

  (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

  		$\displaystyle\int d^{4}x\sqrt{-\bar{g}}\delta R\delta\mathcal{F}\left(\square_{s}\right)\delta R\approx 4\epsilon^{2}\int d^{4}x\sqrt{-\bar{g}}\bar{R}\mathcal{R}\sum_{n=0}^{\infty}f_{n}\sum_{l=0}^{n-1}\bar{\square}_{s}^{l}\delta\square_{s}\bar{\square}_{s}^{n-l-1}\bar{R}\mathcal{R}\,$		(4.11)
  		$\displaystyle\overset{\bar{R}=\mathrm{const}}{\approx}4\epsilon^{2}\int d^{4}x\sqrt{-\bar{g}}\bar{R}\mathcal{R}\delta\square_{s}\mathcal{Z}_{1}\left(\bar{\square}_{s}\right)\bar{R}\mathcal{R}\overset{\bar{R}=\mathrm{const}}{=}0\,.$

  Looking carefully at the above derivation, we can see that passing from the first line to the second one, our assumption of approximating $\bar{R}\approx\mathrm{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 $\mathcal{Z}_{1}\left(\bar{\square}_{s}\right)\mathcal{R}=\mathcal{F}^{(\dagger)}_{R}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}\right)\mathcal{R}=0$ that follows from (2.6). Thus, we do not treat $\bar{R}=\mathrm{const}$ in the second line, rather we apply the background solution $\bar{\square}\bar{R}=M^{2}\bar{R}$ and $\dot{\bar{R}}\approx-2\bar{R}H\epsilon$, 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}/|\dot{H}|$, 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\ll aH$, the three point correlations of $\mathcal{R}$ do not evolve there, too. This means that they keep information of primordial interactions of scalar modes during the period $-\infty<\tau<\tau_{e}$ where $\tau_{e}\sim-\frac{1}{K}$ is a (conformal) time after few e-foldings of the Hubble radius exit [95].

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 $\mathcal{R}$ of AID gravity in the leading order slow-roll approximation is

	$\displaystyle\delta^{(3)}S_{(S)}=$	$\displaystyle 4\epsilon M_{p}^{2}\int d\tau d^{3}x\Bigg{\{}T_{1}\mathcal{R}\nabla\mathcal{R}\cdot\nabla\mathcal{R}+T_{2}\mathcal{R}\mathcal{R}^{\prime 2}+T_{3}\mathcal{H}^{2}\mathcal{R}^{3}$		(4.12)
	$\displaystyle+$	$\displaystyle T_{4}\mathcal{H}\mathcal{R}\mathcal{R}\mathcal{R}^{\prime}+T_{5}\mathcal{H}^{-1}\nabla\mathcal{R}\cdot\nabla\mathcal{R}\mathcal{R}^{\prime}+T_{6}\mathcal{H}^{-1}\mathcal{R}^{\prime 3}+T_{7}\mathcal{H}^{-2}\mathcal{R}^{\prime}\nabla\mathcal{R}\cdot\nabla\mathcal{R}^{\prime}\Bigg{\}}\,,$

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:

	$\displaystyle A_{\mathcal{R}}\left(k_{1},\,k_{2},\,k_{3}\right)=$	$\displaystyle\,-\left(2\epsilon+\frac{3\epsilon^{2}}{4}\right)A_{1}+\left(2\epsilon+\frac{3\epsilon^{2}}{4}\right)A_{2}-\frac{\epsilon^{2}}{2}A_{3}$		(4.13)
		$\displaystyle+\frac{\bar{R}}{M_{p}^{2}}\epsilon^{3}T_{\textrm{NL}}\Bigg{[}\frac{16}{3}\epsilon A_{2}-32A_{4}+2A_{5}-2A_{6}+\frac{16}{3}\epsilon A_{7}\Bigg{]}$
		$\displaystyle+\frac{4\bar{R}^{2}}{M_{p}^{2}\mathcal{M}_{s}^{2}}\epsilon^{4}\mathcal{F}^{(\dagger)}_{R}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{(}A_{2}+\frac{1}{3}A_{7}\Bigg{)}\,,$

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)

	$\displaystyle A_{1}=$	$\displaystyle 2\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}\left[K-\frac{k_{1}k_{2}+k_{2}k_{3}+k_{3}k_{1}}{K}-\frac{k_{1}k_{2}k_{3}}{K^{2}}\right]+\textrm{perms}\,,$		(4.14)
	$\displaystyle A_{2}=$	$\displaystyle\frac{2k_{1}^{2}k_{2}^{2}}{K}+\frac{2k_{1}^{2}k_{2}^{2}k_{3}}{K^{2}}+\textrm{perms}\,,$
	$\displaystyle A_{3}=$	$\displaystyle-\frac{K^{3}}{3}+2Kk_{1}k_{2}+\frac{k_{1}k_{2}k_{3}}{3}+\text{perms}\,$
	$\displaystyle A_{4}=$	$\displaystyle k_{3}^{2}\left[-\frac{4K}{3}-\frac{2k_{1}k_{2}}{K}-\frac{2k_{1}+2k_{2}}{3}\right]+\text{perms}\,,$
	$\displaystyle A_{5}=$	$\displaystyle 2\left(\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}\right)k_{3}^{2}\left[\frac{2}{K}+\frac{2k_{1}+2k_{2}}{K^{2}}+\frac{4k_{1}k_{2}}{K^{3}}\right]+\text{perms}\,,$
	$\displaystyle A_{6}=$	$\displaystyle\frac{12k_{1}^{2}k_{2}^{2}k_{3}^{2}}{K^{3}}\,,$
	$\displaystyle A_{7}=$	$\displaystyle\left(\boldsymbol{k}_{2}\cdot\boldsymbol{k}_{3}\right)k_{1}^{2}k_{3}^{2}\left[-\frac{2}{K^{3}}-\frac{6k_{2}}{K^{4}}\right]+\text{perms}\,,$

where $K=k_{1}+k_{2}+k_{3}$ and

	$\displaystyle T_{\textrm{NL}}$	$\displaystyle=\Bigg{[}\mathcal{F}_{R}\left(\frac{M^{2}}{\mathcal{M}_{s}^{2}}+\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)-\mathcal{F}_{1}\Bigg{]}\Bigg{|}_{K=aH}$		(4.15)
		$\displaystyle\approx\Bigg{[}\mathcal{F}_{R}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)-\mathcal{F}_{1}+\epsilon\frac{\bar{R}}{8\mathcal{M}_{s}^{2}}\mathcal{F}^{(\dagger)}_{R}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{]}\Bigg{|}_{K=aH}\,.$

The non-Gaussianity parameter $f_{NL}$ 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 [60, 96], 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 [96, 97]. 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⁷⁷7Note that the sound speed of curvature perturbation in our case is unity, so our result is fundamentally very different from single field non-canonical models [96].. 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_{NL}$ for squeezed $(k_{1}\to 0,k_{2}=k_{3}=\frac{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_{\textrm{NL}}\to 0$ and $\mathcal{F}_{R}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\to 0$. Especially, we can verify here the result of squeezed shape ($k_{1}\ll k_{2}\sim k_{3}$)

$$f_{NL}\approx\frac{5}{12}\left(1-n_{s}\right)+O\left(\epsilon^{2}\right)\,,$$

(4.16)

which is exactly the Maldacena consistency relation found in the single canonical scalar field inflation [60]. Note here that the Maldacena consistency relation was alternatively derived on kinematic grounds and was proven to be valid in general for single clock slow-roll models of inflation [61, 63]. The proof is based on the fact that in the standard slow-roll inflation the long wavelength mode that is far outside of the horizon effectively implies a renormalization of the background field and the corresponding rescaling of the short wavelength mode [61, 88]. However, it is known that even in the case of a single field slow-roll inflation, the consistency relation can be violated if the second, time dependent, mode of curvature perturbations ${\cal R}$ (often dubbed as the decaying mode) becomes important temporarily, so that the total ${\cal R}$ can not be considered as a constant during this period [98, 99, 100]. In this connection, there is a very instructive calculation of the bi-spectrum for the single field inflationary model [101] that has $n_{s}=1$ exactly, i.e. outside the slow-roll approximation. In this case, contrary to the zero non-Gaussianity parameter value predicted by the Maldacena consistency relation, $f_{NL}$ is non-zero and strongly ($\propto k_{1}^{2}$) scale-dependent [98]. This example shows also that the violation of the consistency condition for the bi-spectrum does not necessarily lead to the appearance of any new features in the spectrum itself.

In our case, the framework differs from standard single clock models. Even though still we have the slow-roll regime, the non-local nature of gravity leads to considerable effects at non-linear level where we have new scale-dependent interactions between different modes depending on the background scalar curvature and the form factor $\mathcal{F}_{R}\left(\square\right)$. In the case of local $R^{2}$ theory interactions between different modes is slow-roll suppressed and leads to the consistency relation [60] but in the case of non-local $R^{2}$-like inflation, the non-linear evolution of the curvature perturbation gets a significant scale dependence due to non-locality, especially around the scales of Hubble radius exit $k\sim aH$ when $\bar{R}\gtrsim\mathcal{M}_{s}^{2}$, leading to violation of consistency relation. This means that we get enhanced scale dependent interactions between the intermediate modes $\mathcal{R}$ with $k\sim aH$ and long wavelength modes with $k\ll aH$. Being more precise, at the second order action level (3.5) we can perform a field redefinition and do canonical normalization that gives the nearly scale invariant power spectrum that matches with the local theory exactly. At the level of the third order action (4.12) we invoke the effects of strong scale dependent non-local interactions (especially those interactions with time derivatives of different curvature modes⁸⁸8Interactions without any time derivatives of curvature perturbation are slow-roll suppressed, therefore, we have dropped them in our computation since we are interested in only leading order contributions.) between different modes in the next to leading order approximation in slow-roll as discussed around (4.7)-(4.11) and in Appendices B and C. It is crucial to emphasize that the key role in the modification of the consistency relation comes from these scale-dependent interaction terms whose effects are negligible in the regime when $\bar{R}\ll\mathcal{M}_{s}^{2}$ but relevant in the regimes of $\bar{R}\gtrsim\mathcal{M}_{s}^{2}$. We can see that the presence of non-local interactions of the different modes in our case gives a result for non-Gaussianity somewhat similar to one in a quasi-single field inflation [70] where heavy fields with masses of the order of Hubble parameter interact during inflation with curvature perturbation. These novel non-local interactions result in various scale dependent shapes (C.40)-(C.45) which do not resembles ones known in the literature [91]. Indeed, these non-standard interactions yield the violation of the consistency relation.

Below we discuss standard shapes of the reduced bi-spectrum which are squeezed, equilateral and orthogonal configurations. As explained above the shapes in question can be discriminated in the non-local $R^{2}$-like inflation in comparison to the original local $R^{2}$ model as well as other general scalar field(s) inflationary scenarios [91, 88, 66] due to non-local nature of gravity. Considering the general structure of the form-factor in (3.3), we obtain three shapes of reduced bi-spectrum $f_{NL}$ in the leading order approximation as

	$\displaystyle f_{NL}^{sq}$	$\displaystyle\approx\Bigg{\{}\frac{5}{12}\left(1-n_{s}\right)-23\epsilon^{2}\left[e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)}-1\right]-\frac{4\bar{R}}{\mathcal{M}_{s}^{2}}\epsilon^{3}e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)}\gamma_{S}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{\}}\Bigg{|}_{k=aH}\,$		(4.17)
	$\displaystyle f_{NL}^{eq}$	$\displaystyle\approx\Bigg{\{}\frac{5}{12}\left(1-n_{s}\right)-49\epsilon^{2}\left[e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{S}^{2}}\right)}-1\right]-\frac{9\bar{R}}{\mathcal{M}_{s}^{2}}\epsilon^{3}e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)}\gamma_{S}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{\}}\Bigg{|}_{k=aH}$
	$\displaystyle f_{NL}^{ortho}$	$\displaystyle\approx\Bigg{\{}\frac{5}{12}\left(1-n_{s}\right)-43\epsilon^{2}\left[e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)}-1\right]-\frac{22\bar{R}}{3\mathcal{M}_{s}^{2}}\epsilon^{3}e^{\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)}\gamma_{S}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{\}}\Bigg{|}_{k=aH}\,.$

From the above result for $f_{NL}$ we can deduce that the reduced bispectrum in the non-local $R^{2}$-like inflation depends not only on the slow-roll parameter $\epsilon$ but also on the background quasi-de Sitter curvature $\bar{R}$ and the non-local interactions that arise due to the perturbative expansion of form factor $\mathcal{F}_{R}\left(\bar{\square}\right)$ as discussed above in this Section. We can also notice that our results for the reduced bispectrum in (4.17) indicate a possible running of reduced bispectrum $f_{NL}$. For the time being we postpone the corresponding study for future investigations. In (4.17), we have presented only the leading order terms in $\epsilon$ 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

$$\Biggl{\{}\frac{\bar{R}}{4\mathcal{M}_{s}^{2}},\,\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right),\,\gamma_{S}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Biggr{\}}\Bigg{|}_{k=aH}$$

(4.18)

From (4.17) we notice that if $\gamma_{S}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Big{|}_{k=aH}=\gamma_{S}^{(\dagger)}\left(\frac{\bar{R}}{4\mathcal{M}_{s}^{2}}\right)\Bigg{|}_{k=aH}=0$, we recover the result of local $R^{2}$ theory in the leading order [88].