Mimetic Inflation

We study inflationary solution in an extension of mimetic gravity with the higher derivative interactions coupled to gravity. Because of the higher derivative interactions the setup is free from the ghost and gradient instabilities while it hosts a number of novel properties. The dispersion relation of scalar perturbations develops quartic momentum correction similar to the setup of ghost inflation. Furthermore, the tilt of tensor perturbations can take either signs with a modified consistency relation between the tilt and the amplitude of tensor perturbations. Despite the presence of higher derivative interactions coupled to gravity, the tensor perturbations propagate with the speed equal to the speed of light as required by the LIGO observations. Furthermore, the higher derivative interactions induce non-trivial interactions in cubic Hamiltonian, generating non-Gaussianities in various shapes such as the equilateral, orthogonal, and squeezed configurations with observable amplitudes.


Introduction
Mimetic gravity is a novel scalar-tensor theory proposed by Chamseddine and Mukhanov [1] as a modification of General Relativity (GR). The idea is to express the physical metric g µν in the Einstein-Hilbert action by performing a conformal transformation g µν = −(g αβ ∂ α φ∂ β φ)g µν from an auxiliary metricg µν in which φ is a scalar field. As a result, the longitudinal mode of gravity becomes dynamical even in the absence of any matter source. The above transformation can also be considered as a singular limit of the general disformal transformation where the transformation is not invertible [2,3]. With the physical metric, the scalar field is subject to the constraint 1 , g µν ∂ µ φ∂ ν φ = −1. (1.1) As a consequence of this constraint, the theory mimics the roles of cold dark matter in cosmic expansion, hence the theory is dubbed as the mimetic dark matter. The original mimetic model was then extended to inflation, dark energy and also theories with non-singular cosmological and black hole solutions [4][5][6]. See also Refs.  for further theoretical developments in mimetic gravity.
The original version of the mimetic theory is free from instabilities [37,38], but there is no nontrivial dynamics for scalar-type fluctuations. In order to circumvent this problem, the higher derivative term ( φ) 2 is added to the original action which generates a dynamical scalar degree of freedom propagating with a nonzero sound speed [4,7]. In addition, the mimetic model with a general higher derivative function in the form f ( φ) has been considered in Refs. [5,6]. However, these extended mimetic setups with a propagating scalar degree of freedom are plagued with the ghost and the gradient instabilities [39][40][41] 2 . To remedy these issues, it was suggested in [49][50][51] to extend the mimetic model further by considering direct couplings of the higher derivative terms to the curvature tensor of the spacetime such as φR, ∇ µ ∇ ν φR µν , φ∇ µ φ∇ ν φR µν and so on. By appropriate choices of these higher derivative couplings one can bypass the problems of the gradient and the ghost instabilities. However, now the background dynamics is more complicated and a simple dark mater solution is not a direct outcome of the analysis.
In this work our goal is to construct inflationary solutions in the extended mimetic setup with the effects of higher derivative couplings taking into account. As we will see the presence of higher derivative couplings to gravity generate new interactions and the analysis of cosmological perturbations become non-trivial. For example, because of these higher derivative interactions the dispersion relation of scalar perturbations receive higher order corrections resembling non-relativistic dispersion relation as in ghost inflation setup [52]. In addition, the predictions for the tensor perturbations are modified with a new consistency condition between the scalar spectral index n s , the sound speed of scalar perturbations c s and the tensor to scalar ratio r t .
Because of the higher derivative interactions, the model predicts novel non-Gaussianity features. The situation here is somewhat similar to the EFT studies of higher derivative corrections to the single field model [53] where large non-Gaussianity of various shapes such as equilateral and orthogonal types can be generated. In addition, similar to models with a nonstandard kinetic energy such as DBI model [54], the sound speed of scalar perturbations play non-trivial roles in generating large non-Gaussianities. The strong observational bounds on primordial non-Gaussianities can be used to constrain the model parameters. More specifically, the amplitude of non-Gaussianity parameter f NL in the squeezed, equilateral and orthogonal configurations from the Planck observations [55,56] are constrained to be We shall use these bounds to constrain model parameters and various couplings. The organization of the paper is as follows. In next Section, we present our setup and construct the background solutions which mimic cold dark matter even in the absence of normal matter. Then we extend these analysis to obtain an inflationary solution. In Section 3 we obtain the power spectrum of the curvature and tensor perturbations and calculate various cosmological observables. In Section 4, we study bispectrum and calculate the non-Gaussianity parameter f NL for local, equilateral and orthogonal configurations numerically, followed by discussions and summaries in Section 5. Many technical analysis of cosmological perturbations associated to power spectra and scalar bispectrum are relegated to Appendices A and B.

Inflationary Solution
In this section we study the background dynamics to obtain a period of inflation in early universe.
As summarized in Introduction, to remedy the ghost and gradient instabilities various higher derivative terms are added to the mimetic setup. Besides the higher derivative terms such as ( φ) 2 and ∇ µ ∇ ν φ∇ µ ∇ ν φ, we also require the higher derivative couplings of of the mimetic field to the curvature of the spacetime such as φR, ∇ µ ∇ ν φR µν , φ∇ µ φ∇ ν φR µν and so on [49][50][51]. Here we restrict ourselves to the simplest case where there is only a direct coupling of the higher derivative term φ ≡ χ to the Ricci scalar as follows, in which M P is the reduced Planck mass 3 . The Lagrangian multiplier λ enforces the constraint Eq. (1.1) [58]. In addition, we have allowed a potential term for the mimetic field which will drive inflation. In this setup P and F are arbitrary smooth function of χ. The former is added to make the scalar perturbation propagating (i.e. inducing a non-zero c s ) while the latter is required to remedy the gradient and the ghost instabilities [49][50][51]. As mentioned before, more complicated function of derivatives of φ and χ such as F 2 (∇ µ ∇ ν φR µν ) and F 3 (χ∇ µ φ∇ ν φR µν ) can also be added along with the simple function F (χ). However, the analysis even in the simplest setup of action (2.1) is complicated enough so we do not consider models with other higher derivative couplings. Before presenting the fields equation one important comment is in order. In the action (2.1) the effective gravitational coupling (effective reduced Planck mass) is actually M P F (χ) 1/2 . We can perform the calculations in the given "Jordan frame" but with a proper interpretation of the physical gravitational coupling. Alternatively, we may perform a metric field redefinition and go to the "Einstein frame" where the gravitational coupling is simply M P . The latter is rather complicated as the model presented in action (2.1) contains various higher derivative terms. Instead, we follow the first approach and work in the original Jordan frame. However, we make a further assumption that at the end of inflation the fields φ becomes trivial with F (χ e ) = 1 and one recovers the standard GR afterwards. This is a simplification made based on the intuitive ground though we do not have a dynamical mechanism to enforce it. We leave it as an open question as how or whether this transition from a mimetic setup to a standard GR setup can be achieved at the end of inflation. With these discussions in mind, we set M P = 1 in the rest of the analysis.
By taking the variation of the action (2.1) with respect to the inverse metric g µν , one obtains the Einstein field equations as F (χ)G µν = T µν where G µν is the Einstein tensor and T µν is the effective energy momentum tensor, given by in which P χ ≡ ∂P/∂χ and so on. In obtaining the above expression we have implemented the mimetic constraint (1.1). Clearly the energy momentum tensor given above can not be cast into the form of the energy momentum tensor of a perfect fluid. Moreover, varying the action (2.1) with respect to the scalar filed φ gives the following modified Klein-Gordon equation, The background cosmological solution is in the form of FRLW universe with the metric in which t and a are the cosmic time and the scale factor respectively. One can check that the background field equations take the following forms and where H =ȧ(t)/a(t) is the Hubble expansion rate. Note that at the background level, the mimetic constraint (1.1) enforcesφ = 1 and correspondingly χ = −3H. Before constructing the inflationary solution, let us for the moment set V = 0 to see how the mimetic setup can yield the dark matte solution. Using Eq. (2.3) the Lagrangian multiplier is obtained as follows in which C is an integration constant. Plugging this into the Friedmann equation (2.5), one obtains The constant C above indicates that a dark matter type solution can exist. However, in order for this conclusion to be valid we require the combination in the big bracket in Eq. (2.8) to be a constant. Since this combination plays important roles in our analysis below let us define (2.9) Using this definition of the function K and taking the time derivative of Eq. (2.8) once more, we obtain For a dark matter solution with a ∼ t 2/3 , the left hand side of Eq. (2.10) vanishes so indeed to have a dark matter solution we require the function K to be constant. In this case the Friedmann equation simplifies to 3H 2 =C/a 3 = ρ where ρ is the effective energy density and C ≡ 2C/K. To have a consistent solution one requires thatC > 0 while there is no restriction on the signs of C and K separately at this level. For example, in the original mimetic setup with F = −1 and P = 0, one obtains K = −1 and the dark matter solution is a direct outcome of the analysis. In conclusion, while the functions P (χ) and F (χ) are arbitrary, but in order to obtain a dark matter candidate in this setup one requires the combination K defined in Eq. (2.9) to be a constant. Now we consider the case when V (φ) = V (t) = 0 in order to obtain inflationary solution. Starting with the second Einstein equation (2.6) and using the definition of K to eliminate the combinations containing F χχ and P χχ in favours ofK we obtain In particular, if we set V = 0, we recover Eq. (2.10) as expected. Now noting that we can integrate Eq. (2.11) to obtain in which as in previous case C is a constant.
In an inflationary background one can neglect the constant term C in Eq. (2.13) as it is diluted rapidly. In addition, if V is nearly flat as in conventional slow-roll scenarios then we can obtain a phase of near dS spacetime with H 2 −V /K. Now we see the curious effect that in order to obtain an inflationary background we require the signs of V and K to be opposite. However, as we shall see in nest section, in order to have healthy scalar and tensor perturbations we require K > 0. As a result, to obtain an inflationary solution in this setup we need a negative potential. This should be compared with the analysis in [4] where F = 1, P (χ) = γ 2 χ 2 and V > 0. For these values of F and P we obtain K = −1 + 3γ/2. On the other hand, the sound speed of scalar perturbations in the model of [4] is c 2 s = γ/(2 − 3γ) so K has opposite sign compared to c 2 s (with γ > 0). Now in order to avoid the gradient instability one requires c 2 s > 0 so in their inflationary solution they need K < 0. But as we shall see in next section the sign of the quadratic action for the scalar perturbations is proportional to the sign of K so a negative K indicates the propagation of ghost as pointed out in details in [40]. As we mentioned above, this problem arises because in the analysis of [4] they considered V > 0 to construct inflationary solution.
Although the mimetic field φ is not an ordinary "rolling" scalar field in the sense that it appears with a constraint in the setup, but the requirement of a negative potential looks unexpected. However, we remind that negative potentials have been employed in the past in other contexts such as in contracting universes [59][60][61][62][63], see also [64,65].
To construct a specific inflationary setup, we consider the inverted quadratic potential as follows Inflation occurs when t < 0 while the hot big bang phase follows inflation for t > 0. We obtain the inflationary phase as the field rolls up the negative potential toward the origin. In addition, we assume that the potential vanishes for t > 0, so as discussed below Eq. (2.10), the seeds of observed dark matter can be obtained in this setup when inflation ends. So far our analysis were general and we did not specify the forms of the functions F (χ) and P (χ), unless we require the combination K to be positive. From now one, we further demand that K is a constant which simplifies the construction of the inflationary solution greatly. Now by introducing the new variable y ≡ a 3 2 , Eq. (2.11) becomes a linear differential which is similar to Eq. (31) in Ref. [4]. The inflationary branch of the solution is given by is the modified Bessel function of the second kind. At large negative t, one finds that the scale factor grows as whereas it is proportional to (t) 2/3 for positive t after inflation which is the indication of a dark matter dominated universe. Of course, we have to include reheating where radiation should be generated after inflation. The behaviours of the corresponding slow-roll parameters, For the future uses, let us define the slow-roll parameters associated with the background functions such as H, F etc as follows Note that the minus sign is chosen when X = H while for other background functions we chose the plus sign. Using the scale factor (2.17) in the inflationary phase, we find the following relations among the Hubble slow-roll parameters: which are satisfied for N ∼ 50 − 60 number of e-folds before the end of inflation. As mentioned before, our functions F (χ) and P (χ) are arbitrary except that we have imposed that the combination K defined in Eq. (2.9) to be constant. For example, the following pair of the polynomial functions satisfy the constraint (2.9) with K = −1 + 3γ/2. In general case, the polynomial functions F (χ) = n=0 f n χ n , and P (χ) = n=0 p n χ n can be considered with p 0 = 0 and p 1 an arbitrary constant. If the rest of coefficients f n and p n satisfy the relations p 2 = f 0 3 and One open question in this setup is the issue of reheating after inflation. To be consistent with the big bang cosmology, the inflationary phase has to be followed by a hot radiation dominated background. In conventional slow-roll models this is achieved via the (p)reheating mechanism in which the inflaton field transfers its energy to the Standard Model (SM) particles and fields while oscillating in its global minimum. In our mimetic scenario the field φ is not a rolling field in the usual sense but instead it is a space-filling field with the profile φ = t. So in order to achieve reheating one has to modify the current setup and couple the mimetic field to the SM fields one way or another. This is an open question which deserves a separate study elsewhere. We also comment that in the current setup with the potential (2.14) a dark matter solution is inherited in the solution for t > 0 so one may only need reheating to generate the host radiation while the dark matter can come from the mimetic source.

Primordial Power Spectra
In this section we calculate the power spectra of the curvature and tensor perturbations. For this purpose we calculate the quadratic actions associated to these perturbations.
The details of the analysis of the quadratic actions are presented in Appendix A. The quadratic action for the comoving curvature perturbation R and the tensor perturbations γ ij is obtained to be in which we have defined the parameters ϑ and σ as and during inflation the sound speed of scalar perturbations c 2 s as Note that ϑ is dimensionless while σ has the dimension of length square.
In order for the perturbations to be free from the ghost and gradient instabilities we require that all three parameters ϑ, c 2 s and σ 2 to be positive. Correspondingly we require In particular note that if K < 0 then the scalar perturbations develop ghost instability. This is the reason why we needed to couple the higher derivative terms to gravity to cure the ghost and gradient instabilities in the original setup of mimetic gravity [49][50][51].

Scalar power spectrum
The quadratic action for the scalar perturbations from the quadratic action (3.1) in Fourier space can be written as in which the prime indicates the derivative with respect to the conformal time dτ = dt/a(t) and we have defined the canonically normalized field u ≡ zR with z ≡ √ ϑa. In a near de Sitter background where ϑ, c s and σ are approximately constant 4 we have z /z 2/τ 2 and the corresponding mode function equation is given by The above equation indicates that we are dealing with a modified dispersion relation. With K > 0 we have σ 2 > 0 and the dispersion relation (3.6) is known as the Corley-Jacobson dispersion relation which was studied for investigating the black holes physics [67,68] and for the effects of trans-Planckian physics on cosmological perturbations [69,70]. In addition, this type of dispersion relation occurs in ghost inflation [52] where a timelike scalar field fills the entire spacetime with the profile φ = t as in our mimetic setup. Such modified dispersion relations indicate the violation of Lorentz invariance in the UV limit. However, for low physical momentum when k a c s σ , the linear dispersion relation is recovered. We can define the scale at which the modification to the linear dispersion relation becomes important as Λ ≡ c s /σ. For the physical momentum larger than this scale, k phy Λ, the quartic contribution to the dispersion relation becomes important. For future purpose we introduce the parameter ν via which quantifies the ratio of the sound-Hubble horizon parameter over the momentum scale Λ around when the behaviour of the dispersion relation changes. For the models in which ν 1 the dispersion relation is a linear relation, ω ∝ k, as in standard slow-roll models, while for ν 1 the mode function u k is described by the non-relativistic dispersion relation ω ∝ k 2 .
Imposing the adiabatic vacuum initial conditions, the mode function of the comoving curvature perturbation is obtained to be [71][72][73] where a k and a † k are the creation and the annihilation operators as usual and W i 4ν , 3 4 is the Whittaker function.
With the help of the above mode function, it is easy to calculate the super horizon (c s kτ → 0) limit of the power-spectrum for comoving curvature perturbation. Taking into account the asymptotic behaviour of Whittaker function, i.e.
, the curvature perturbations power spectrum on superhorizon scales is given by (3.10) Let us now discuss about the asymptotic behaviour of the the power spectrum in small and large ν limits. In the limit ν 1 5 , we find out (3.11) Since in the limit ν 1 we have a relativistic dispersion relation, one expects that the power spectrum in this limit resembles that of standard slow-roll inflation. Indeed, if we formally identify the coefficient ϑ in the quadratic action (3.1) with the corresponding factor in the action of slow-roll models [75], ϑ ↔ 2 H /c 2 s , then the power spectrum in Eq. (3.11) reduces to the standard result P R = H 2 /8π 2 H c s in slow-roll models. On the other hand, in the limit ν 1 the quartic term k 4 dominates in ω 2 and the dispersion relation becomes non-relativistic as in the model of ghost inflation [52]. In this limit the power spectrum (3.10) reduces to (3.12) Identifying a suitable choice of the parameters ϑ and σ with the corresponding parameters in [52] we reproduce the power spectrum for ghost inflation as well.
Having calculated the curvature perturbation power spectrum, we can also calculate the spectral index n s as where the subscript * shows the time of horizon crossing for the mode of interest k and we have used our slow-roll notation (2.19) for the background variables X = c s , ϑ, g(ν). In order to have an almost scale invariant power spectrum, one requires the four parameters H , ϑ , cs , and g to be very small.

Tensor power spectrum
To calculate the power spectrum of tensor perturbations, let us first expand the tensor modes of the quadratic action (3.1) in terms of their polarization tensors e + ij and e × ij as γ ij = +,× γ λ e λ ij where λ = +, × and e λ ij are symmetric, transverse and traceless tensors. Moreover, using the normalization condition, e λ ij e λ ij = 2δ λλ , we obtain the second-order action for the tensor modes in Fourier space as follows wherez 2 ≡ F (χ) a 2 /2. In order for the perturbation to be stable, we require that F > 0.
Interestingly, from the above action we see that the tensor modes propagate with the speed equal to unity, c T = 1, i.e. the tensor perturbations propagate with the speed of light. This is because we considered the special case of higher derivative coupling to gravity in the form of F (χ). However, it is well-known that for general higher derivative interactions with gravity, c T is not equal to speed of light. These types of modified gravity theories are under strong constraints from the LIGO observations which require that |c T −1| < 5×10 −15 [76][77][78]. For example, in our setup if we allow more general higher derivative interactions such as the curvature independent quadratic higher derivative terms ∇ µ ∇ ν φ∇ µ ∇ ν φ and the curvature dependent cubic higher derivative terms φ∇ µ φ∇ ν φR µν and ∇ µ ∇ ν φR µν then c T = 1 [47].
Upon defining the canonically normalized field associated with γ λ by v λ ≡zγ λ and imposing the the Minkowski (Bunch-Davies) initial condition, the mode function is obtained to be Defining the power spectrum of the gravitational tensor modes via we obtain where both H and F are evaluated at the time of horizon crossing. Compared to conventional models of inflation, we see the additional factor 1/F in tensor power spectrum. This is understandable if one notes that, naively speaking, we have rescaled the gravitational coupling M 2 P → M 2 P F in the starting action Eq. (2.1). The spectral index of P γ is also given by where F is the slow-roll parameter associated with F as defined in Eq. (2.19) which is given by whereν ≡ K/F and we have used Eq. (3.3) in the last step. In particular, we see that n t depends on c s andν after plugging the slow parameter F from Eq. (3.19) into Eq. (3.18). In Fig. 2 we have presented the predictions for n t for some values of (c s ,ν) in the parameter space. Interestingly, we see that in some regions of parameter space n t > 0, i.e. the tensor power spectrum is blue-tilted. This is unlike the conventional slow-roll models which generally predict a red-tilted tensor power spectrum. As is the case in our model, the detection of a blue-tilted tensor perturbations cannot rule out inflation automatically [59,79,80].
As long as we assumeν O( √ H ), then Eq. (3.19) guarantees that F − H in the subluminal regime with 0 < c s < 1. It means that the function F (χ) changes very slowly during slow-roll inflation. Therefore, we can consider it approximately as a constant during inflation. To estimate this value, let us first define the tensor to scalar ratio as follows, . (3.20) Then, by restoring M P in the scalar and tensor power spectra and using the current observational constraint on inflationary parameters [55], i.e. r t 0.056 and P R 2.1 × 10 −9 , the value of F at horizon crossing can be estimated as which implies that we need to choose K 10 −5 to satisfy the conditionν O( √ H ). As mentioned before, the slow roll approximation F ≈ − H is valid only in the region where the comoving curvature perturbation propagates with c s < 1. The superluminal propagation speed with c s > 1 is not a problem per se as it does not directly violate causality on the background [81,82]. However, we restrict ourselves to scalar perturbations with subluminal speeds.
Using Eqs. (3.19) and (3.18), we can also obtain the following generalized consistency relation between r t and n t , We see that the consistency relation in conventional models of inflation [83], r t = −8c s n t , is modified in our model due to the mimetic constraint. In Fig. 3 we have presented the predictions for r t (ν,ν) for c s = 1. The white areas correspond to the regions of parameter space which are not allowed due to the observational bound r t ≤ 0.056 [55]. By choosing smaller values of c s the allowed regions become more extended.

Primordial Bispectra
In this section, we calculate the three-point correlation of the scalar perturbations RRR and look at the amplitudes and shapes of non-Gaussianity in various limits. Utilizing the standard methods, the expectation value of the three point correlation is given by [84] where H int is the interaction Hamiltonian which is calculated from expanding the Lagrangian (2.1) up to 3rd orders in curvature perturbations, given in (A. 16), with H int = −L 3 . Moreover, k i are the wave vectors and τ i is the initial time when the inflationary perturbations are deep inside the Hubble radius. Since during a quasi-de Sitter expansion τ −1/(aH), it is a good approximation to calculate the integral in the limit τ i → −∞ and τ e → 0.
In the Fourier space, we can write the three-point correlation function of curvature perturbations as in which k i = |k i | and B R (k 1 , k 2 , k 3 ) is called the bispectrum 6 which can be parameterized as where A is called the amplitude of bispectrum. Finally, the non-linearity parameter f NL associated with the amplitude of bispectrum is defined by the following relation As we see from Eq. (A.16), our interaction Hamiltonian contains 22 independent terms (interactions). These complicated interactions originate from the higher derivative terms in F (χ) and P (χ). Each of them induce different shapes and amplitudes of non-Gaussianities.
As examples, let us calculate the bispectrum for the following two terms of the cubic action (A. 16), which also exist in the model of ghost inflation with the modified dispersion relation ω 2 ∝ k 4 . 6 Because of the translational invariance, the total momentum K ≡ k1 + k2 + k3 is conserved.
The bispectrum for each term in Eq. (4.5) is evaluated using the mode function of R given in Eq. (3.9) as follows Using the explicit expression for the wave function (3.9) and substituting the above results into Eq. (4.3) for A, we obtain the following expressions for the amplitudes A (2) and A (8) associated with each interaction: with p 1 = 1 and p 2 = p 3 = 0, and Here the function I p 1 ,p 2 ,p 3 n 1 ,n 2 is defined via 11) and the upper index p i denotes the order of derivative with respect to the function variables. For example h and so on. The amplitudes for all other interactions are listed in Appendix B.
To study the shape function of the above amplitudes, in Figs. 4 and 5 we have presented the 3D plot of r −1 2 r −1 3 A(1, r 2 , r 3 ) as a function of r 2 ≡ k 2 /k 1 and r 3 ≡ k 3 /k 2 for ν = {1, 10, 50}. The plots are produced numerically, after rotating the contour of integration over τ along the direction ∝ −(1 + i) so that they converge exponentially. We see that A (2) and A (8) roughly have similar shapes and amplitudes and both roughly peak at the equilateral limit k 1 = k 2 = k 3 = k. In addition, the variation of ν has no significant effects on the shapes.
In Table 1 we list the shape of each contribution presented in Appendix B. One can see that most of the non-Gaussianity shapes peak at the equilateral limit where all three modes have comparable wavelengths. However, some shapes are close to the orthogonal shape and   the local shape which has a peak in the squeezed limit. For example, as shown in Fig. 6, A (4) and A (10) peak in the squeezed triangle limit (k 3 k 1 k 2 ) and in orthogonal triangle limit (k 3 = k 2 = k 1 /2), respectively.
Combining the contributions from all interactions listed in Appendix B, the total non-

Shape Local Equi Local Equi Equi Local Equi Ortho Equi Equi
Amplitude A (13) A (14) A (15) A (16) A (17) A (18) A (19) A (20) A (21) A (22) Shape Equi Local Equi Equi Equi Equi Equi Equi Ortho Equi Gaussianity parameter f NL is given by Correspondingly, we can calculate f NL numerically for squeezed (k 1 = k 2 = k, k 3 → 0), equilateral (k 1 = k 2 = k 3 = k) and orthogonal (k 1 = k, k 2 = k 3 = k/2) shapes. In Figs. 7, 8, and 9, f NL is presented in the various range of ν in the squeezed, equilateral and orthogonal configurations. It is worth mentioning that f NL is controlled by three parameters, the sound speed c s , the scalar to tensor ratio r t and ν. In the left hand panels of these figures, f NL can take the observationally allowed values in some range of ν by varying c s while r t = 0.01 is held fixed. A similar conclusion holds in the right hand panels where we fix c s = 1 and vary r t . Generally, f NL increases by reducing c s and r t . One can find corners of parameter space which yield to acceptable amplitudes for f NL as required by observations in Eq. (1.2).
For further studies of bispectrum and its expansion in terms of slow-roll parameters see Appendix B.1.

Summaries and Conclusions
In this paper we have studied inflationary solution in an extension of mimetic gravity with higher derivative interactions coupled to gravity. It is known that the original mimetic setup is plagued with the ghost and gradient instabilities. These instabilities can be removed with the help of higher derivative interactions coupled to gravity. There are a number of options to include higher derivative corrections. In this paper we have studied the simplest higher derivative correction in the form F (χ)R with χ ≡ φ. It would be interesting to extend the current analysis to include other higher derivative terms coupled to gravity such as F 2 (∇ µ ∇ ν φR µν ), F 3 (χ∇ µ φ∇ ν φR µν ) etc. In addition, in order for the scalar perturbations to become dynamical with a non-zero sound speed c s we have included the term P (χ) as well. One curious effect in our analysis is that in order to obtain an inflationary solution we have to work with a negative potential. This conclusion is a consequence of the fact that we deal with a constrained theory. More specifically, in oder for the quadratic actions of the scalar and tensor perturbation to be free from instabilities the higher derivative functions F (χ) and P (χ) are subject to certain conditions which cause the potential to be negative. Inflation is achieved while the field rolls up the potential towards V = 0. While a negative potential may be considered problematic a priori but our analysis show that the setup shows no pathologies either at the background or at the perturbation level.
While the background yields a period of slow-roll inflation the cosmological perturbations in this setup have novel behaviours. Because of the higher derivative interactions the dispersion relation associated with the scalar perturbations receives higher order momentum corrections as in the model of ghost inflation. Furthermore, the tilt of tensor perturbations can take either signs in contrast to conventional inflation models [84,85]. In addition, we obtain a new consistency relation between r t and n t which involves c s and other model parameters encoding the higher derivative interactions. Despite the presence of higher derivative corrections the tensor perturbations propagate with the speed equal to speed of light as strongly implied by the LIGO observations.
We also studied the predictions of this setup for the amplitudes and shapes of non-Gaussianities. Because of higher derivative interactions, various types of interactions are This result simplifies our following calculations considerably.

A.1 Linear perturbations: Quadratic action
Plugging the above perturbations into the action (2.1) and after some integration by parts, the quadratic action in comoving gauge for R and γ ij is obtained as follows where we have used Eq. (2.9) to simplify the result in terms of the function K. Note that in the these analysis there is no assumption on the function K so it is kept general. It is evident that the ψ mode is a non dynamical degrees of freedom which can be integrated out from the action. Varying Eq. (A.4) with respect to ψ, we find Substituting the above result into the action (A.4) and after some integration by parts the quadratic Lagrangian in comoving gauge for R and γ ij is obtained to be in which we have defined and during inflation, the sound speed of the scalar perturbations c 2 s as As discussed in the main text, in order to avoid the ghost and gradient instabilities we demand that ϑ > 0, c 2 s > 0 and σ 2 > 0. Applying these conditions on Eqs. (A.7) and (A.8), we obtain the following constraints

A.2 Nonlinear scalar perturbations: Cubic action
In this section, we calculate the cubic action for the scalar perturbations which is used to calculate the bispectrum.
Expanding the action (2.1) up to third order, the cubic action is given by where, using the definition of K in Eq. (2.9), the coefficientsf i are given bỹ The next step is to eliminate ψ in the above action by utilizing Eq. (A.5). To do this, let us define ∂ 2 Ξ ≡ Q 2Ṙ in which Q 2 ≡ 3ν 2 /(1 +ν 2 ) withν = K/F . This definition provides us with the contributions proportional to the linear differential equation of R, i.e., in the final cubic action. Substituting the relation into the action (A.10) and performing a lot of integrations by parts and dropping the total derivative terms 7 , the corresponding cubic Lagrangian is obtained to be 16) in which the coefficient in front of δL 2 /δR| 1 is Here ∂ −2 is the inverse Laplacian and κ ≡ F / H + H / H . Clearly, all contributions in F include time and spatial derivatives of R which vanishes in the large-scale limit (k → 0). When we calculate the bispectrum, we neglect the last term in cubic Lagrangian (A.16) relative to those coming from other terms. The other coefficients are given by of the cubic Lagrangian (A.16) as follows: For the case of p 1 = 1 and p 2 = p 3 = 0, we have Finally, for p 1 = p 2 = 1 and p 3 = 0, the rest of expressions are given by | ijl |I p i ,p j ,p l 1,0 k l (B.13) | ijl |I p i ,p j ,p l 0,−2 k i k j (B.14) (B.20)

B.1 Expansion in terms of slow-roll parameters
In order to derive a simple expression for f NL in the equilateral configuration to the order of the slow roll parameter H , let us first consider all terms defined in Eq. (2.19) to be much smaller than unity andν O( √ H ). Then the amplitude of non-Gaussianities can be written as A (j) = B j S (j) , (B.21) where B j are the coefficients coming in front of each shape function S (j) for the amplitude A (j) listed in previous subsection, for example B 1 = f 1 H 4 /(c 8 s P 2 R ϑ 3 ). With this decomposition, the shape coefficients can be expanded as shown in  Table 2. The expansion of the shape function coefficient B (i) in terms of the slow-roll parameter H .
Here we have defined η ≡ It worth mentioning that one can not discard the sub-leading orders of the slow-roll parameter relative to the leading order, because integral functions I p,q,r n,m andν are running with ν. It is interesting that the relative error in this approximation is This means that the expansion coefficients presented in Table. 2 are near to their exact values with high accuracies. Having these shape coefficients in hand, one can calculate f NL in other configurations by using the amplitudes presented in Appendix. B.