Spinning Guest Fields during Inflation: Leftover Signatures

Emanuela Dimastrogiovanni Matteo Fasiello and A. Emir Gümrükçüoğlu

Abstract

We consider the possibility of extra spinning particles during inflation, focussing on the spin-2 case. Our analysis relies on the well-known fully non-linear formulation of interacting spin-2 theories. We explore the parameter space of the corresponding inflationary Lagrangian and identify regions therein exhibiting signatures within reach of upcoming CMB probes. We provide a thorough study of the early and late-time dynamics ensuring that stability conditions are met throughout the cosmic evolution. We characterise in particular the gravitational wave spectrum and three-point function finding a local-type non-Gaussianity whose amplitude may be within the sensitivity range of both the LiteBIRD and CMB-S4 experiments.

1 Introduction

The existence of an early phase of accelerated expansion, inflation, is by now a centrepiece of our understanding of how the universe begun [1] and has come to assume its current structure. A scalar field with a sufficiently flat potential can drive inflation long enough to make sense of key observations, such as those of the cosmic microwave background (CMB) anisotropies. Within the single-field slow-roll paradigm, one may find several consistent realisations of a technically natural small inflaton mass (see e.g. [2]) and verify complete agreement with current observational bounds on inflation. These successes notwithstanding, ours remains a broad-brush picture of inflation. Key questions remain unanswered, including fundamental ones: what is the energy scale of inflation? what is the early phase field content?

The available constraints on inflationary parameters such as the tensor-to-scalar ratio $r$ and non-Gaussianity $f_{\rm NL}$ leave the door open to a plethora of intriguing possibilities, an embarrassment of riches one would like to sacrifice in favour of more a well-defined inflationary dynamics. The unprecedented array of cosmological probes scheduled for launch in the coming years promises to propel us in this very direction. The projected constraints on cosmological parameters will cross qualitative thresholds, including those that will help us solve the single vs multi-field dichotomy. One should mention here the target of $\sigma_{r}\sim 10^{-3}$ within reach for CMB-S4 experiments [3] and the LiteBIRD satellite [4], which will rule in (out) celebrated models such as Starobinsky inflation [5]. The constraints on (scalar) non-Gaussianity from Euclid [6], LSST [7] and SKA [8, 9] will be just as consequential: these may go as far as $\sigma_{f_{\rm NL}}\sim 1$ and down an order of magnitude or so [10] with 21cm cosmology. In typical single-field slow-roll realisations non-Gaussianities are much smaller than unity. On this basis, a near future detection of primordial non-Gaussianity may provide the “smoking gun” for multi-field inflation.

Remarkably, the last decade has witnessed the rise of new probes of gravity. Following the first direct gravitational wave (GW) detection in 2015, these probes have led to spectacular new insights into gravitational wave astronomy and hold the potential to transform also our understanding of early universe physics. Ground (and, soon, space-) based interferometers access very small scales thus opening up, together with pulsar timing arrays, an observational window on inflation that is complementary to both CMB and large scale structure (LSS) measurements. For example, the detection of a primordial GW signal at small scales would strongly point to multi-field inflationary mechanisms. Further characterisation of the GW signal may reveal much more specific information about inflation. A case in point is the possibility of a chiral GW spectrum, often associated to Chern-Simons-type couplings [11, 12].

The opportunities afforded by a multi-probe and multi-scale search for imprints of early universe physics call for a full characterisation of the signatures associated to the inflationary particle content, going well beyond the single-field hypothesis. The reasons for doing so are manifold. From the top-down perspective, it is clear that UV finite theories such as string theory point to a rich inflationary field content [13]. From moduli fields to gauge fields to Kaluza Klein modes, all can be 4-dimensional end products of higher dimensional constructions. A multi-field mechanism is not just a “natural” scenario but also an eminently testable one. Besides the discovery potential stemming from exploring primordial GW scale dependence and chirality, non-Gaussianity provides the most powerful handle on the inflationary particle content. From the momentum scaling, angular dependence, and possible periodic features of squeezed non-Gaussianities we can infer crucial information on the mass, spin, and couplings of the inflationary content [14, 15].

There is a rich and interesting literature on the model-building aspects and the signatures of extra scalar and vector fields during inflation [16]. In this work we shall focus instead on the possibility of higher spin particles, and in particular spin-2 fields. Our starting point is the relatively recent formulation [17, 18, 19] of a consistent fully non-linear theory of interacting spin two theories. These models have been widely used in attempts to deliver a technically natural mechanism to drive the universe late-time acceleration [20]. In our present work, we adopt a different perspective, namely that of studying the cosmology associated to the presence of an extra spin-2 field $f$ that couples to gravity (i.e. the usual massless spin-2 field $g$ that couples to matter) in a consistent manner. Crucially, already the results in [21] indicated that in a theory of interacting spin-2 particles there may be at most one massless spin-2 field. That place having been taken up by the massless field $g$ , the theory with both the $g$ and $f$ field will necessarily contain massive eigenstates.

We provide a thorough study of the cosmology of consistent interacting spin-2 theories in the inflationary context. The presence of different tensor sectors, such as those due to multiple spin-2 particles during inflation, offers the intriguing possibility of sourcing gravitational waves at linear order, in contradistinction to single-field scenarios. As we shall see, large GW non-Gaussianities are perhaps the most notable signature associated with the field content we explore in this manuscript.

This paper is organised as follows. In Section 2 we study the background solution and the perturbative stability of our set-up. We provide extensive details for the benefit of the reader who may not be familiar with the intricacies stemming from requiring a stable FRW cosmology for a theory with multiple spin-2 fields. Section 3 is devoted to the calculation of the GW power spectrum and the three-point function amplitude and scale dependence. In Section 4 we identify the allowed parameter space of the model and provide in Section 5 a concrete example with corresponding signatures. In Section 6 we offer our conclusions and point to future work. For completeness, we report in Appendix A the additional details necessary to reproduce our calculations.

2 The theory of extra spin–2 field with gravitational interactions

In this section we review the theoretical framework at the heart of our analysis and discuss the conditions that arise from a consistent cosmological evolution. We are after a theory of cosmic inflation driven by a scalar field but whose overall particle content is enriched by an extra spin-2 field. Fortunately, in describing such an early universe set-up, we can immediately benefit from an extensive literature on spin-2 particles at late times. The background being that of Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime, we will be able to exploit here much of the late universe results and “technology”.
We start with the notion that an additional interactive spin-2 particle in the inflationary context (i.e. with a scalar and a massless spin-2 field already in place) is necessarily a massive one [21]. The fully non-linear consistent formulation of a spin-2 massive theory was first put forward in [22, 17]. In our set-up the fiducial metric of [17] becomes dynamical, making the Hassan-Rosen Lagrangian [18] together with a “matter” sector Lagrangian $\mathcal{L}_{\phi}$ that describes the slowly rolling inflaton our ideal starting point:

S=\int d^{4}x\left[\frac{M_{g}^{2}}{2}\sqrt{-g}\,R[g]+\frac{M_{f}^{2}}{2}\sqrt{-f}\,R[f]-m^{2}M^{2}\sqrt{-g}\sum_{n=0}^{4}\beta_{n}\,e_{n}(\sqrt{g^{-1}f})+\sqrt{-g}\,\mathcal{L}_{\phi}\right]\;.

(2.1)

In the above, the matter follows the geodesics of the tensor $g_{\mu\nu}$ which we thus interpret as the physical metric. The second dynamical tensor field $f_{\mu\nu}$ represents the additional spin–2 field. The coupling between the two tensors $g$ and $f$ is given by the Lorentz invariant de Rham-Gabadadze-Tolley (dRGT) terms, which are constructed by the elementary polynomials $e_{n}$ . These are defined via the determinant expansion $\det(\mathbb{1}+X)=\sum_{n=0}^{4}e_{n}(X)$ , as

$\displaystyle e_{0}(X)$	$\displaystyle\equiv 1\,,$
$\displaystyle e_{1}(X)$	$\displaystyle\equiv[X]\,,$
$\displaystyle e_{2}(X)$	$\displaystyle\equiv\frac{1}{2!}\,\left([X]^{2}-[X^{2}]\right)\,,$
$\displaystyle e_{3}(X)$	$\displaystyle\equiv\frac{1}{3!}\,\left([X]^{3}-3\,[X]\,[X^{2}]+2\,[X^{3}]\right)\,,$
$\displaystyle e_{4}(X)$	$\displaystyle\equiv\frac{1}{4!}\,\left([X]^{4}+8\,[X]\,[X^{3}]-6\,[X]^{2}[X^{2}]+3\,[X^{2}]^{2}-6\,[X^{4}]\right)\,,$	(2.2)

where square brackets around tensorial quantities represent the trace operation. Following the convention of Ref. [23], we also define

\frac{1}{M^{2}}\equiv\frac{1}{M_{g}^{2}}+\frac{1}{M_{f}^{2}}\,,

(2.3)

such that the two-tensor interaction term is symmetric under the $g\leftrightarrow f$ exchange.

2.1 Review of background cosmology

For the background configuration, we aim at a flat FLRW cosmology in the $g$ -metric. We also restrict the $f$ –metric to be homogeneous and isotropic with a vanishing spatial curvature. The corresponding line elements are

	$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=-N^{2}dt^{2}+a^{2}\delta_{ij}dx^{i}dx^{j}\,,$
	$\displaystyle f_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=\xi^{2}\left(-N^{2}\tilde{c}^{2}dt^{2}+a^{2}\delta_{ij}dx^{i}dx^{j}\right)\,,$		(2.4)

where $\xi$ quantifies the ratio of the scale factors in the two metrics, while $\tilde{c}$ defines the light cone on the $f$ –metric. In order to determine the background dynamics, we keep explicit the freedom on time reparametrization. We will later fix the cosmological time by setting $N=1$ .

Combining the zero-order contributions from kinetic (A.3), interaction (A.7) and matter (A.12) terms, we can write down the total mini-superspace action as

S_{\rm mss}=V\,\int dt\,a^{3}N\left[-\frac{3\,M_{g}^{2}\dot{a}^{2}}{a^{2}N^{2}}-\frac{3\,M_{f}^{2}(\xi\,\dot{a}+a\,\dot{\xi})^{2}}{a^{2}\tilde{c}\,N^{2}}-m^{2}M^{2}(\rho_{m,g}+\tilde{c}\,\xi^{4}\rho_{m,f})+P\right]\,,

(2.5)

where $P$ represents the analogue pressure for the inflaton field. In the rest of the paper we use the following functions to represent convenient combinations of the $\beta_{n}$ parameters ¹¹1The definition of function $\Gamma$ in Ref. [24] differs from the one here by an additional term $(\tilde{c}-1)Q$ . The two definitions coincide in the limit $\tilde{c}=1$ .

$\displaystyle\rho_{m,g}\equiv$	$\displaystyle\beta_{0}+3\,\beta_{1}\xi+3\,\beta_{2}\xi^{2}+\beta_{3}\xi^{3}\,,$
$\displaystyle\rho_{m,f}\equiv$	$\displaystyle\beta_{4}+\frac{3\,\beta_{3}}{\xi}+\frac{3\,\beta_{2}}{\xi^{2}}+\frac{\beta_{1}}{\xi^{3}}\,,$
$\displaystyle\Gamma\equiv$	$\displaystyle\xi\,(\beta_{1}+2\,\beta_{2}\xi+\beta_{3}\xi^{2})\,,$
$\displaystyle Q\equiv$	$\displaystyle\xi^{2}(\beta_{2}+\beta_{3}\xi)\,.$	(2.6)

Making use of Eqs.(A.11), we vary the mini-superspace action (2.5) with respect to $N$ , $\tilde{c}$ , $a$ , $\xi$ and $\phi$ . We obtain he following set of equations of motion:

$\displaystyle 3\,M_{g}^{2}H^{2}$	$\displaystyle=m^{2}M^{2}\rho_{m,g}+\rho\,,$	(2.7)
$\displaystyle 3\,M_{f}^{2}H_{f}^{2}$	$\displaystyle=m^{2}M^{2}\rho_{m,f}\,,$	(2.8)
$\displaystyle 2\,M_{g}^{2}\,\dot{H}$	$\displaystyle=m^{2}M^{2}(\tilde{c}-1)\Gamma-(\rho+P)\,,$	(2.9)
$\displaystyle 2\,M_{f}^{2}\,\dot{H}_{f}$	$\displaystyle=-\frac{m^{2}M^{2}}{\xi^{3}}(\tilde{c}-1)\Gamma\,,$	(2.10)
$\displaystyle\dot{\rho}$	$\displaystyle=-3\,H\,(\rho+P)\,,$	(2.11)

where we have set $N=1$ and defined the expansion rates for the two metrics as

H\equiv\frac{\dot{a}}{a}\,,\qquad H_{f}\equiv\frac{1}{\tilde{c}\,\xi}\,\left(H+\frac{\dot{\xi}}{\xi}\right)\,.

(2.12)

Additionally, the contracted Bianchi identity for either metric gives rise to the constraint:

\Gamma\,(H-H_{f}\xi)=0\,.

(2.13)

This equation has two solutions. The first one, $\Gamma=0$ , is known to lead to a non-linear ghost instability [25, 26], a result that holds true also in theories with multiple spin-2 fields [27, 24]. We will then restrict our analysis to the case where

H=H_{f}\,\xi\,.

(2.14)

The constraint implies two consistency relations. The first arises from requiring that (2.8) and (2.7) are consistent under (2.14):

\rho_{m,g}-\frac{M_{g}^{2}\xi^{2}}{M_{f}^{2}}\,\rho_{m,f}=-\frac{\rho}{m^{2}M^{2}}\,,

(2.15)

which can be seen as the equation that determines $\xi$ . The second one arises from the consistency of (2.10) with (2.9),

2(\tilde{c}-1)W=\frac{P+\rho}{M_{g}^{2}}\,,

(2.16)

where we defined

W\equiv\frac{m^{2}(1+\kappa\,\xi^{2})\,\Gamma}{2(1+\kappa)\xi^{2}}-H^{2}\,,

(2.17)

and $\kappa\equiv M_{f}^{2}/M_{g}^{2}$ . The function $W$ is the Higuchi function, since $W>0$ is the dRGT generalisation of the Higuchi bound [28, 29, 24]. The consistency equation (2.16) is an algebraic equation that determines $\tilde{c}$ . Using the background equations, we can write the solution for $\tilde{c}$ in terms of the equation of state of matter $w\equiv P/\rho$ and $\xi$ as

\tilde{c}-1=(1+w)\,\frac{3\,\xi^{2}(\xi^{2}\rho_{m,f}-\kappa\,\rho_{m,g})}{3\,(1+\kappa\,\,\xi^{2})\Gamma-2\,\xi^{4}\rho_{m,f}}\,.

(2.18)

Finally, using (2.12) and (2.14), we can express the time variation of $\xi$ in terms of these functions via

\frac{\dot{\xi}}{H\,\xi}=\tilde{c}-1=(1+w)\,\frac{3\,\xi^{2}(\xi^{2}\rho_{m,f}-\kappa\,\rho_{m,g})}{3\,(1+\kappa\,\xi^{2})\Gamma-2\,\xi^{4}\rho_{m,f}}\,.

(2.19)

At this point, we can in principle replace $\xi$ and $\tilde{c}$ by solving Eqs.(2.15) and (2.16) to remove all references to the metric $f$ . However, (2.15) is a fourth order algebraic equation for $\xi$ and different solutions have different implications for the cosmological evolution. The requirement of a sensible background evolution which is stable against small perturbations is known to severely limit the cosmological solution in the late Universe [30, 27, 24]. In the following section, we revisit these requirements by focussing on an inflaton dominated universe.

2.2 Perturbative stability

For a given dynamical degree of freedom, the free Lagrangian in Fourier space, quadratic in the field, is formally

\mathcal{L}_{k}\ni\mathcal{K}_{\psi}\,\left(|\dot{\psi}_{k}|^{2}-(c_{\psi}^{2}\,k^{2}+m_{\psi}^{2})|\psi|^{2}\right)\,.

(2.20)

For small perturbations, there are in general three stability requirements, which protect the background from ghost, gradient and tachyonic instabilities. In (2.20), these correspond respectively to a positive kinetic term $\mathcal{K}_{\psi}>0$ , a real propagation speed $c_{\psi}^{2}>0$ and a real mass $m_{\psi}^{2}>0$ .

In a cosmology with multiple spin-2 fields, the dynamical metric $f_{\mu\nu}$ brings in a new tensor perturbation adding two tensorial degrees of freedom. However, due to the two–metric interaction term breaking the two diffeomorphism symmetries down to one, a vector (i.e. two additional degrees of freedom) and a scalar mode arise in the gravitational sector. Moreover, we are including a scalar field in the matter sector. As a result, the dynamical and independent degrees of freedom comprise 4 tensor, 2 vector and 2 scalar perturbations. Each of these degrees of freedom need to pass the perturbative stability test around the cosmological background.

One challenge is that both the tensor and scalar sectors in our setup are coupled systems with time-dependent coefficients. In particular, it is not possible to exactly diagonalise the Lagrangian for the scalar perturbations to obtain the mass eigenstates, although the kinetic terms and propagation speeds can be determined unambiguously. Fortunately, the instabilities generated by the imaginary part of the eigenmasses would be mild in our context, therefore we will disregard the conditions arising from avoiding tachyonic instabilities ²²2The mass of gravitational perturbations are typically of order $H$ or smaller. Therefore these instabilities would take longer than the age of the universe to develop.. For the remaining instabilities, we adapt the conditions to avoid ghost and gradient instabilities in all sectors that have been obtained previously in Ref.[24] ³³3The main difference with Ref.[24] arises from the normalisation of the interaction term, which can be accommodated by rescaling $m^{2}\to m^{2}\kappa/(1+\kappa)$ when changing notations. .

For the tensor perturbations, only the kinetic term for the $f$ –metric perturbations brings a stability condition (see Eq.(A.14))

\mathcal{K}_{\rm tensor}\propto\tilde{c}>0\,.

(2.21)

For the vector modes, positivity of the kinetic term in the UV and reality of propagation speed impose [24]

\mathcal{K}_{\rm vector}\propto\frac{\Gamma}{\tilde{c}+1}>0\,,\qquad c^{2}_{\rm vector}\equiv 1+\frac{P+\rho}{4\,M_{g}^{2}W}+\frac{Q\,(P+\rho)}{2\,M_{g}^{2}W\,\Gamma}\,\left(1+\frac{P+\rho}{4\,M_{g}^{2}W}\right)>0\,.

(2.22)

Finally for the scalar modes, avoiding a ghost in the UV demands that [24]

\mathcal{K}_{\rm matter}\propto\frac{\rho+P}{c_{s}^{2}}>0\,,\qquad\mathcal{K}_{\rm scalar}\propto\Gamma\,W>0\,,

(2.23)

where $c_{s}^{2}\equiv dP/d\rho$ for a scalar field is defined in (A.9). For stable gradients, the reality of propagation speeds leads to [24]

c_{\rm matter}^{2}=c_{s}^{2}>0\,,\qquad c_{\rm scalar}^{2}\equiv 1+\frac{2\,Q\,(P+\rho)}{3\,M_{g}^{2}W\,\Gamma}+\frac{H^{2}(2\,Q-\Gamma)(P+\rho)}{6\,M_{g}^{2}W^{2}\Gamma}>0\,.

(2.24)

Thus, throughout the evolution, we need to make sure that the conditions (2.21)–(2.24) are not violated. In the present study, the matter sector consists of a scalar field with a canonical kinetic term, rolling down a positive potential. As a result $c_{s}^{2}=1$ and $\rho+P>0$ are automatically satisfied, while the equation of state is bound to the range $-1\leq w\leq 1$ by construction. In this case, Eq.(2.16) implies that as long as $W>0$ is satisfied, we have $\tilde{c}>1$ , automatically fulfilling the requirement for the stability of tensor perturbations (2.21). Moreover, the propagation speed for the vector mode can be written as

c_{\rm vector}^{2}=\frac{W}{2}\,\left(1+\frac{P+\rho}{4\,M_{g}^{2}W}\right)\left[\frac{1+3\,c_{\rm scalar}^{2}}{H^{2}+2\,W}+\frac{2\,H^{2}}{W\,(H^{2}+2\,W)}\,\left(1+\frac{P+\rho}{4\,M_{g}^{2}W}\right)\right]\,,

(2.25)

which is manifestly positive if the conditions $W>0$ , $P+\rho>0$ and $c_{\rm scalar}^{2}>0$ are satisfied.

In summary, for a matter sector made up of a canonical scalar with a positive potential, the independent stability conditions that are required to be true throughout the evolution are

\Gamma>0\,,\quad W>0\,,\quad c_{\rm scalar}^{2}>0\,.

(2.26)

2.3 Branching of the evolution

The equation for $\xi$ (2.15), in terms of the $\beta_{i}$ , can be written as

\frac{\beta_{1}}{\kappa}\,\frac{1}{\xi}+\left(-\beta_{0}+\frac{3\,\beta_{2}}{\kappa}\right)+3\,\left(-\beta_{1}+\frac{\beta_{3}}{\kappa}\right)\,\xi+\left(-3\,\beta_{2}+\frac{\beta_{4}}{\kappa}\right)\,\xi^{2}-\beta_{3}\xi^{3}=\frac{1+\kappa}{\kappa}\,\frac{\rho}{m^{2}M_{g}^{2}}\,.

(2.27)

In general, this equation has four roots. We find it useful to discuss some properties of the $\xi$ solution in two distinct asymptotic limits that depend on the ratio $\rho/(m^{2}M_{g}^{2})$ . In the late time asymptotics, we are in the low energy regime, where $\rho\ll m^{2}M_{g}^{2}$ and $\xi$ tends to a constant value $\xi\to\xi_{c}$ . Conversely, in the early time asymptotics, we are in the high energy regime such that $\rho\gg m^{2}M_{g}^{2}$ .

As first shown in Ref.[30], the stability of the background solution in the high energy regime is highly sensitive to the specific evolution of $\xi$ as well as the equation of state of the dominant matter fluid. In this section, we catalogue the various evolution branches and determine how the stability conditions manifest themselves in each case. Although the early behaviour poses a roadblock to many dark energy applications [27, 31, 32, 33], we show that there is one evolution branch that leads to a healthy solution in the context of inflation. We then review the late time attractor solution where we can interpret the second metric as an external spin–2 field on top of standard inflationary universe.

2.3.1 Early times (high energy)

The different branches arise from the early asymptotics $\rho\gg m^{2}M_{g}^{2}$ . As we go back in time, the right hand side of (2.27) turns on, thus the left hand side should grow accordingly. We then define two main branches in this case:

•

Branch A: $\lim_{a\to 0}\xi=0$ , i.e. $\xi(t)<\xi_{c}$ with $\dot{\xi}>0$ ;
•

Branch B: $\lim_{a\to 0}\xi=\infty$ , i.e. $\xi(t)>\xi_{c}$ with $\dot{\xi}<0$ .

At this point, we already have enough information about Branch B to determine its stability properties. Using Eq.(2.16), we replace $\tilde{c}$ in Eq.(2.19) to obtain

\dot{\xi}=\frac{H\,\xi\,(P+\rho)}{2\,W\,M_{g}^{2}}\,.

(2.28)

For the matter source that satisfies the null energy condition, which is the case for the scalar inflaton in our construction, we see that $\dot{\xi}<0$ implies $W<0$ . In other words, Branch B suffers from the Higuchi instability. We therefore focus only on Branch A for the discussion of the high energy limit.

In Branch A, $\xi^{-1}$ term dominates the Eq.(2.27) at early times. Keeping up to next-to-leading order terms, we have

\xi^{A}=\frac{\beta_{1}}{1+\kappa}\,\frac{m^{2}M_{g}^{2}}{\rho}\left[1+\frac{m^{2}M_{g}^{2}(3\,\beta_{2}-\kappa\,\beta_{0})}{(1+\kappa)\rho}+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)^{2}\right]\,,

(2.29)

where the positivity of $\xi$ requires

\beta_{1}>0\,.

(2.30)

Thus, for negative $\beta_{1}$ , this branch would not give a physical solution and only the unstable Branch B solutions would be available.

The remaining equations of motion, when expanded for large density, yield

$\displaystyle\tilde{c}-1=$	$\displaystyle 3\,(1+w)\left[1+\frac{m^{2}M_{g}^{2}(3\,\beta_{2}-\kappa\,\beta_{0})}{(1+\kappa)\rho}+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)^{2}\right]\,,$
$\displaystyle 3\,H^{2}=$	$\displaystyle\frac{\rho}{M_{g}^{2}}\left[1+\frac{m^{2}M_{g}^{2}\kappa\,\beta_{0}}{(1+\kappa)\rho}+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)^{2}\right]\,,$
$\displaystyle 2\,\dot{H}=$	$\displaystyle-\frac{\rho+P}{M_{g}^{2}}\left[1-\frac{3\,m^{4}\,M_{g}^{4}\kappa\,\beta_{1}^{2}}{(1+\kappa)^{2}\rho^{2}}+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)^{3}\right]\,.$	(2.31)

In this limit, the key quantities that appear in the stability conditions (2.26) become

$\displaystyle\Gamma=$	$\displaystyle\frac{m^{2}\,M_{g}^{2}\,\beta_{1}^{2}}{(\kappa+1)\,\rho}+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)^{2}\,,$
$\displaystyle W=$	$\displaystyle\frac{\rho}{6\,M_{g}^{2}}+m^{2}\mathcal{O}(1)\,,$
$\displaystyle c^{2}_{\rm scalar}=$	$\displaystyle-(1+2\,w)+\mathcal{O}\left(\frac{m^{2}M_{g}^{2}}{\rho}\right)\,.$	(2.32)

The condition $W>0$ is automatically satisfied for this branch, as we predicted in Eq.(2.28). Combined with the matter stability condition $|w|<1$ , the range for the equation of state that keeps $c^{2}_{\rm scalar}$ positive is

-1<w<-\frac{1}{2}\,.

(2.33)

This branch has sensible Friedmann equations and has no early-time ghost instability. Ref.[27] found that this branch contains a gradient instability in the dark energy context, since for radiation and matter dominated universe, one has $c_{\rm scalar}^{2}<0$ at early times.⁴⁴4This instability may be tamed by non-linear terms [34, 35, 36]. This conclusion of course applies only to fluids with $w>-1/2$ . However, in the context of inflation, we can assume that the beginning of inflation is extended to the energy scales that correspond to the cut-off of the effective theory. Provided that the transition to the low energy regime occurs while inflation is in progress, i.e. $w\sim-1$ , this instability is never invoked.

2.3.2 Late times (low energy)

At late times, the right hand side of (2.27) decreases and the solution for $\xi$ converges to a constant $\xi=\xi_{c}$ defined by

\frac{\beta_{1}}{\kappa}\,\frac{1}{\xi}+\left(-\beta_{0}+\frac{3\,\beta_{2}}{\kappa}\right)+3\,\left(-\beta_{1}+\frac{\beta_{3}}{\kappa}\right)\,\xi+\left(-3\,\beta_{2}+\frac{\beta_{4}}{\kappa}\right)\,\xi^{2}-\beta_{3}\xi^{3}\Bigg{|}_{\xi=\xi_{c}}=0\,.

(2.34)

Although there are four different end-points for the $\xi$ evolution, the late time approach to these values is universal. Defining

\rho_{m,g~{}c}=\rho_{m,g}|_{\xi=\xi_{c}}\,,\qquad\rho_{m,f~{}c}=\rho_{m,f}|_{\xi=\xi_{c}}\,,\qquad\Gamma_{c}=\Gamma|_{\xi=\xi_{c}}\,,\qquad Q_{c}=Q|_{\xi=\xi_{c}}\,,

(2.35)

and expanding for $\rho\ll m^{2}M_{g}^{2}$

	$\displaystyle\xi=$	$\displaystyle\xi_{c}\left[1-\left(\frac{(1+\kappa)\xi_{c}^{2}}{3\,(1+\kappa\,\xi_{c}^{2})\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\,\rho_{m,g~{}c}}\right)\frac{\rho}{m^{2}M_{g}^{2}}+\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)^{2}\right]\,,$
	$\displaystyle\tilde{c}=$	$\displaystyle 1+\frac{\rho+P}{m^{2}M_{g}^{2}}\left[\frac{3(1+\kappa)\xi_{c}^{2}}{3\,(1+\kappa\,X_{\xi}^{2})\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\,\rho_{m,g~{}c}}+\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)\right]\,.$		(2.36)

The equations of motion for the $g$ –metric become

	$\displaystyle 3\,H^{2}=$	$\displaystyle\frac{m^{2}\kappa\,\rho_{m,g~{}c}}{1+\kappa}+\frac{(3\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\rho_{m,g~{}c})}{3\,(1+\kappa\,\xi_{c}^{2})\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\,\rho_{m,g~{}c}}\,\frac{\rho}{M_{g}^{2}}+m^{2}\,\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)^{2}\,,$
	$\displaystyle 2\,\dot{H}=$	$\displaystyle-\frac{(3\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\rho_{m,g~{}c})}{3\,(1+\kappa\,\xi_{c}^{2})\,\Gamma_{c}-2\,\kappa\,\xi_{c}^{2}\,\rho_{m,g~{}c}}\,\frac{\rho+P}{M_{g}^{2}}+m^{2}\,\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)^{2}\,.$		(2.37)

At late times, some of the functions in the stability conditions (2.26) are

	$\displaystyle W=$	$\displaystyle\frac{m^{2}}{1+\kappa}\left(\frac{\Gamma_{c}(1+\kappa\,\xi_{c}^{2})}{2\,\xi_{c}^{2}}-\frac{\kappa\,\rho_{m,g~{}c}}{3}\right)+m^{2}\,\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)\,,$
	$\displaystyle c^{2}_{{\rm scalar}}=$	$\displaystyle 1+\mathcal{O}\left(\frac{\rho}{m^{2}M_{g}^{2}}\right)\,.$		(2.38)

We see that the gradient instability in the scalar sector never appears, while the no-ghost conditions $W>0$ and $\Gamma>0$ are parameter dependent and can be satisfied by appropriate choices of $\Gamma_{c}$ and $\rho_{m,g~{}c}$ .

In the remainder of the paper, we consider a scenario where the transition from the high to low energy regime occurs during inflation, while the largest observable scale is deep inside the horizon. In the calculations of the correlators, we will therefore take

\xi\simeq\xi_{c}\,,\qquad\tilde{c}\simeq 1\,.

(2.39)

The second of these is automatic, while the specific value of $\xi_{c}$ can be chosen by fixing one of the $\beta_{n}$ parameters.

All late-time instabilities are avoided if

\Gamma_{c}>0\,,\qquad\frac{\Gamma_{c}(1+\kappa\,\xi_{c}^{2})}{2\,\xi_{c}^{2}}-\frac{\kappa\,\rho_{m,g~{}c}}{3}>0\,.

(2.40)

It is important to note that in this low energy regime with $\rho\ll m^{2}M_{g}^{2}$ , the energy contribution from the 2–metric interaction term can be parametrically larger than the inflaton density, which may lead to de Sitter evolution to that goes on indefinitely. To avoid this, we set to zero the contribution to Friedmann equation:

\rho_{m,g~{}c}=0\,.

(2.41)

Evaluating Eq.(2.15) in the late time asymptotics, the above selection implies that we also need to fix $\rho_{m,f~{}c}=0$ . Both of these can be realised at the price of fixing two of the $\beta_{n}$ parameters. As a result of this choice, the background equations at late times (2.37) become:

	$\displaystyle 3\,H^{2}\simeq$	$\displaystyle\frac{\rho}{M_{eff}^{2}}\,,$
	$\displaystyle 2\,\dot{H}\simeq$	$\displaystyle-\frac{\rho+P}{M_{eff}^{2}}\,.$		(2.42)

Thus at the background level, the two metrics are decoupled and the resulting cosmology is simply the usual inflationary universe in General Relativity, with an effective Planck mass defined as

M_{eff}^{2}\equiv(1+\kappa\,\xi_{c}^{2})\,M_{g}^{2}\,.

(2.43)

This construction allows us to interpret the second metric as an external spin–2 field. In this case, all of the stability conditions in Eq. (2.26) are satisfied at late times provided that $\Gamma_{c}>0$ .

2.4 Action for tensors in the low energy regime

In this Section, we present the action up to cubic order in the low energy regime discussed in Sec. 2.3.2. Excluding the background portion, we decompose the action as

S=\int d^{4}x\,a^{3}\left(\mathcal{L}_{\rm tot}^{(2)}+\mathcal{L}_{\rm tot}^{(3)}+\dots\right)\,.

(2.44)

The full expressions are given in Appendix A. Here, we present them in the low energy limit by fixing $\xi=\xi_{c}$ and $\tilde{c}=1$ . From here on, we suppress the subscript $c$ that represents the constant values reached at late times.

At quadratic order in tensor perturbations, the dynamics of the free tensor modes are described by (A.14)

	$\displaystyle{\mathcal{L}}^{(2)}=$	$\displaystyle\frac{M_{g}^{2}}{8}\,\Bigg{[}\dot{h}_{ij}\dot{h}^{ij}-\frac{1}{a^{2}}\partial_{k}h^{ij}\partial^{k}h_{ij}+\kappa\,\xi^{2}\,\left(\dot{\gamma}_{ij}\dot{\gamma}^{ij}-\frac{1}{a^{2}}\partial_{k}\gamma^{ij}\partial^{k}\gamma_{ij}\right)$
		$\displaystyle\qquad\qquad\qquad\qquad-\frac{m^{2}\kappa\,\Gamma}{1+\kappa}\,(h_{ij}h^{ij}-2\,h^{ij}\gamma_{ij}+\gamma^{ij}\gamma_{ij})\Bigg{]}\,,$		(2.45)

where we defined $\kappa\equiv M_{f}^{2}/M_{g}^{2}$ . We see that there is a constant interaction term proportional to $\Gamma$ .
The cubic interactions that involve the tensor modes are (A.15)

$\displaystyle{\mathcal{L}}^{(3)}_{\rm tot}=$	$\displaystyle-\frac{M_{g}^{2}}{4}h_{ij}\left[\dot{h}^{ik}\dot{h}_{k}^{\;\;j}-\frac{1}{2\,a^{2}}\Big{(}\partial^{i}h^{kl}\partial^{j}h{}_{kl}-2\,\partial^{i}h_{kl}\partial^{l}h^{jk}+2\,\partial_{l}h^{i}_{\;\;k}\,\partial^{l}h^{jk}\Big{)}\right]$
	$\displaystyle-\frac{M_{f}^{2}}{4}\xi^{2}\,\gamma_{ij}\left[\dot{\gamma}^{ik}\dot{\gamma}_{k}^{\;\;j}-\frac{1}{2\,a^{2}}\Big{(}\partial^{i}\gamma^{kl}\partial^{j}\gamma{}_{kl}-2\,\partial^{i}\gamma_{kl}\partial^{l}\gamma^{jk}+2\,\partial_{l}\gamma^{i}_{\;\;k}\,\partial^{l}\gamma^{jk}\Big{)}\right]$
	$\displaystyle+\frac{m^{2}M^{2}}{48}\left(\Gamma-2\,Q\right)[(h-\gamma)^{3}]+\frac{m^{2}M^{2}}{8}\Gamma\,[(h-\gamma)^{2}(h+\gamma)]\,.$	(2.46)

The non-derivative cubic interactions are controlled by combinations of $Q$ and $\Gamma$ .

In Sec.3, we will be using these expressions to calculate the two-point and three-point correlation functions.

3 Correlators amplitude and shape

We now proceed with the computation of the power spectrum and bispectrum for $h$ , given by:

	$\displaystyle\langle h^{\lambda_{1}}_{\textbf{k}_{1}}h^{\lambda_{2}}_{\textbf{k}_{2}}\rangle=(2\pi)^{3}\delta^{(3)}(\textbf{k}_{1}+\textbf{k}_{2})\delta_{\lambda_{1}\lambda_{2}}P_{h}(k)\,,$		(3.1)
	$\displaystyle\langle h^{\lambda_{1}}_{\textbf{k}_{1}}h^{\lambda_{2}}_{\textbf{k}_{2}}h^{\lambda_{3}}_{\textbf{k}_{3}}\rangle=(2\pi)^{3}\delta^{(3)}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})B_{h}^{\lambda_{1}\lambda_{2}\lambda_{3}}(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3})\,.$		(3.2)

3.1 Power spectrum

The starting point for the power spectrum computation is the Lagrangian given in Eq. (2.4). From here one obtains the equations of motion for the canonically normalised fields, $h_{c}\equiv aM_{g}h$ and $\gamma_{c}\equiv aM_{f}\xi\gamma$ ,

	$\displaystyle h^{{}^{\prime\prime}}_{c}+\left(k^{2}+a^{2}m_{h}^{2}-\frac{a^{{}^{\prime\prime}}}{a}\right)h_{c}=a^{2}m_{X}^{2}\gamma_{c}\,,$		(3.3)
	$\displaystyle\gamma^{{}^{\prime\prime}}_{c}+\left(k^{2}+a^{2}m_{\gamma}^{2}-\frac{z^{{}^{\prime\prime}}}{z}\right)\gamma_{c}=a^{2}m_{X}^{2}h_{c}\,.$		(3.4)

where $s\equiv aM_{g}$ and $z\equiv aM_{f}\xi$ . Using Eq. (2.3) and the definition of $\kappa$ , the masses of the two modes are found to be equal to

m_{h}^{2}\equiv\frac{m^{2}\kappa\,\Gamma}{\kappa+1}\,,\qquad m_{\gamma}^{2}\equiv\frac{m_{h}^{2}}{\kappa\,\xi^{2}}\,,

(3.5)

whereas the “mass” scale for the cross term is given by

m_{X}^{2}\equiv\frac{m_{h}^{2}}{\sqrt{\kappa}\,\xi}\,.

(3.6)

The mixing term can be treated as a perturbation on top of the free Lagrangian as long as the following condition is in place

\frac{m^{2}\,\sqrt{\kappa}\,\Gamma}{H^{2}(1+\kappa)\,\xi}\ll 1

(3.7)

or, equivalently, $m_{h}^{2}/H^{2}\ll\sqrt{\kappa}\xi$ . From now on we will assume condition (3.7), as this allows us to use the standard mode functions for free fields and introduce the small interaction via the in-in formalism. In addition, we must take the Higuchi bound into account. From the definition in (2.17), the latter reads

\frac{m^{2}(1+\kappa\,\xi^{2})\Gamma}{(1+\kappa)\,\xi^{2}}>2\,H^{2}

(3.8)

or, equivalently, as $m^{2}_{h}/H^{2}>2\kappa\xi^{2}/(1+\kappa\xi^{2})$ , which is compatible with Eq.(3.7) for

\sqrt{\kappa}\,\xi\ll 1\,.

(3.9)

In this regime, the three scales $m_{h}^{2}$ , $m_{X}^{2}$ and $m_{\gamma}^{2}$ are organised according to the following hierarchy

m_{h}^{2}\ll m_{X}^{2}\ll m_{\gamma}^{2}\,.

(3.10)

Notice also that, for $\sqrt{\kappa}\xi\ll 1$ , Eq. (3.7) ensures that $h$ is fairly light. As one would expect (see also our explicit result in Eq. (3.14)), this is crucial for the mode that couples to matter fields not to be suppressed at late times.

As we will explicitly verify below, the leading order contribution to $P_{h}$ is given by the diagram on the left-hand side of Fig. 1. This is expected given that the 2-vertex interaction is taken to be small by construction.

Refer to caption — Figure 1: Left: the standard tensor two-point function from vacuum fluctuations. Right: the diagram corresponding to the power spectrum being “sourced” by the $h-\gamma$ mixing.

Let us first derive the mode-function for $h$ from its free quadratic Lagrangian. Working at next-to-leading order in slow roll, the equation of motion reads

\displaystyle h_{k}^{{}^{\prime\prime}}-2\left(\frac{1+\epsilon_{*}}{\tau}\right)h_{k}^{{}^{\prime}}+\left(k^{2}+\frac{m_{h}^{2}}{H^{2}_{*}\tau^{2}}\right)h_{k}=0\,,

(3.11)

where the suffix “*” denotes a quantity calculated at the reference time $\tau_{*}$ ⁵⁵5Such reference point is typically chosen to be the (conformal) time at which the longest mode that is observationally relevant today exits the horizon.. The solution for $h$ reads

h_{k}(\tau)=(-\tau)^{3/2+\epsilon_{*}}\cdot\,\mathcal{A}\cdot\,H^{(1)}_{\nu_{T}}(-k\tau)\,,

(3.12)

where $\nu_{T}\equiv\sqrt{9/4+3\epsilon-m_{h}^{2}/H_{*}^{2}}\simeq 3/2+\epsilon-m_{h}^{2}/(3H_{*}^{2})$ and

\mathcal{A}\equiv\frac{H_{*}(1-\epsilon_{*})}{M_{g}}\sqrt{\pi}(-\tau_{*})^{-\epsilon_{*}}e^{i\left(\nu_{T}\cdot\pi/2+\pi/4\right)}\,.

(3.13)

The mode-function in Eq. (3.12) may be expanded at late times to give

\displaystyle h_{k}(\tau)|_{k\tau\rightarrow 0}

\displaystyle=\frac{\sqrt{2}\,H_{*}}{M_{g}k^{3/2}}(-i)\cdot e^{i\left(\nu_{T}\cdot\pi/2+\pi/4\right)}\cdot(-k\tau_{*})^{3/2-\nu_{T}}\cdot\left(\frac{\tau}{\tau_{*}}\right)^{\frac{m_{h}^{2}}{3H^{2}_{*}}}\,.

(3.14)

The corresponding contribution to the power spectrum then amounts to

P_{h}^{(0)}=\frac{2\,H^{2}_{*}}{M_{g}^{2}k^{3}}\cdot(-k\tau_{*})^{3-2\nu_{T}}\cdot\left(\frac{\tau}{\tau_{*}}\right)^{\frac{2m_{h}^{2}}{3H^{2}_{*}}}\,.

(3.15)

In this expression, $\tau/\tau_{*}$ corresponds to the number of e-folds, $N_{*}$ , between the horizon exit of the mode $k_{*}\equiv-1/\tau_{*}$ and the end of inflation.

Let us now compute the other contribution to the tree-level two-point function (depicted in the right panel of Fig. 1), which we will indicate as $P_{h}^{\text{mix}}$ ,

\displaystyle P_{h}^{\text{mix}}=A^{\text{mix}}+B^{\text{mix}}\,,

(3.16)

where

	$\displaystyle A^{\text{mix}}\equiv\frac{\left(m^{2}M^{2}\Gamma\right)^{2}}{32}h_{k_{1}}(\tau)h_{k_{1}}^{*}(\tau)\Big{\\|}\int_{-\infty}^{\tau}d\tau_{1}\,a^{4}(\tau_{1})h_{k_{1}}(\tau_{1})\gamma_{k_{1}}(\tau_{1})\Big{\\|}^{2}\,,$		(3.17)
	$\displaystyle B^{\text{mix}}\equiv-\frac{\left(m^{2}M^{2}\Gamma\right)^{2}}{16}\text{Re}\left[h_{k_{1}}^{2}(\tau)\int_{-\infty}^{\tau}d\tau_{1}\,a^{4}(\tau_{1})h^{}_{k_{1}}(\tau_{1})\gamma_{k_{1}}(\tau_{1})\int_{-\infty}^{\tau_{1}}d\tau_{2}\,a^{4}(\tau_{2})h_{k_{1}}^{}(\tau_{2})\gamma_{k_{1}}^{*}(\tau_{2})\right]\,.$

The mode function of $\gamma$ is obtained from the free-field equation of motion derived from Eq. (2.4) and reads, this time to leading order in slow-roll,

\gamma_{k}(\tau)=(-\tau)^{3/2}\cdot\,\frac{H}{M_{g}\sqrt{\kappa}\xi}\sqrt{\pi}\,e^{i\left(\tilde{\nu_{T}}\cdot\pi/2+\pi/4\right)}\cdot\,H^{(1)}_{\tilde{\nu}_{T}}(-k\tau)\,,

(3.19)

where $\tilde{\nu}_{T}\equiv\sqrt{9/4-m_{\gamma}^{2}/H^{2}}$ , and we assumed $m_{\gamma}\leq 3H/2$ . After a suitable change of variables ( $x\equiv-k_{1}\tau$ ) and more algebraic manipulations, one finds

	$\displaystyle A^{\text{mix}}=P_{h}^{(0)}\cdot\frac{\pi^{2}}{32}\left(\frac{m^{2}M^{2}\Gamma}{H^{2}M_{g}^{2}}\right)^{2}\cdot\frac{1}{\kappa\xi^{2}}\\|\mathcal{I}_{A}\\|^{2}\,,$		(3.20)
	$\displaystyle B^{\text{mix}}=P_{h}^{(0)}\cdot\frac{\pi^{2}}{16}\left(\frac{m^{2}M^{2}\Gamma}{H^{2}M_{g}^{2}}\right)^{2}\cdot\frac{1}{\kappa\xi^{2}}\text{Re}\left[\mathcal{I}_{B}\right]\,,$		(3.21)

where

			$\displaystyle\mathcal{I}_{A}\equiv\int_{x}^{+\infty}\frac{dx_{1}}{x_{1}}\,H_{\nu_{T}}^{(1)}(x_{1})H_{\tilde{\nu}_{T}}^{(1)}(x_{1})\,,$		(3.22)
			$\displaystyle\mathcal{I}_{B}\equiv\int_{x}^{+\infty}\frac{dx_{1}}{x_{1}}\,H_{\nu_{T}}^{(1)}(x_{1})H_{\tilde{\nu}_{T}}^{(1)}(x_{1})\int_{x_{1}}^{+\infty}\frac{dx_{2}}{x_{2}}\,H_{\nu_{T}}^{(1)}(x_{2})H_{\tilde{\nu}_{T}}^{(1)*}(x_{2})\,.$		(3.23)

The total (i.e. accounting for both polarisations) tree-level power spectrum, $P_{\rm tot}\equiv 2\left(P_{h}^{(0)}+P_{h}^{\text{mix}}\right)$ , reads

	$\displaystyle P_{\rm tot}$	$\displaystyle=$	$\displaystyle 2P_{h}^{(0)}\left[1+\left(\frac{m^{2}M^{2}\Gamma}{H^{2}M_{g}^{2}}\right)^{2}\cdot\frac{1}{\kappa\xi^{2}}\frac{\pi^{2}}{32}\left(\\|\mathcal{I}_{A}\\|^{2}+2\text{Re}[\mathcal{I}_{B}]\right)\right]$		(3.24)
			$\displaystyle\equiv 2P_{h}^{(0)}\left[1+\Theta\cdot\mathcal{C}(\nu_{T},\tilde{\nu}_{T})\right]\,,$		(3.24)

where $\mathcal{C}(\nu_{T},\tilde{\nu}_{T})\equiv\|\mathcal{I}_{A}\|^{2}+2\text{Re}[\mathcal{I}_{B}]$ and, using Eqs. (3.6)-(3.7), one has $\Theta\equiv\left(\frac{\pi\,m_{X}^{2}}{4\sqrt{2}\,H^{2}}\right)^{2}\ll 1$ in our chosen regime. The function $\mathcal{C}$ takes on small values (see Fig. 2), hence $P_{h}^{\text{mix}}\ll P_{h}^{(0)}$ . The spectral index for $P_{\text{tot}}$ is given, to leading order, by the one obtained in Eq. (3.14), as we verify below:

\displaystyle n_{T}=\frac{d\ln\mathcal{P}_{\text{tot}}}{d\ln k}=\frac{d\ln\left[2\mathcal{P}_{h}^{(0)}\left(1+\Theta\cdot\mathcal{C}\right)\right]}{d\ln k}=\frac{d\ln\left[\mathcal{P}_{h}^{(0)}\right]}{d\ln k}+\frac{d\ln\left[1+\Theta\cdot\mathcal{C}\right]}{d\ln k}\,,

(3.25)

where we have introduced the reduced power spectra defined as $\mathcal{P}\equiv P\cdot k^{3}/(2\pi^{2})$ . One finds

\displaystyle\frac{d\ln\left[1+\Theta\cdot\mathcal{C}\right]}{d\ln k}\simeq 4\epsilon\,\Theta\cdot\mathcal{C}\,,

(3.26)

which is negligible compared to the contribution from $\mathcal{P}_{h}^{(0)}$ , leading again to $n_{T}\simeq 3-2\nu_{T}$ . Notice that the latter can also be rewritten as $n_{T}\simeq\frac{2m_{h}^{2}}{3H_{*}^{2}}-2\epsilon$ , hence the tensor power spectrum is blue-(red)-tilted depending whether the mass term is larger(smaller) than the slow-roll parameter. A close inspection of our results so far reveals that as $m_{h}$ varies from smaller to larger values (and viceversa for $\epsilon$ ), the amplitude of the tensor-to-scalar ratio, $r\equiv\mathcal{P}_{\text{tot}}/\mathcal{P}_{\zeta}$ (with $\mathcal{P}_{\zeta}=H^{2}/(8\pi^{2}\epsilon M_{g}^{2})$ ), decreases whereas $n_{T}$ increases:

r_{*}=16\,\epsilon_{*}\left(\frac{\tau}{\tau_{*}}\right)^{\frac{2m_{h}^{2}}{3H^{2}_{*}}}\simeq 16\,\epsilon_{*}\,e^{-\frac{33\,m_{h}^{2}}{H_{*}^{2}}}\,,\quad\quad n_{T_{*}}\simeq\frac{2m_{h}^{2}}{3H_{*}^{2}}-2\epsilon_{*}\,,

(3.27)

having set $N_{*}\simeq 50$ . As expected, one recovers the results of standard single-field inflation, $r\simeq 16\,\epsilon$ and $n_{T}\simeq-2\epsilon$ , by setting $m_{h}=0$ in the quantities appearing in Eq. (3.27).
Because of the opposite effects of $m_{h}$ (and $\epsilon$ ) on the tensor-to-scalar ratio and on the spectral tilt, the model is able to easily accomodate fairly large values of $r$ while at the same time predicting an $n_{T}$ not too far from zero. This ensures that the tensor power spectrum is potentially detectable at CMB scales by upcoming probes such as the LiteBIRD satellite [4] and the new generation of ground-based experiments [37, 38, 39, 40], while at interferometer scales it could be visible by detectors such as the proposed Big-Bang Observer [41] or DECIGO [42].

3.2 Bispectrum

The total bispectrum of $h$ at tree-level is given by the sum of the four contributions shown in Fig. 3:

B^{\lambda_{1}\lambda_{2}\lambda_{3}}_{h}=\mathcal{A}^{\lambda_{1}\lambda_{2}\lambda_{3}}\left(B_{h}^{[1]}+B_{h}^{[2]}+B_{h}^{[3]}+B_{h}^{[4]}\right)\,,

(3.28)

where $\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}\equiv\epsilon_{ij}^{\lambda_{1}}(\hat{k}_{1})\epsilon_{jk}^{\lambda_{2}}(\hat{k}_{2})\epsilon_{ki}^{\lambda_{3}}(\hat{k}_{3})$ can be expressed as follows

	$\displaystyle\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}$	$\displaystyle=$	$\displaystyle-\frac{1}{64\,x_{2}^{2}\,x_{3}^{2}}\left[2x_{2}+\left(1+x_{2}^{2}-x_{3}^{2}\right)\lambda_{1}\lambda_{2}\right]\left[-2x_{3}+\left(-1+x_{2}^{2}-x_{3}^{2}\right)\lambda_{1}\lambda_{3}\right]$		(3.29)
			$\displaystyle\quad\times\left[2x_{2}x_{3}+\left(-1+x_{2}^{2}+x_{3}^{2}\right)\lambda_{2}\lambda_{3}\right]\,,$		(3.29)

with $x_{i}\equiv k_{i}/k_{1}$ , and $\lambda_{i}=\pm 1$ respectively for R,L helicities. The terms in the cubic Lagrangian responsible for each of the $B_{h}^{[i]}$ are given in the third line of Eq. (2.46)⁶⁶6The derivative interactions (first line of (2.46)) contribute to the tree-level bispectrum with a diagram as in Fig. 3-[1]. However, these interactions are not new: they are already present in minimal single-field slow-roll models. Further, we may not rely on the flexibility of the theory parameter space (e.g. the ratio $Q/\Gamma$ ) to enhance their contribution., which we write schematically as

\mathcal{L}^{(3)}=\alpha_{hhh}\,h^{3}+\alpha_{hh\gamma}\,h^{2}\gamma+\alpha_{h\gamma\gamma}\,h\gamma^{2}+\alpha_{\gamma\gamma\gamma}\,\gamma^{3}\,,

(3.30)

with

	$\displaystyle\alpha_{hhh}\equiv\frac{m^{2}M^{2}}{48}\left(7\,\Gamma-2\,Q\right)\,,\quad\quad\alpha_{hh\gamma}\equiv-\frac{m^{2}M^{2}}{16}\left(3\,\Gamma-2\,Q\right)\,,$
	$\displaystyle\alpha_{h\gamma\gamma}\equiv-\frac{m^{2}M^{2}}{16}\left(\Gamma+2\,Q\right)\,,\quad\quad\alpha_{\gamma\gamma\gamma}\equiv\frac{m^{2}M^{2}}{48}\left(5\,\Gamma+2\,Q\right)\,,$		(3.31)

where the explicit expressions for $\Gamma$ and $Q$ are in Eqs. (2.6). One can obtain the parametric dependence of the amplitude for the four contributions as follows:

	$\displaystyle f_{\text{NL}}^{[1]}\propto\epsilon^{2}\,\frac{m_{h}^{2}}{H^{2}}\left(7-\frac{2Q}{\Gamma}\right)\,,\quad\quad\quad\quad\quad\quad f_{\text{NL}}^{[2]}\propto\epsilon^{2}\,\frac{m_{h}^{2}}{H^{2}}\left(3-\frac{2Q}{\Gamma}\right)\left(\frac{m_{h}/H}{\sqrt{\kappa}\xi}\right)^{2}\,,$
	$\displaystyle f_{\text{NL}}^{[3]}\propto\epsilon^{2}\,\frac{m_{h}^{2}}{H^{2}}\left(1+\frac{2Q}{\Gamma}\right)\left(\frac{m_{h}/H}{\sqrt{\kappa}\xi}\right)^{4}\,,\quad f_{\text{NL}}^{[4]}\propto\epsilon^{2}\,\frac{m_{h}^{2}}{H^{2}}\left(5+\frac{2Q}{\Gamma}\right)\left(\frac{m_{h}/H}{\sqrt{\kappa}\xi}\right)^{6}\,,$

where $f_{\text{NL}}^{[i]}\propto B^{[i]}_{h}/P_{\zeta}^{2}$ , $P_{\zeta}^{2}$ being the power spectrum of curvature fluctuations (see Eq. (3.43) for the exact definition of $f_{\text{NL}}$ in the specific momentum configuration that will be relevant for our final results). The above expressions for $f_{\text{NL}}^{[2]}$ , $f_{\text{NL}}^{[3]}$ and $f_{\text{NL}}^{[4]}$ have been obtained under the assumption $m_{\gamma}\leq 3H/2$ . In this regime, and under the conditions on $m_{h}/H$ and $\sqrt{\kappa}\xi$ spelled out at the beginning of Sec. 3.1 (see also Sec. 4 for a full recap of these relations), one finds that $m_{h}/H\simeq\mathcal{O}(\sqrt{\kappa}\xi)$ . All the quantities in Eq. (3.2) have therefore similar magnitudes. The complementary regime, $m_{\gamma}>3H/2$ , on the other hand, would result in oscillatory factors in the mode functions which suppress the corresponding amplitudes for $B_{h}^{[2,3,4]}$ (see e.g. [43]).

In the following, we will perform a full derivation of $B_{h}^{[1]}$ , i.e. the contribution from the first diagram in Fig. 3. Our computations above show that it provides a very good indication of the overall level of non-Gaussianity in our model. The tree-level diagram is given by:

\langle h_{\textbf{k}_{1}}^{\lambda_{1}}(\tau)h_{\textbf{k}_{2}}^{\lambda_{2}}(\tau)h_{\textbf{k}_{3}}^{\lambda_{3}}(\tau)\rangle^{[1]}=-i\int_{-\infty}^{\tau}d\tau^{\prime}\,a(\tau^{\prime})\langle\left[h_{\textbf{k}_{1}}^{\lambda_{1}}(\tau)h_{\textbf{k}_{2}}^{\lambda_{2}}(\tau)h_{\textbf{k}_{3}}^{\lambda_{3}}(\tau),H^{(3)}(\tau^{\prime})\right]\rangle\,,

(3.33)

where the brackets indicate a commutator and $H^{(3)}$ is the non-derivative interaction Hamiltonian,

H^{(3)}(\tau^{\prime})=-\alpha_{hhh}\int d^{3}\textbf{x}^{\prime}\,a^{3}(\tau^{\prime})\,h^{3}(\textbf{x}^{\prime},\tau^{\prime})\,.

(3.34)

The latter can be expanded in Fourier space as

	$\displaystyle H^{(3)}(\tau^{\prime})=$		$\displaystyle-\alpha_{hhh}\int d^{3}\textbf{x}^{\prime}\,a^{3}(\tau^{\prime})\,\int\frac{d^{3}k_{1}^{{}^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}k_{2}^{{}^{\prime}}}{(2\pi)^{3}}\int\frac{d^{3}k_{3}^{{}^{\prime}}}{(2\pi)^{3}}e^{-i\textbf{x}^{\prime}\cdot\left(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3}\right)}$
			$\displaystyle\times\,h_{ij,\textbf{k}_{1}}(\tau^{\prime})h_{jk,\textbf{k}_{2}}(\tau^{\prime})h_{ki,\textbf{k}_{3}}(\tau^{\prime})\,,$

where

h_{ij,\textbf{k}_{1}}(\tau^{\prime})=\sum_{\lambda^{\prime}}\epsilon_{ij}^{\lambda^{\prime}}(\hat{k}_{1})\left[a_{\textbf{k}_{1}}^{\lambda^{\prime}}h_{k_{1}}(\tau^{\prime})+a_{-\textbf{k}_{1}}^{\lambda^{\prime}\dagger}h_{k_{1}}^{*}(\tau^{\prime})\right]\,,

(3.35)

with commutation relations $\left[a_{\textbf{k}_{1}^{{}^{\prime}}}^{\lambda^{\prime}},a_{-\textbf{k}_{1}^{{}^{\prime\prime}}}^{\lambda^{\prime\prime}\dagger}\right]=(2\pi)^{3}\delta_{\lambda^{\prime}\lambda^{\prime\prime}}\delta^{(3)}(\textbf{k}_{1}^{{}^{\prime}}+\textbf{k}_{1}^{{}^{\prime\prime}})$ . Similarly, for $\gamma$ one has

\gamma_{ij,\textbf{k}_{1}}(\tau^{\prime})=\sum_{\lambda^{\prime}}\epsilon_{ij}^{\lambda^{\prime}}(\hat{k}_{1})\left[b_{\textbf{k}_{1}}^{\lambda^{\prime}}\gamma_{k_{1}}(\tau^{\prime})+b_{-\textbf{k}_{1}}^{\lambda^{\prime}\dagger}\gamma_{k_{1}}^{*}(\tau^{\prime})\right]\,,

(3.36)

and $\left[b_{\textbf{k}_{1}^{{}^{\prime}}}^{\lambda^{\prime}},b_{-\textbf{k}_{1}^{{}^{\prime\prime}}}^{\lambda^{\prime\prime}\dagger}\right]=(2\pi)^{3}\delta_{\lambda^{\prime}\lambda^{\prime\prime}}\delta^{(3)}(\textbf{k}_{1}^{{}^{\prime}}+\textbf{k}_{1}^{{}^{\prime\prime}})$ . After performing the Wick contractions and a few simplifications, the bispectrum reads

	$\displaystyle\langle h_{\textbf{k}_{1}}^{\lambda_{1}}h_{\textbf{k}_{2}}^{\lambda_{2}}h_{\textbf{k}_{3}}^{\lambda_{3}}\rangle^{[1]}$	$\displaystyle=$	$\displaystyle(2\pi)^{3}\delta^{(3)}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})\,12\,\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}$
			$\displaystyle\cdot\text{Im}\left[-\alpha_{hhh}\,h_{k_{1}}(\tau)h_{k_{2}}(\tau)h_{k_{3}}(\tau)\int_{-\infty}^{\tau}d\tau^{{}^{\prime}}\,a^{4}(\tau^{{}^{\prime}})\cdot h_{k_{1}}^{}(\tau^{{}^{\prime}})h_{k_{2}}^{}(\tau^{{}^{\prime}})h_{k_{3}}^{*}(\tau^{{}^{\prime}})\right]\,,$

where the factor of 6 accounts for the number of permutations. Let us now perform a change of variables and replace the mode functions in Eq. (3.2). We define $x\equiv-k_{1}\tau$ , $x^{\prime}\equiv-k_{1}\tau^{\prime}$ . One then finds

$\displaystyle\langle h_{\textbf{k}_{1}}^{\lambda_{1}}h_{\textbf{k}_{2}}^{\lambda_{2}}h_{\textbf{k}_{3}}^{\lambda_{3}}\rangle^{[1]}$	$\displaystyle=$	$\displaystyle(2\pi)^{3}\delta^{(3)}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})\,48\sqrt{2}\,\pi^{3/2}\,\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}\left(-\frac{\alpha_{hhh}H^{2}_{}}{M_{g}^{6}}\right)\left(\frac{x}{x_{}}\right)^{3\epsilon_{*}}x^{9/2-3\nu_{T}}$	(3.38)
		$\displaystyle\cdot\text{Re}\left[\int_{x}^{+\infty}dx^{\prime}\,\sqrt{x^{\prime}}\,\left(\frac{x^{\prime}}{x_{}}\right)^{-\epsilon_{}}H_{\nu_{T}}^{(1)}\left(\frac{}{}x^{\prime}\right)H_{\nu_{T}}^{(1)}\left(\frac{k_{2}}{k_{1}}x^{\prime}\right)H_{\nu_{T}}^{(1)}\left(\frac{k_{3}}{k_{1}}x^{\prime}\right)\right]$
		$\displaystyle\cdot\left(\frac{k_{2}}{k_{1}}\right)^{3-\nu_{T}}\left(\frac{k_{3}}{k_{1}}\right)^{3-\nu_{T}}\frac{k_{1}^{3}}{(k_{1}k_{2}k_{3})^{3}}\,.$

In the regime where $m_{h}\ll 3H/2$ , the bispectrum can be computed analytically and one obtains

\displaystyle\langle h_{\textbf{k}_{1}}^{\lambda_{1}}h_{\textbf{k}_{2}}^{\lambda_{2}}h_{\textbf{k}_{3}}^{\lambda_{3}}\rangle^{[1]}

\displaystyle=

\displaystyle(2\pi)^{3}\delta^{(3)}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}\frac{H_{*}^{4}}{M_{g}^{4}}\cdot\frac{m_{h}^{2}}{H_{*}^{2}}\cdot\left(7-\frac{2Q}{\Gamma}\right)\cdot\mathcal{F}(k_{1},k_{2},k_{3})\,,

to leading order in both $m_{h}/H$ and the slow-roll parameter, with

\displaystyle\mathcal{F}(k_{1},k_{2},k_{3})\equiv-\frac{2}{3\,k_{1}^{3}k_{2}^{3}k_{3}^{3}}\left[\left(\frac{32}{3}-8\gamma\right)\sum_{i=1}^{3}k_{i}^{3}+8\sum_{i\not=j=1}^{3}k_{i}^{2}k_{j}-8\,\prod_{i=1}^{3}k_{i}-8\,\sum_{i=1}^{3}k_{i}^{3}\cdot\ln\left(\frac{k_{t}\,x}{k_{1}}\right)\right]\,.

The quantity $\gamma$ is here Euler’s constant and $k_{t}\equiv k_{1}+k_{2}+k_{3}$ . One may verify that this bispectrum receives its largest contribution in the configuration where one of the momenta is much smaller than the other two, e.g. $k_{3}\ll k_{1}\simeq k_{2}$ (see Sec. 3.3 for a full discussion on the shape function). In this limit, one finds

\displaystyle B_{h}^{[1]}\simeq\frac{H_{*}^{4}}{M_{g}^{4}}\cdot\frac{m_{h}^{2}}{H_{*}^{2}}\cdot\left(7-\frac{2Q}{\Gamma}\right)\left(\frac{32\,(3\gamma+\ln(8)-7+\ln(x))}{9\,k_{1}^{3}k_{3}^{3}}\right)\simeq-\frac{H_{*}^{4}}{M_{g}^{4}}\cdot\frac{m_{h}^{2}}{H_{*}^{2}}\left(7-\frac{2Q}{\Gamma}\right)\cdot\frac{545}{k_{1}^{3}k_{3}^{3}}

and

\displaystyle\mathcal{A}_{\lambda_{1}\lambda_{2}\lambda_{3}}\simeq\begin{cases}\frac{1}{4}&\,\,\text{for $\,\,\lambda_{1},\,\lambda_{2},\,\lambda_{3}=\text{RRR},\,\text{LLL},\,\text{RRL},\,\text{LLR}$}\\ \frac{1}{16}\,\frac{k_{3}^{2}}{k_{2}^{2}}&\,\,\text{for $\,\,\lambda_{1},\,\lambda_{2},\,\lambda_{3}=\text{LRR},\,\text{RLL},\,\text{RLR},\,\text{LRL}$\,.}\end{cases}

(3.42)

The numerical coefficient in Eq. (3.2) was obtained [44] by equating $x=e^{-N_{*}}$ (i.e. choosing $k_{1}=k_{*}$ ), with $N_{*}\simeq 50$ .

Upon defining the non-linearity parameter in the squeezed limit as [45]

f_{\text{NL}}^{\lambda_{1}\lambda_{2}\lambda_{3}}\equiv\frac{B_{h}^{\lambda_{1}\lambda_{2}\lambda_{3}}(k_{1},\,k_{2},\,k_{3})|_{k_{3}\ll k_{1}\simeq k_{2}}}{(12/5)P_{\zeta}(k_{1})P_{\zeta}(k_{3})}\,,

(3.43)

one finds for the leading contributions (i.e. those due to the RRR, RRL, LLR and LLL helicity combinations):

f_{\text{NL}}^{[1]}\simeq-908\,\epsilon^{2}\,\frac{m_{h}^{2}}{H^{2}}\,\left(7-\frac{2Q}{\Gamma}\right)\,.

(3.44)

We can see from Eqs. (3.44) and (3.27) that the amplitude of non-Gaussianity is sensitive to both the tensor-to-scalar ratio and the mass. One can safely maximise $r$ by setting it equal to $0.05$ ⁷⁷7The current bound from CMB data, obtained under the assumption of single field slow roll inflation is $0.056$ [46]; this bound weakens slightly in models, such as the one under study, where the relation $n_{T}=-r/8$ no longer holds.. The quantity $m_{h}^{2}/H_{*}^{2}$ can be chosen as large as $0.003$ ⁸⁸8This stems from requiring $\frac{2m_{h}^{2}}{3H_{*}^{2}}\ln(x)\ll 1$ , a condition necessary to prevent perturbation theory from breaking down (see e.g. [44]).. With these values for $r$ and $m_{h}^{2}/H^{2}_{*}$ in place, one obtains the following set of predictions:

f_{\text{NL}}^{[1]}\simeq-3\times 10^{-5}\,\left(7-\frac{2Q}{\Gamma}\right)\,,\quad\epsilon\simeq 0.0034\,,\quad n_{T}\simeq-0.0049\,,\quad H_{*}\simeq 5.8\times 10^{13}\text{GeV}\,,

(3.45)

where the value of $H_{*}$ was obtained from the amplitude of the scalar power spectrum, $\mathcal{P}_{\zeta}\simeq 2.1\times 10^{-9}$ [46]. An experiment like LiteBIRD would be able to measure the tensor bispectrum with an error $\Delta F_{\text{NL}}\simeq 1$ [45], hence we would need a coefficient $\left|7-\frac{2Q}{\Gamma}\right|\gtrsim 10^{5}$ , for the bispectrum in our model to be observable.

3.3 GW Bispectrum shape-function

We plot in Fig. 4 the bispectrum as a function of momenta. As per the usual convention, the shape-function is obtained by setting one of the momenta to $1$ (given our symmetric interaction term, this can be done without loss of generality) i.e. by setting $x_{i}=k_{i}/k_{1}$ , and multiplying the bispectrum by the factor $x_{1}^{2}x_{2}^{2}x_{3}^{2}$ .

One can see that, for the parameter space of interest, the shape-function is peaked in correspondence of the $k_{3}\ll k_{2}$ corner of the domain (in plotting the shape function, we have made the choice $k_{3}\leq k_{2}$ , without loss of generality). This is of course very reminiscent of the so-called local template. Such finding is expected on account of two notions: (i) the contribution we are considering is that of a very light field; (ii) the specific interactions at hand are strictly non-derivative. Moving towards a massive wave-function and/or considering derivative interactions will tilt the shape towards the equilateral (or flattened, orthogonal) configuration. This interpolation between different shapes is found, for example, also in the scalar sector of quasi single-field inflation [43].

The fact that the shape-function of (increasingly) derivative interactions favours equilateral, orthogonal or flat shapes (all have in common a peak that is far from the $k_{3}\ll k_{2}\sim 1$ corner is intuitively clear. Derivatives in $k$ space generate positive powers of momenta, which in turn suppress the shape-function in the squeezed limit (where the local template peaks). For the sake of completeness⁹⁹9Naturally, an external field with mass of order Hubble would decay within a few e-folds. We are interested here only in confirming the expectation that a massive wave-function will “move” the shape towards the equilateral configuration as compared to its light or nearly massless counterpart. The explanation of such an effect is due to the $(k\tau)^{\Theta(1)}$ extra behaviour that a massive field wave-function has w.r.t. a light particle at late times. Most of the contribution from cosmological correlators is generated at horizon crossing so that one can write $(k\tau)^{\Theta(1)}\sim(k/k_{\rm fixed})^{\Theta(1)}$ . These extra factors propagate all the way to the shape function and suppress the signal at small k values, i.e. in the squeezed configuration. , we consider in Fig. 5 the case of a significantly more massive field $h$ .

In reference to literature on the subject of the shape of gravitational non-Gaussianities, we should say that our model evades the conclusions in [49] in that the results therein pertain to multiple scalar fields whilst we are dealing with multiple spin-2 particles. The presence of the latter is crucial for the local shape to emerge. We are dealing with non-derivative interactions of interacting spin-2 fields whose massive eigenstates comprise a very light/massless field and a very massive one (with respect to the Hubble rate $H$ ). The shape we find is one of the possibilities beyond the so-called partially massless case. Such additional possibility for the “shapes of gravity” was also listed by [50] in a somewhat similar context to ours.

4 Allowed parameter space

In this section, we revisit and discuss all the conditions on the parameter space stemming from the regime of validity of our approximations, theoretical consistency and observational bounds.

4.1 Theoretical consistency

Validity of approximations and Higuchi bound

We start by defining two key dimensionless parameters:

x\equiv\sqrt{\kappa}\,\xi\,,\qquad\mu^{2}\equiv\frac{m_{h}^{2}}{H^{2}}=\frac{m^{2}\kappa\,\Gamma}{(\kappa+1)\,H^{2}}\,.

(4.1)

With these definitions, the conditions (3.7)–(3.9) can be written in the following compact form

\frac{2\,x^{2}}{1+x^{2}}<\mu^{2}\ll x\ll 1\,.\\

(4.2)

The values of $x$ and $\mu$ allowed by this inequality are shown in Fig.6.¹⁰¹⁰10Note that, were we to impose real values for $\tilde{\nu}_{T}\equiv\sqrt{9/4-m_{\gamma}^{2}/H^{2}}$ , the available region of the $(\mu,x)$ plane would be dramatically reduced, with $\mu/x\in[\sqrt{2},3/2]$ . In this study however, we allow also for an imaginary $\tilde{\nu}_{T}$ for simplicity. In this last case, the correlators are unaffected by the more massive tensor degree of freedom as it decays away, hence the first diagram in Fig. 3 remains the dominant contribution to the bispectrum.

Conservative cutoff of the EFT description

Ours is an interacting spin-2 model that is often dubbed a $\Lambda_{3}$ theory because its naïve cutoff is of the type $(m^{2}M_{f})^{1/3}$ or $(m^{2}M_{g})^{1/3}$ , whichever quantity is the smallest. Our constraints here are conservative because, depending on the background, the cutoff is typically revisited towards higher values [17].

Assuming $\kappa<1$ , we require that the energy scale at which we are operating during inflation, i.e. $H$ , does not exceed the cutoff, i.e. ¹¹¹¹11The assumption of $\kappa<1$ provides the most conservative conditions. If $\kappa>1$ , the lowest scale becomes $(m^{2}M_{g})^{1/3}$ and the condition to remain within the consistent EFT regime becomes weaker.

H<(m^{2}M_{f})^{1/3}\,.

(4.3)

Since we will later fix $H$ in order to (almost) saturate the bound on the tensor-to-scalar ratio $r$ [46], this condition can be interpreted as a bound on $m$

m^{2}>\frac{H^{3}}{M_{g}\,\sqrt{\kappa}}\,.

(4.4)

Existence of the low energy regime

The cosmological solutions for two interacting spin–2 fields satisfy $\tilde{c}=1$ and $\dot{\xi}=0$ in the regime where

\rho\ll m^{2}M^{2}\,.

(4.5)

In order to reach this regime already during inflation, we ought to impose another lower bound on $m$ , namely

m^{2}\gg\frac{3\,H^{2}\,(1+\kappa)\,(1+x^{2})}{\kappa}\;.

(4.6)

This condition is stronger than the condition on the validity of EFT (4.4) so long as $H<M_{p}$ and $x\ll 1$ , which, as we shall see, are both desirable and well justified in our scenario. By using the definition of $\mu$ (4.1), we can reformulate Eq.(4.6) to give an upper bound on a combination of the interaction parameters $\beta_{n}$ as

\Gamma\ll\frac{\mu^{2}}{3}\,,

(4.7)

where we approximated $1+x^{2}\simeq 1$ in accordance with (5.1), a regime we are going to discuss soon.

Perturbative Stability

In order to ensure the absence of scalar and vector instabilities, as detailed in (2.40) and (2.30), one ought to impose

\Gamma>0\,,\qquad\beta_{1}>0\,.

(4.8)

We remind the reader that in order to prevent an effective cosmological constant arising from the two-metric interaction, we fix the mass parameters so that

\rho_{m,g}=\rho_{m,f}=0\,.

(4.9)

4.2 Consistency with observations

Gravitational coupling constants

The gravitational constant for cosmology (2.43) differs from the gravitational constant for Newton’s law, which is determined solely by the coupling constant $M_{g}$ (see e.g.[51]). Since the value of $G_{N}$ is determined observationally with $\mathcal{O}(10^{-5})$ relative uncertainty [52], we fix

M_{g}=M_{p}\,.

(4.10)

In cosmology however, the effective coupling constant could affect the freeze-out time of weak interactions and introduce a modification to the primordial Helium abundance [53]

\Delta Y_{p}=0.08\,\left(\frac{G_{C}}{G_{N}}-1\right)\,.

(4.11)

Recent observations of $Y_{p}$ [54, 55, 56] have an uncertainty $|\Delta Y_{p}|<0.01\,(99.7\%\,\,{\rm C.L.})$ ¹²¹²12We use the uncertainty suggested by the PDG in Ref.[57], which constrains the deviation of the Planck masses at

\left|\frac{M_{p}^{2}}{M_{eff}^{2}}-1\right|=\frac{x^{2}}{1+x^{2}}<\frac{1}{8}\,,

(4.12)

which is consistent with the results of the recent study [58]. The region where this condition is satisfied is

x<\frac{1}{\sqrt{7}}\,,

(4.13)

which is already implemented via the last inequality of (4.2).

Solar system tests

The scalar polarisation of the graviton can affect local tests of gravity. It is possible to evade such tests whenever the Vainshtein screening mechanism is active. Screening is very effective in the following regime [59, 24, 33]

\left|\frac{d\,\log\Gamma}{d\,\log\xi}\right|\gg 1\,,

(4.14)

a condition that simply gives

\left|\frac{Q}{\Gamma}-1\right|\gg 1\,.

(4.15)

In the regimes where the tensor bispectrum (3.44) is observable, i.e. $f_{\text{NL}}^{[1]}\sim\mathcal{O}(1)$ , the ratio $|Q/\Gamma|$ is already very large so that the local tests of gravity are impervious to the scalar graviton without the need for any extra assumption.

GW observations

LIGO observations of gravitational waves from binary mergers have set an upper bound on the mass of the tensor mode in theories with interacting spin-2 particles. If we were to assume that the full spin-2 field content active during inflation survives all the way to today, such constraint would affect our parameter space via [60]

m_{h}<10^{-22}{\rm eV}\sim 10^{-49}M_{p}\,.

(4.16)

Using the relation $H\simeq 10^{15}\,\sqrt{\epsilon}$ , obtained from the amplitude of the scalar power spectrum, and rewriting (4.16) as a constraint on $\mu$ , one would have

\mu\sqrt{\epsilon}<{10^{-46}}\,.

(4.17)

Such bound would be very demanding on our parameter space, either on $\mu$ (i.e. the mass of the graviton) or on $\epsilon$ . In the latter case we would have a very flat inflationary potential, to the detriment of any observational prospect for the GW spectrum. On the other hand, applying the LIGO bound would necessarily mean presupposing that the additional “ $f$ ” spin-2 sector has survived until the present time without interacting with any other field. It is perhaps more realistic to assume instead that this massive spin-2 particle does interact with other field content and may decay some time after inflation.

One potential realisation for such a scenario is to allow the spin-2 interaction parameter $m$ to vary with the dynamics of a scalar field. By dialling down the interaction in the post-inflationary universe, it is possible to have tensors that are massive at early times and are hardly affected by the LIGO bound at the present time. Such a mechanism was introduced in Ref. [61] in the context of massive spin-2 theories with varying mass [62]. With these considerations in mind, in the remainder of the paper we will disregard the constraints on the tensor mass from late-time gravitational wave observations.

5 Discussion: summary of the constraints and a concrete example

We summarise the constraints obeyed by the parameters of the model under the additional requirement that the theory generate observable GW non-Gaussianity. The conditions arising from the validity of approximations (4.2) and BBN (4.13), are simultaneously satisfied for

2\,x^{2}<\mu^{2}\ll x\ll 1/\sqrt{7}\,.

(5.1)

Before moving on to the conditions on the $\beta_{n}$ parameters, we first make use of specific relations to reduce the number of free parameters. Since the absence of an effective cosmological constant at late times (4.9) imposes $\rho_{m,g}=\rho_{m,f}=0$ , we can fix two of the parameters using the definitions in Eq.(2.6)

	$\displaystyle\beta_{2}\xi^{2}=$	$\displaystyle-\frac{1}{6}\,\left(3\,\beta_{0}+8\,\beta_{1}\xi-\beta_{4}\xi^{4}\right)\,,$
	$\displaystyle\beta_{3}\xi^{3}=$	$\displaystyle\frac{1}{2}\,\left(\beta_{0}+2\,\beta_{1}\xi-\beta_{4}\xi^{4}\right)\,.$		(5.2)

With this specification, $\beta_{0}$ and $\beta_{4}$ can be written as:

	$\displaystyle\beta_{0}=$	$\displaystyle Q-2\,\Gamma-\beta_{1}\xi\,,$
	$\displaystyle\beta_{4}\xi^{4}=$	$\displaystyle-3\,Q-\beta_{1}\xi\,,$		(5.3)

thus reducing the freedom in the $\beta_{n}$ parameters to the one in $\beta_{1}$ , $Q$ and $\Gamma$ . The constraints on these quantities can be summarised by combining the validity of the late-time attractor condition (4.7), the perturbative stability conditions (4.8), and the requirement for an observable tensor bispectrum (3.44):

0<\Gamma\ll\frac{\mu^{2}}{3}\,,\qquad\beta_{1}>0\,,\qquad\left|\frac{Q}{\Gamma}\right|\gtrsim\frac{1}{r^{2}\,\mu^{2}}\,.

(5.4)

Let us provide a concrete example of an interesting parameter space region in our model. We consider the values discussed in (3.45) with $r=0.05$ and the mass at $\mu^{2}=3\times 10^{-3}$ . Next, in compliance with (5.1), we fix $x=10^{-2}$ . True to Eq.(5.4), we consider the value $\Gamma=10^{-5}$ . For the other two constrained parameters, we take $|Q/\Gamma|=10^{5}$ and $\beta_{1}=1$ . In order to arrive at specific numbers, we set $\kappa=1/9$ , which implies $\xi=0.03$ . The hierarchy between $Q$ and $\Gamma$ is then obtained by judiciously choosing the $\beta_{n}$ parameters. There is one subtlety associated to the choice of any value for $\xi$ ¹³¹³13Although Eqs.(5.2) and (5.3) always produce parameters that do satisfy the late-time equation (2.34), one ought to verify that the solution used for the late-time value of $\xi$ is the first root to be reached. This is because the early time evolution starts from $\xi\sim 0$ and therefore only the smallest positive root of Eq.(2.34) will be dynamically accessible. The reason is that when Eq.(2.34) is satisfied, $\xi$ stops evolving. We find that, upon requiring $Q<0$ , a set of parameters consistent with Eq.(4.9) can be readily obtained.. With that taken care of, we may obtain $Q/\Gamma=-10^{5}$ by for example employing the following $\beta_{n}$ values: $\beta_{0}=-100302\times 10^{-4}\;;\beta_{1}\xi=3\times 10^{-2}\;;\beta_{2}\xi^{2}=99701\times 10^{-4}\;;\beta_{3}\xi^{3}=-199701\times 10^{-4}\;;\beta_{4}\xi^{4}=2997\times 10^{-2}\,$ , a choice for which the smallest positive root of Eq.(4.9) is indeed precisely $\xi=0.03$ .

It is important to remark here that the specific bispectrum contribution we calculate does break so-called consistency relations (CRs) and is therefore directly physical. In order to see that this is indeed the case, we can rely on a parametric proof. It suffices to show that the tensor power spectrum may not contain a certain “free” (or -adjustable- within a interval of values set by e.g. observational or unitarity constraints) parameter that appears instead in the squeezed GW bispectrum. Now, in our case it is the parameter $Q$ , see e.g. Eq. (3.31), (as opposed to the parameter $\Gamma$ , which also regulates the bispectrum but does not share this property) that plays this role: it appears in the bispectrum but, upon analysing the quadratic tensor Lagrangian (as well its quadratic corrections via the in-in formalism, see Eqs. 3.20 and 3.21), it becomes clear it may not appear in the power spectrum.

The parameter space region corresponding to the set of values we just analysed will deliver GW non-Gaussianity that is within reach for both the LiteBIRD satellite mission and CMB-S4 experiments [45]. Next generation CMB probes provide an excellent handle on the early universe particle content. The same is in general true for probes of gravity at much smaller scales, including laser interferometer. For the purposes of testing the model at hand though, it is large scales that will deliver the most stringent constraints.

6 Conclusions

Current and next generation cosmological probes will put our ideas on the early (and late) universe to the test. In this paper we have focused on an inflationary model with a rich field content, going beyond the minimal single-field slow-roll scenario. It is well-known that we may test for a rich inflationary content across all scales. At large scales, CMB experiments as well as astronomical surveys of the galaxy distribution promise to zoom in on inflationary (self)interactions by setting strong constraints on non-Gaussianity parameters.

The improvements with respect to existing bounds will be especially remarkable for gravitational wave non-Gaussianities, with a projected $\sigma_{f^{\rm tens}_{\rm NL}}\sim 1$ from both LiteBIRD and CMB-S4 experiments. Large scale structure constraints (e.g. from Euclid, SKA, LSST) will also reach $\sigma_{f_{\rm NL}}\sim 1$ for the scalar sector. Remarkably, testing gravitational waves at intermediate and small frequencies is also guaranteed to teach us a lot about the early universe. There exist several cosmological sources of primordial gravitational waves, from those produced during pre-heating dynamics, to those associated to first order phase transitions, to GWs due to energy loss from cosmic strings. However, perhaps no cosmological mechanism may claim GWs as a universal prediction more than inflation.

Given that the minimal single-field slow-roll paradigm is typically associated with a slightly red-tilted GW spectrum, a mere GW detection¹⁴¹⁴14 That is, the detection of a GW event whose source may be clearly identified as being of cosmological rather than astrophysical origin. in the intermediate (i.e. within reach of pulsar timing arrays) and small frequency range (i.e. in the LISA, LIGO, Einstein Telescope, Taiji, Cosmic Explorer and BBO band) is generally seen as a smoking gun for multi-field dynamics.

The model we studied in this work however, does not exhibit a GW spectrum that is detectable at small scales. To identify and test its most salient features, we relied on its large scale signatures. The field content consists of minimal ingredients (massless spin-2 + scalar) with the addition of an “extra” spin-2 field. Spinning fields are particularly interesting in view, for example, of their non-analytic momentum scaling and extra angular behaviour associated to their (squeezed) non-Gaussian signal [15]. Starting with spin-2 particles, it is also possible to source GW already at linear level.

The cosmological imprints of inflationary fields with multiple spin-2 content have been tackled in the literature from several points of view. It is possible, for example [63, 64, 65, 66], to rely on an EFT description that assumes a well-behaved FRW background and starts from a Lagrangian from fluctuations atop the classical solution. A more conservative approach, more in the spirit of our present work, consists in employing the only available fully non-linear theory of interacting spin-2 fields in the inflationary context (see e.g. [67, 23]).

Our work in this manuscript has gone beyond the results in [67, 23] in manifold ways. One of the key differences is our use of the full FRW solutions (vis-à-vis the de Sitter approximation found in previous works) for the Lagrangian, which revealed a much wider parameter space, sufficiently flexible to grant detectable levels of GW non-Gaussianity. Our thorough study of constraints stemming from theoretical consistency of the FRW solution is also a novel addition in the inflationary context.

It remains to be seen whether, upon considering non-minimal coupling to the inflaton, the same field content may produce a GW signal that is testable also at small scales. This is an interesting direction we plan to pursue in the future, full in the knowledge that direct consistent coupling to “matter” is possible only under rather stringent conditions [68, 69] in the context of fully non-linear models.

Acknowledgments

MF would like to acknowledge support from the “Atracción de Talento” grant 2019-T1/TIC15784. AEG is supported by a Dennis Sciama Fellowship at the University of Portsmouth.

Appendix A Perturbed action up to cubic order

Using FLRW ansatze for both background metrics (2.4), we introduce tensor perturbations via:

	$\displaystyle ds_{g}^{2}$	$\displaystyle=-N^{2}dt^{2}+a^{2}\left(\delta_{ij}+2\,h_{ij}\right)dx^{i}dx^{j}\,,$
	$\displaystyle ds_{f}^{2}$	$\displaystyle=\xi^{2}\left[-N^{2}\tilde{c}^{2}dt^{2}+a^{2}\left(\delta_{ij}+2\,\gamma_{ij}\right)dx^{i}dx^{j}\right]\,,$		(A.1)

where $h_{ij}$ and $\gamma_{ij}$ are transverse and traceless with respect to the 3-d metric $\delta_{ij}$ .

A.1 Kinetic terms

We calculate the Einstein-Hilbert terms for the two metrics as

	$\displaystyle\sqrt{-g}R[g]=$	$\displaystyle a^{3}\left[\mathcal{L}_{g{\rm-kin}}^{(0)}+\mathcal{L}_{g{\rm-kin}}^{(2)}+\mathcal{L}_{g{\rm-kin}}^{(3)}+\dots\right]\,,$
	$\displaystyle\sqrt{-f}R[f]=$	$\displaystyle a^{3}\xi^{4}\tilde{c}\left[\mathcal{L}_{f{\rm-kin}}^{(0)}+\mathcal{L}_{f{\rm-kin}}^{(2)}+\mathcal{L}_{f{\rm-kin}}^{(3)}+\dots\right]\,,$		(A.2)

where

$\displaystyle\mathcal{L}_{g{\rm-kin}}^{(0)}\equiv$	$\displaystyle-\frac{6\,\dot{a}^{2}}{a^{2}N}\,,$
$\displaystyle\mathcal{L}_{g{\rm-kin}}^{(2)}\equiv$	$\displaystyle\frac{1}{4}\,\dot{h}{}_{ij}\dot{h}^{ij}-\left(\dot{H}+\frac{3}{2}\,H^{2}\right)h{}_{ij}h^{ij}-\frac{1}{4\,a^{2}}\,\partial_{k}h^{ij}\partial^{k}h{}_{ij}\,,$
$\displaystyle\mathcal{L}_{g{\rm-kin}}^{(3)}\equiv$	$\displaystyle-\frac{1}{2}\,h^{ij}\dot{h}{}_{i}{}^{k}\dot{h}{}_{jk}+\left(\frac{2}{3}\,\dot{H}+H^{2}\right)\,h{}_{i}{}^{k}h^{ij}h{}_{jk}+\frac{1}{4\,a^{2}}\,h^{ij}\left[\partial_{i}h^{kl}\partial_{j}h{}_{kl}-2\,\partial_{j}h{}_{kl}\partial^{l}h{}_{i}{}^{k}+2\,\partial_{l}h{}_{jk}\,\partial^{l}h{}_{i}{}^{k}\right]\,,$
$\displaystyle\mathcal{L}_{f{\rm-kin}}^{(0)}\equiv$	$\displaystyle-\frac{6\,(\xi\,\dot{a}+a\,\dot{\xi})^{2}}{a^{2}\tilde{c}^{2}N\,\xi^{4}}\,,$
$\displaystyle\mathcal{L}_{f{\rm-kin}}^{(2)}\equiv$	$\displaystyle\frac{1}{4\,\tilde{c}^{2}\xi^{2}}\dot{\gamma}{}_{ij}\dot{\gamma}^{ij}-\left(\frac{\dot{H}_{f}}{\tilde{c}\,\xi}+\frac{3}{2}\,H_{f}^{2}\right)\gamma{}_{ij}\gamma^{ij}-\frac{1}{4\,a^{2}\xi^{2}}\,\partial_{k}\gamma^{ij}\partial^{k}\gamma{}_{ij}\,,$
$\displaystyle\mathcal{L}_{f{\rm-kin}}^{(3)}\equiv$	$\displaystyle-\frac{1}{2\,\tilde{c}^{2}\xi^{2}}\gamma^{ij}\dot{\gamma}{}_{i}{}^{k}\dot{\gamma}{}_{jk}+\left(\frac{2}{3}\,\frac{\dot{H}_{f}}{\tilde{c}\,\xi}+H_{f}^{2}\right)\gamma{}_{i}{}^{k}\gamma^{ij}\gamma{}_{jk}$
	$\displaystyle+\frac{1}{4\,a^{2}\xi^{2}}\,\gamma^{ij}\left[\partial_{i}\gamma^{kl}\partial_{j}\gamma{}_{kl}-2\,\partial_{j}\gamma{}_{kl}\partial^{l}\gamma{}_{i}{}^{k}+2\,\partial_{l}\gamma{}_{jk}\,\partial^{l}\gamma{}_{i}{}^{k}\right]\,,$	(A.3)

where the above relations were obtained by adding appropriate boundary terms and we used the Cayley-Hamilton theorem for a traceless $3\times 3$ tensor,

(h^{3})_{ij}=\frac{[h^{2}]}{2}\,h_{ij}+\frac{[h^{3}]}{3}\,\delta_{ij}\,,

(A.4)

with square brackets denoting the trace operation. A similar expression applies to $\gamma_{ij}$ . We also defined the Hubble rates for the two metrics as

H\equiv\frac{\dot{a}}{a\,N}\,,\qquad H_{f}\equiv\frac{1}{\tilde{c}\,\xi}\left(\frac{\dot{a}}{a\,N}+\frac{\dot{\xi}}{\xi\,N}\right)\,.

(A.5)

We note that the lapse function $N$ is kept only for the background action, as in the main text we calculate the equations of motion in the mini-superspace approximation. For the rest of the terms, we fix the time coordinate by choosing $N=1$ .

Although we are after the $f$ – $g$ interaction terms, the expansion-induced self-interaction terms from the Einstein Hilbert terms may become relevant since $h_{ij}$ and $\gamma_{ij}$ are not mass eigenstates. Moreover, we expect (at least) part of the two-metric interaction to simplify on shell when combined with the kinetic terms.

A.2 2–metric interaction terms

We now calculate the contribution from the $g$ – $f$ interaction. We define,

\sqrt{-g}\sum_{n=0}^{4}\beta_{n}\,e_{n}(\sqrt{g^{-}1f})=a^{3}\left[\mathcal{L}_{{\rm int}}^{(0)}+\mathcal{L}_{{\rm int}}^{(2)}+\mathcal{L}_{{\rm int}}^{(3)}+\dots\right]\,,

(A.6)

where

$\displaystyle\mathcal{L}_{{\rm int}}^{(0)}=$	$\displaystyle N\,\left(\rho_{m,g}+\tilde{c}\,\xi^{4}\,\rho_{m,f}\right)\,,$
$\displaystyle\mathcal{L}_{{\rm int}}^{(2)}=$	$\displaystyle-\frac{(\tilde{c}-1)}{4}\,\Gamma\,\left([h^{2}]-[\gamma^{2}]\right)-\frac{1}{4}\,\left(\rho_{m,g}\,[h^{2}]+\tilde{c}\,\xi^{4}\rho_{m,f}[\gamma^{2}]\right)+\frac{[\Gamma+(\tilde{c}-1)Q]}{8}[(h-\gamma)^{2}]\,,$
$\displaystyle\mathcal{L}_{{\rm int}}^{(3)}=$	$\displaystyle\frac{\beta_{1}(\tilde{c}-1)\,\xi}{24}\,[(h-\gamma)^{3}]-\frac{\Gamma}{8}\,\left(\frac{1}{6}\,[(h-\gamma)^{3}]+[(h-\gamma)^{2}.(h+\gamma)]-(\tilde{c}-1)\,[(h-\gamma).(h+\gamma)^{2}]\right)$
	$\displaystyle+\frac{1}{6}\left(\rho_{m,g}[h^{3}]+\tilde{c}\,\xi^{4}\rho_{m,f}[\gamma^{3}]\right)+\frac{Q}{8}\,\left(\frac{3\,\tilde{c}-1}{6}\,[(h-\gamma)^{3}]-(\tilde{c}-1)[(h-\gamma)^{2}.(h+\gamma)]\right)\,.$	(A.7)

In the above, we used the definitions (2.6) to simplify certain combinations of the $\beta_{n}$ parameters. As in the previous section, we keep the lapse $N$ only for the background contribution and set it to unity for higher order terms.

A.3 Matter sector

For the matter sector, we use a single scalar field with the general action

S=\int d^{4}x\sqrt{-g}\,P(X,\phi)\,,

(A.8)

where $X\equiv-\frac{1}{2}\partial_{\mu}\phi\,\partial^{\mu}\phi$ , which reduces to $X=\frac{\dot{\phi}^{2}}{2\,N^{2}}$ when only tensor perturbations are included. In analogy with a perfect fluid, $P(X,\phi)$ corresponds to the pressure while the energy density and sound speed are defined as:

\rho=2\,X\,P_{,X}-P\,,\qquad c_{s}^{2}=\frac{P_{,X}}{2\,X\,P_{,XX}+P_{,X}}\,,

(A.9)

with

P_{,X}\equiv\frac{\partial\,P(X,\phi)}{\partial X}\,,\qquad P_{,XX}\equiv\frac{\partial^{2}\,P(X,\phi)}{\partial X^{2}}\,,\qquad P_{,\phi}\equiv\frac{\partial\,P(X,\phi)}{\partial\phi}\,.

(A.10)

Below are some relations that we use to derive the background equations in Sec. 2.1:

$\displaystyle\frac{\partial P}{\partial\dot{\phi}}=$	$\displaystyle\frac{P+\rho}{\dot{\phi}}\,,$
$\displaystyle\frac{\partial P}{\partial N}=$	$\displaystyle-\frac{P+\rho}{N}\,,$
$\displaystyle\frac{1}{N}\,\frac{d}{dt}\left(\frac{\dot{\phi}}{N}\right)=$	$\displaystyle\frac{c_{s}^{2}\,\dot{\phi}}{N^{2}(P+\rho)}\left(\dot{\rho}-\dot{\phi}\,\rho\right)\,,$
$\displaystyle\dot{P}=$	$\displaystyle\dot{\phi}\,P_{,\phi}+c_{s}^{2}\left(\dot{\rho}-\dot{\phi}\,\rho\right)\,.$	(A.11)

We thus have

\sqrt{-g}\,P(X,\phi)=a^{3}N\left(1-\frac{1}{4}\,[h^{2}]+\frac{1}{6}\,[h^{3}]+\dots\right)P\,.

(A.12)

A.4 Total action

We now combine all the contributions from the kinetic (A.3), interaction (A.7) and (A.9) and calculate the total action up to cubic order in tensor perturbations

S=\int d^{4}x\,a^{3}\left(\mathcal{L}_{\rm tot}^{(0)}+\mathcal{L}_{\rm tot}^{(2)}+\mathcal{L}_{\rm tot}^{(3)}+\dots\right)\,,

(A.13)

where the explicit expression for the background is presented in Eq. (2.5) and we have

$\displaystyle{\mathcal{L}}^{(2)}_{\rm tot}=$	$\displaystyle\frac{M_{g}^{2}}{8}\,\left(\dot{h}_{ij}\dot{h}^{ij}-\frac{1}{a^{2}}\partial_{k}h^{ij}\partial^{k}h_{ij}\right)+\frac{M_{f}^{2}\xi^{2}\tilde{c}}{8}\,\left(\frac{1}{\tilde{c}^{2}}\dot{\gamma}_{ij}\dot{\gamma}^{ij}-\frac{1}{a^{2}}\partial_{k}\gamma^{ij}\partial^{k}\gamma_{ij}\right)$
	$\displaystyle-\frac{m^{2}M^{2}}{8}\,\left[\Gamma+(\tilde{c}-1)\,Q\right]\,(h_{ij}-\gamma_{ij})(h^{ij}-\gamma^{ij})\,.$	(A.14)
$\displaystyle{\mathcal{L}}^{(3)}_{\rm tot}=$	$\displaystyle-\frac{M_{g}^{2}}{4}h_{ij}\left[\dot{h}^{ik}\dot{h}_{k}^{\;\;j}-\frac{1}{2\,a^{2}}\Big{(}\partial^{i}h^{kl}\partial^{j}h{}_{kl}-2\,\partial^{i}h_{kl}\partial^{l}h^{jk}+2\,\partial_{l}h^{i}_{\;\;k}\,\partial^{l}h^{jk}\Big{)}\right]$
	$\displaystyle-\frac{M_{f}^{2}\xi^{2}\,\tilde{c}}{4}\,\gamma_{ij}\left[\frac{1}{\tilde{c}^{2}}\,\dot{\gamma}^{ik}\dot{\gamma}_{k}^{\;\;j}-\frac{1}{2\,a^{2}}\Big{(}\partial^{i}\gamma^{kl}\partial^{j}\gamma{}_{kl}-2\,\partial^{i}\gamma_{kl}\partial^{l}\gamma^{jk}+2\,\partial_{l}\gamma^{i}_{\;\;k}\,\partial^{l}\gamma^{jk}\Big{)}\right]$
	$\displaystyle+\frac{m^{2}M^{2}}{48}\left(\Gamma-(\tilde{c}+1)\,Q+2\,\xi^{2}(\tilde{c}-1)\beta_{2}\right)[(h-\gamma)^{3}]$
	$\displaystyle+\frac{m^{2}M^{2}}{8}\left(\Gamma+(\tilde{c}-1)\,Q\right)\,[(h-\gamma)^{2}(h+\gamma)]\,.$	(A.15)

References

[1] D. H. Lyth and A. Riotto, Particle physics models of inflation and the cosmological density perturbation, Phys. Rept. 314 (1999) 1 [hep-ph/9807278].
[2] E. Pajer and M. Peloso, A review of Axion Inflation in the era of Planck, Class. Quant. Grav. 30 (2013) 214002 [1305.3557].
[3] CMB-S4 collaboration, CMB-S4 Science Book, First Edition, 1610.02743.
[4] A. Lee, P. A. R. Ade, Y. Akiba, D. Alonso, K. Arnold, J. Aumont et al., Litebird: an all-sky cosmic microwave background probe of inflation, Bulletin of the AAS 51 (2019) .
[5] A. A. Starobinsky, A New Type of Isotropic Cosmological Models Without Singularity, Phys. Lett. B 91 (1980) 99.
[6] Euclid Theory Working Group collaboration, Cosmology and fundamental physics with the Euclid satellite, Living Rev. Rel. 16 (2013) 6 [1206.1225].
[7] LSST Science, LSST Project collaboration, LSST Science Book, Version 2.0, 0912.0201.
[8] S. Camera, M. G. Santos and R. Maartens, Probing primordial non-Gaussianity with SKA galaxy redshift surveys: a fully relativistic analysis, Mon. Not. Roy. Astron. Soc. 448 (2015) 1035 [1409.8286].
[9] SKA collaboration, Cosmology with Phase 1 of the Square Kilometre Array: Red Book 2018: Technical specifications and performance forecasts, Publ. Astron. Soc. Austral. 37 (2020) e007 [1811.02743].
[10] J. B. Muñoz, Y. Ali-Haïmoud and M. Kamionkowski, Primordial non-gaussianity from the bispectrum of 21-cm fluctuations in the dark ages, Phys. Rev. D 92 (2015) 083508 [1506.04152].
[11] M. M. Anber and L. Sorbo, Naturally inflating on steep potentials through electromagnetic dissipation, Phys. Rev. D 81 (2010) 043534 [0908.4089].
[12] N. Bartolo and G. Orlando, Parity breaking signatures from a Chern-Simons coupling during inflation: the case of non-Gaussian gravitational waves, JCAP 07 (2017) 034 [1706.04627].
[13] D. Baumann and L. McAllister, Inflation and String Theory, Cambridge Monographs on Mathematical Physics. Cambridge University Press, 5, 2015, 10.1017/CBO9781316105733, [1404.2601].
[14] A. Kehagias and A. Riotto, High Energy Physics Signatures from Inflation and Conformal Symmetry of de Sitter, Fortsch. Phys. 63 (2015) 531 [1501.03515].
[15] N. Arkani-Hamed and J. Maldacena, Cosmological Collider Physics, 1503.08043.
[16] J. Martin, C. Ringeval and V. Vennin, Encyclopædia Inflationaris, Phys. Dark Univ. 5-6 (2014) 75 [1303.3787].
[17] C. de Rham, G. Gabadadze and A. J. Tolley, Resummation of Massive Gravity, Phys. Rev. Lett. 106 (2011) 231101 [1011.1232].
[18] S. F. Hassan and R. A. Rosen, Bimetric Gravity from Ghost-free Massive Gravity, JHEP 02 (2012) 126 [1109.3515].
[19] K. Hinterbichler and R. A. Rosen, Interacting Spin-2 Fields, JHEP 07 (2012) 047 [1203.5783].
[20] C. de Rham, Massive Gravity, Living Rev. Rel. 17 (2014) 7 [1401.4173].
[21] N. Boulanger, T. Damour, L. Gualtieri and M. Henneaux, Inconsistency of interacting, multigraviton theories, Nucl. Phys. B 597 (2001) 127 [hep-th/0007220].
[22] C. de Rham and G. Gabadadze, Generalization of the Fierz-Pauli Action, Phys. Rev. D 82 (2010) 044020 [1007.0443].
[23] E. Dimastrogiovanni, M. Fasiello and G. Tasinato, Probing the inflationary particle content: extra spin-2 field, JCAP 08 (2018) 016 [1806.00850].
[24] A. De Felice, A. E. Gümrükçüoğlu, S. Mukohyama, N. Tanahashi and T. Tanaka, Viable cosmology in bimetric theory, JCAP 06 (2014) 037 [1404.0008].
[25] A. E. Gumrukcuoglu, C. Lin and S. Mukohyama, Cosmological perturbations of self-accelerating universe in nonlinear massive gravity, JCAP 03 (2012) 006 [1111.4107].
[26] A. De Felice, A. E. Gumrukcuoglu and S. Mukohyama, Massive gravity: nonlinear instability of the homogeneous and isotropic universe, Phys. Rev. Lett. 109 (2012) 171101 [1206.2080].
[27] D. Comelli, M. Crisostomi and L. Pilo, Perturbations in Massive Gravity Cosmology, JHEP 06 (2012) 085 [1202.1986].
[28] M. Fasiello and A. J. Tolley, Cosmological perturbations in Massive Gravity and the Higuchi bound, JCAP 11 (2012) 035 [1206.3852].
[29] M. Fasiello and A. J. Tolley, Cosmological Stability Bound in Massive Gravity and Bigravity, JCAP 12 (2013) 002 [1308.1647].
[30] D. Comelli, M. Crisostomi, F. Nesti and L. Pilo, FRW Cosmology in Ghost Free Massive Gravity, JHEP 03 (2012) 067 [1111.1983].
[31] D. Comelli, M. Crisostomi and L. Pilo, FRW Cosmological Perturbations in Massive Bigravity, Phys. Rev. D 90 (2014) 084003 [1403.5679].
[32] F. Koennig, Y. Akrami, L. Amendola, M. Motta and A. R. Solomon, Stable and unstable cosmological models in bimetric massive gravity, Phys. Rev. D 90 (2014) 124014 [1407.4331].
[33] M. Kenna-Allison, A. E. Gümrükçüoǧlu and K. Koyama, Viability of bigravity cosmology, Phys. Rev. D 99 (2019) 104032 [1812.05496].
[34] K. Aoki, K.-i. Maeda and R. Namba, Stability of the Early Universe in Bigravity Theory, Phys. Rev. D 92 (2015) 044054 [1506.04543].
[35] E. Mortsell and J. Enander, Scalar instabilities in bimetric gravity: The Vainshtein mechanism and structure formation, JCAP 10 (2015) 044 [1506.04977].
[36] M. Högås, F. Torsello and E. Mörtsell, On the stability of bimetric structure formation, JCAP 04 (2020) 046 [1910.01651].
[37] H. Hui et al., BICEP Array: a multi-frequency degree-scale CMB polarimeter, Proc. SPIE Int. Soc. Opt. Eng. 10708 (2018) 1070807 [1808.00568].
[38] SPT-3G collaboration, SPT-3G: A Next-Generation Cosmic Microwave Background Polarization Experiment on the South Pole Telescope, Proc. SPIE Int. Soc. Opt. Eng. 9153 (2014) 91531P [1407.2973].
[39] Simons Observatory collaboration, The Simons Observatory: Science goals and forecasts, JCAP 02 (2019) 056 [1808.07445].
[40] CMB-S4 collaboration, CMB-S4 Science Book, First Edition, 1610.02743.
[41] J. Crowder and N. J. Cornish, Beyond LISA: Exploring future gravitational wave missions, Phys. Rev. D 72 (2005) 083005 [gr-qc/0506015].
[42] S. Kawamura et al., The Japanese space gravitational wave antenna DECIGO, Class. Quant. Grav. 23 (2006) S125.
[43] X. Chen and Y. Wang, Quasi-Single Field Inflation and Non-Gaussianities, JCAP 04 (2010) 027 [0911.3380].
[44] S. Endlich, B. Horn, A. Nicolis and J. Wang, Squeezed limit of the solid inflation three-point function, Phys. Rev. D 90 (2014) 063506 [1307.8114].
[45] M. Shiraishi, Tensor Non-Gaussianity Search: Current Status and Future Prospects, Front. Astron. Space Sci. 6 (2019) 49 [1905.12485].
[46] Planck collaboration, Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10 [1807.06211].
[47] M. Celoria and S. Matarrese, Primordial Non-Gaussianity, Proc. Int. Sch. Phys. Fermi 200 (2020) 179 [1812.08197].
[48] D. Babich, P. Creminelli and M. Zaldarriaga, The Shape of non-Gaussianities, JCAP 08 (2004) 009 [astro-ph/0405356].
[49] L. Bordin, P. Creminelli, M. Mirbabayi and J. Noreña, Tensor Squeezed Limits and the Higuchi Bound, JCAP 09 (2016) 041 [1605.08424].
[50] G. Goon, K. Hinterbichler, A. Joyce and M. Trodden, Shapes of gravity: Tensor non-Gaussianity and massive spin-2 fields, JHEP 10 (2019) 182 [1812.07571].
[51] M. Fasiello and R. H. Ribeiro, Mild bounds on bigravity from primordial gravitational waves, JCAP 07 (2015) 027 [1505.00404].
[52] Particle Data Group collaboration, Review of Particle Physics, PTEP 2020 (2020) 083C01.
[53] S. M. Carroll and E. A. Lim, Lorentz-violating vector fields slow the universe down, Phys. Rev. D 70 (2004) 123525 [hep-th/0407149].
[54] Y. I. Izotov, T. X. Thuan and N. G. Guseva, A new determination of the primordial He abundance using the He i $\lambda$ 10830 Å emission line: cosmological implications, Mon. Not. Roy. Astron. Soc. 445 (2014) 778 [1408.6953].
[55] E. Aver, K. A. Olive and E. D. Skillman, The effects of He I $\lambda$ 10830 on helium abundance determinations, JCAP 07 (2015) 011 [1503.08146].
[56] M. Valerdi, A. Peimbert, M. Peimbert and A. Sixtos, Determination of the Primordial Helium Abundance Based on NGC 346, an H ii Region of the Small Magellanic Cloud, Astrophys. J. 876 (2019) 98 [1904.01594].
[57] Particle Data Group collaboration, Review of Particle Physics, PTEP 2020 (2020) 083C01.
[58] M. Högås and E. Mörtsell, Constraints on bimetric gravity from Big Bang nucleosynthesis, 2106.09030.
[59] A. De Felice, T. Nakamura and T. Tanaka, Possible existence of viable models of bi-gravity with detectable graviton oscillations by gravitational wave detectors, PTEP 2014 (2014) 043E01 [1304.3920].
[60] LIGO Scientific, Virgo collaboration, Tests of general relativity with GW150914, Phys. Rev. Lett. 116 (2016) 221101 [1602.03841].
[61] K. Hinterbichler, J. Stokes and M. Trodden, Cosmologies of extended massive gravity, Phys. Lett. B 725 (2013) 1 [1301.4993].
[62] Q.-G. Huang, Y.-S. Piao and S.-Y. Zhou, Mass-Varying Massive Gravity, Phys. Rev. D 86 (2012) 124014 [1206.5678].
[63] L. Bordin, P. Creminelli, A. Khmelnitsky and L. Senatore, Light Particles with Spin in Inflation, JCAP 10 (2018) 013 [1806.10587].
[64] L. Iacconi, M. Fasiello, H. Assadullahi, E. Dimastrogiovanni and D. Wands, Interferometer Constraints on the Inflationary Field Content, JCAP 03 (2020) 031 [1910.12921].
[65] L. Iacconi, M. Fasiello, H. Assadullahi and D. Wands, Small-scale Tests of Inflation, JCAP 12 (2020) 005 [2008.00452].
[66] A. Malhotra, E. Dimastrogiovanni, M. Fasiello and M. Shiraishi, Cross-correlations as a Diagnostic Tool for Primordial Gravitational Waves, JCAP 03 (2021) 088 [2012.03498].
[67] M. Biagetti, E. Dimastrogiovanni and M. Fasiello, Possible signatures of the inflationary particle content: spin-2 fields, JCAP 10 (2017) 038 [1708.01587].
[68] C. de Rham, L. Heisenberg and R. H. Ribeiro, On couplings to matter in massive (bi-)gravity, Class. Quant. Grav. 32 (2015) 035022 [1408.1678].
[69] C. de Rham, L. Heisenberg and R. H. Ribeiro, Ghosts and matter couplings in massive gravity, bigravity and multigravity, Phys. Rev. D 90 (2014) 124042 [1409.3834].