Primordial Non-Gaussianity in Supersolid Inflation

We study primordial non-gaussianity in supersolid inflation. The dynamics of supersolid is formulated in terms of an effective field theory based on four scalar fields with a shift symmetric action minimally coupled with gravity. In the scalar sector, there are two phonon-like excitations with a kinetic mixing stemming from the completely spontaneous breaking of diffeomorphism. In a squeezed configuration, $f_{\text{NL}}$ of scalar perturbations is angle dependent and not proportional to slow-roll parameters showing a blunt violation of the Maldacena consistency relation. Contrary to solid inflation, the violation persists even after an angular average and generically the amount of non-gaussianity is significant. During inflation, non-gaussianity in the TSS and TTS sector is enhanced in the same region of the parameters space where the secondary production of gravitational waves is sizeable enough to enter in the sensitivity region of LISA, while the scalar $f_{\text{NL}}$ is still within the current experimental limits.


Introduction
We have minimal knowledge of the Universe before radiation domination. The most convincing solution of the horizon and flatness problems of the hot big bang model is to assume that the Universe had gone through an early phase of accelerated expansion driven by some sort of "matter" (inflaton). The only glimpse of the inflationary phase is the seed of primordial perturbations that gave rise to structure formation in the Universe via gravitational instability. Interestingly, many of the inflationary phase features are determined by the spontaneous symmetry breaking during inflation. In the simplest case of single clock inflation [1] 4-dimensional diffeomorphisms are broken down to 3-dimensional diffeomorphisms of the hypersurface φ =constant, where φ is the inflaton field; the Weinberg theorem holds [2], perturbations are adiabatic, and the detailed theoretical predictions for primordial non-gaussianity [3,4] are not affected by the reheating of the Universe. The downside is that primordial non-gaussianity turns out to be very small, and the amplitude of the stochastic background gravitational waves generated during inflation is tiny and out of reach for LISA. When more fields are present, symmetry breaking pattern can be completely different, changing the predictions significantly for primordial non-gaussianity. In the present work, we propose an effective theory description based on up to four scalar fields that allow us to study the symmetry breaking pattern during inflation systematically [5,6]. When diffeomorphisms are completely broken down to a global group required for the existence of de Sitter background space-time, the fluctuations of the scalar fields around their VEVs can be associated with the phonon-like excitations of a self-gravitating medium with the properties of a supersolid [7][8][9]; the breaking of spatial diffs by a solid-like medium was considered in [10]. Such a medium also provides a mass to the graviton [5,6,[11][12][13] via gravitational Higgs mechanism. A similar breaking pattern during inflation was already considered in [14] with results similar to solid inflation. Our analysis shows that in fact the supersolid case is intrinsically rather different from solid inflation, as will be discussed in the rest of the paper. In the scalar sector, two independent modes are present and mix from the beginning, in the Bunch-Davis vacuum, till the end of inflation. Consequently, perturbations are not anymore purely adiabatic.
The presence of non-adiabatic perturbations and anisotropic stress of the solid component of the medium leads to violation of the Weinberg theorem. In this sense, the option of inflating and forgetting is not available and reheating must be taken into account to determine the seed of primordial perturbations to be used at large scales as initial conditions for the standard radiation dominated phase of the Universe's evolution. The curvature of constant energy density hypersurface ζ does not coincide at superhorizon scales with R which represents the curvature of hypersurface orthogonal to medium velocity, and both of them are not perfectly conserved at superhorizon scales; further, the Maldacena consistency relation [3] is not satisfied. One can show that [9] in the instantaneous reheating approximation, the seed of primordial perturbations are almost adiabatic in agreement with the most recent CMB data [15]. On a more phenomenological side, we point out the existence of a region in the parameter space of supersolid inflation such that the primordial tensor power spectrum (PS) and the related stochastic background of gravitational waves could be significantly enhanced with a blue-tilted spectral index via secondary production due to the non-linear coupling among scalars and gravitons [9]. Simultaneously, in the same region, non-gaussianity related to tensor fields can be noticeably enhanced, keeping the scalar prediction within the current experimental constraints. Finally, in solid and probably also in supersolid inflation, as discussed in [16] one should check how fast an anisotropic background gets to a Freedman-Robertson-Walker solution. In particular a statistical anisotropy can be generated [17] due to infrared modes that modify the background; such an effect is sizeable when inflation lasts for more than the minimal number of e-folds needed to generate the CMB anisotropies. In this work, we suppose a suitable duration of inflation that allows to transmit the correct amount of isotropy, and constraints on parameters will be obtained considering the current experimental limits on non-gaussianity only. We leave the detailed analysis of this interesting feature for future work.
The outline of the paper is the following. After a brief introduction to the effective description of a supersolid given in section 2, the parameters entering in the quadratic and cubic action are discussed in 3, while section 4 is devoted to the discussion of the power spectra, resuming the main results found in [9]. Sections 5, 6, 7 and 8 contain the analysis of primordial non-gaussianity. Finally, in section 9 we briefly analyze the phenomenological implications given by supersolid inflation. Our conclusions are drawn in section 10.

Supersolids
Single field inflation is the simplest choice in the vast menu of inflationary models. It successfully addresses all the drawbacks of the hot big bang model and predicts a tiny level of primordial non-gaussianity (PNG) and very low tensor to scalar ratio in agreement with the lower bound from CMB [15,18]. Why then should one consider more complicated models?
The answer is related to the very different PNG predictions. In single clock inflation, the inflaton background value ϕ(t) spontaneously breaks the 4d-diffeomorphisms of general relativity (GR) down to 3d-diffeomorphisms of the hypersurface ϕ =const. The symmetry breaking pattern plays a crucial role in many aspects of inflation and primordial nongaussianity. For instance, in the squeezed limit, the residual symmetry group determines the form of f NL [3,[19][20][21][22][23][24]; see also [25] for a recent discussion. Things are different when 4d-diffeomorphisms are completely broken, leaving only a global symmetry as a leftover. The minimal number of scalar fields, that can implement such a scenario, is four: ϕ A , A = 0, 1, 2, 3 with a background value, during inflation, of the form ϕ 0 =φ(t), ϕ l = x l , l = 1, 2, 3 ,ḡ µν = a(t) 2 η µν . (2.1) To allow a FLRW (Freedman Lemaitre Robertson Walker) background solution, the action describing the scalar fields dynamics needs to be symmetric under internal rotations and shift transformations ϕ l → R l m ϕ m , R R R t R R R = 1 , ϕ A → ϕ A + c A A = 0, 1, 2, 3 . (2. 2) The building block matrix is used to write down the action where (2.5) and u µ plays the role of a normalized time-like four-velocity such that u µ ∂ µ ϕ l = 0 The action (2.4) can be interpreted as the relativistic generalization of the low-energy effective Lagrangian describing homogeneous and isotropic supersolids at zero-temperature [7,26]. Such an action is the most general at leading order in a derivative expansion compatible with the given internal symmetries, see [9] for more details. For the benefit of readers we point out that, besides the presence of ϕ 0 , our choice for the independent operators is slightly different from [10] 1 . From the energy momentum tensor (EMT) of the scalars' action, one can infer the energy density and pressure The action (2.4) reduces to the one of a solid when the following operators are sent to zero: y, χ, w X , w Y , w Z ; that is equivalent to remove the scalar ϕ 0 .

Mass Parameters and Cubic Vertices
At the background level, the EMT of the scalar fields is the one of a perfect fluid with energy density and pressure given bȳ with U and its derivatives evaluated at the background value of the independent operators The simplest way to identify the parameters entering in the quadratic and cubic expansion of (2.4) is to work in unitary gauge where the scalar fields are frozen to their background values and all perturbations are in the metric g µν = a 2 (η µν + h µν ). Schematically, the expansion of the Lagrangian of the scalars has the following structure From SO(3) symmetry, we can define the following independent operators.
• Linear level: h 00 , h ii and the associated parameters, the energy densityρ and the pressurep.
In particular, the quadratic Lagrangian for scalars is given by evaluated at the background values of the operators (3.2). The EMT conservation at the background level is equivalent tō Thanks to the above relations, it is easy to realize that for a generic U , at background level, the following quantity is conserved that we call background entropy per particle, see [8,9,13]. We will consider inflation in a slow-roll (SR) regime for which It is useful to introduce the following parametrization forσ and the masses M î in terms of the Hubble parameter H = a −1 H (almost constant during slow-roll), and slow varying in time. The M i proportionality to ∼ (ρ + p) is required by dynamical stability (4.3, 4.4), while the slow varying condition ensures the plane-wave solutions on subhorizon scales, as discussed in [9]. We will often use c 2 L defined by instead of c 2 2 and we will substitute c 2 3 in terms of the adiabatic sound speed c 2 s that during SR is approximately equal to -1 (see also (4.4)) (3.13) The cubic order Lagrangian L 3 in the unitary gauge, at the leading order in the SR expansion, results where the {λ i } parameters at the leading order in SR are defined by For the spin two field, in Fourier space, we use the decomposition h ij (k) = p p ij (k) h p k in terms of the circular polarization tensors p ij (k) (p = +, ×) that are traceless ( p ii (k) = 0 ) and transverse (k i p ij (k) = 0 ). The scalars 2 Ψ and F have algebraic equations of motion and can be integrated out in favor of the proper Goldstone modes π 0 and π L . In particular, working in a quasi de Sitter background in SR regime, they result to be proportional to so that, at leading order, they can be neglected both in the quadratic and cubic action, and no operator containing such fields will be shown. Expanding (2.4) at the quadratic order in the fluctuations we get where ∆ = ∂ 2 i . In the scalar sector, there are two independent propagating scalar modes that, following [9], can be quantized by diagonalizing the quadratic action and implementing the corresponding Bunch-Davies (BD) vacuum. It is important to point out that the parametrization (3.10) is more than a convenient choice but it is forced by stability and slow-roll. Stability derived from (4.2) requires that [13] on the other hand, in quasi de Sitter,p +ρ ∼ H 2 , thus M 1 ∼ H 2 . Now, the adiabatic speed of sound (3.13) is explicitly given bȳ and during slow-roll is very close to −1. Thus, unless M 4 , M 2 and M 3 are precisely tuned, the choice (3.10) is the only possible. Finally, notice that if π 0 should play the role of an inflaton field 3 , likewise π L , then M 0 ∼ H 2 . The quadratic action (4.2) in the scalar sector describes two massless modes with a kinetic mixing. To remove the mixing, one should set c 2 b = −1 and c 2 1 = 2 c 2 0 ; unfortunately stability imposes that M 1 and then c 2 1 to be negative, while M 0 and then c 2 0 to be positive. Thus, the mixing between longitudinal phonons cannot be undone by a tuning of the parameters. The diagonalisation of (4.2) is non-trivial and it is discussed in [9]. The result is two independent longitudinal phonons with speed of sound c s1 and c s2 .
An essential consequence of (3.10) is that by consistency, imposing that the {M i } are time independent at the leading order in slow-roll andσ is strictly constant, we get some useful relations among the parameters. Namely (4.5) Finally we can express c 2 0 and c 2 1 in terms of the diagonal sound speed parameters c s1 , c s2 and c L according to [9] Thus, the number of free parameters is given by the two sound speeds (c s1 , c s2 ), the "graviton mass" 4 c 2 2 , the mass ratio c 2 b and the entropic parameterσ plus the ten parameters {λ i } of the cubic interactions constrained by the five relations (4.5). Thus we end up with the following 10 independent parameters: {σ, c 2 s1 , c 2 s2 , c 2 2 , c 2 b , λ 2, 6, 7, 9, 10 }. An explicit 3 If M0 ∼ H 2 then π0 is a spectator field and physics is completely different. 4 the parameter c 2 L is equivalently to c 2 2 by eq.(4.20).
example of a Lagrangian describing a supersolid where the constraints (4.5) are automatically incorporated is the following we have defined the combination of operators The term UΛ in U can be called a background Λ-medium due to the fact that, at the background level, has the equation of state of a cosmological constant The perfect fluid V (b) part has to be tuned to the specific functional form to ensure the presence of a quasi de Sitter phase with an almost constant factor Thus, such a rather general Lagrangian automatically satisfies eqs.(4.5) and generates a slow-roll regime with an almost constant . More details are given in Appendix F. Note that in [27] and [28] we defined a bit different potentials named simply Λ-Media, characterised by the fact that ρ Λ + p Λ = 0 worth exactly, not only on their background value 5 . An example of a supersolid that it is also a Λ-Medium, in our notations, is Let us come back to the description of the perturbations PS. One of the key features in (4.2) is that π L − π 0 mix both at super and sub horizon scales 6 . As a result, any scalar field ξ, resulting from a combination of π fields, can be written in Fourier space as a linear combination of two independent annihilation (creation) operators a are two classical solutions of the linearized equations of motion for ξ with initial conditions at early conformal 5 Another realization of an exact Λ-Medium is given in [29]. 6 This is one of the important differences with the analysis in [14] is that the kinetic mixing was treated as small. In addition, the authors setφ = a which leads to a conserved background EMT only if c 2 b = 0. Such a value of c 2 b is rather peculiar, as we will see in what follows. The correct implementation of the Stückelberg trick at the background level requires a non-trivial background for ϕ 0 satisfying (3.7). In [30] the analysis is similar but a non-minimal coupling between the scalars and curvature is also considered. time t → −∞ set by the choice of the BD vacuum. Thus, the two-point function and the related linear power spectrum for ξ reads ξ . (4.11) As discussed in [9,28], reheating 7 is likely to process primordial perturbations produced during inflation such that the initial conditions for the hot phase are determined by the curvature of constant number density of lattice sites n , given by the operator b. The corresponding curvature perturbation related to the b operator is provided in the spatiallyflat gauge at the linear order by Another important scalar field is given by the curvature perturbation of the constant ϕ 0 hypersurface, namely The two above perturbations are closely related to the Goldstone-like excitations π 0 and π L , that are the fluctuations around the VEVs responsible for the complete breaking of four-dimensional diffeomorphisms during inflation. To figure out better their roles, notice that R π 0 , when the single-field breaking pattern is considered (absence of ϕ a fields), tends to coincide with R, the curvature of constant fluid velocity surface. Contrary, when only the solid breaking pattern is present (absence of ϕ 0 ), the ζ n field tends to coincide with the ζ field, the curvature of constant energy density surfaces. Finally, let us defineR This last has no fundamental meaning, but it is useful to minimize the number of cubic interactions. All scalar perturbations can be written in terms of ζ n , R π 0 and their time derivatives, whose linear PS are given by [9] P (1) together with the cross-correlation where expressions (4.15, 4.16, 4.17) are manifestly symmetric under c s1 ←→ c s2 exchange. We have also introduced the power spectrum P of canonical single field inflation given by (4.18) These relations are strictly valid for c 2 b = 0, −1 only, where analytic solutions for the ζ n /R π modes can be found (see Appendix D, eq. (D.5)). Thus, from here below, we will consider the c b parameter fixed; while the diagonal sound speeds are conventionally ordered such that c 2 s2 < c 2 s1 ; moreover stability requires that c 2 s2 < c 2 L < c 2 s1 . In particular, in the rest of this work, we will focus on the c 2 b = −1 case. Our choice can be ascribed to the higher R π 0 sensitivity to low values of c s2 , as it can be easily inferred from fig. 1. As deeply argued in the next sections, this will enhance the tensor NG, 8 in the limit of small c s2 . We would like to clarify that formal results exposed in Appendixes D, E and table 4 will be still valid for c 2 b = 0, and differences rise only once the c s2 expansion is implemented. In order to simplify the PNG notation, we define the following dimensionless ratios: specifying how far from P ζn the single-mode component P (j) ζn or the cross-correlations ζ n -R π 0 could be. Considering that, during an almost instantaneous reheating, ζ n is continuous during the transition and can be regarded as the seed of primordial adiabatic perturbations [9], its PS is constrained by the CMB experiments and assumes the form: (4.20) 8 As well as the tensor power spectrum as argued in [9] whereĉ L is an effective longitudinal sound speed 9 which can be read off from (4.15). One can show [9] that when 10 <ĉ −5 L < 100, the supersolid is dynamical stable, c s1 and c s2 are subluminal and the ζ n power spectrum is well within the PLANCK constraints [15] by taking 10 −11 ≤ P ≤ 10 −10 . For illustrative purposes we setĉ L = 1/2 in our figures, which maximizes the allowed (c s1 , c s2 ) region. It is important to note that when c s2 is much smaller than c s1 andĉ L , and P ζn is fixed, the power spectrum (4.16) of R π 0 and the cross-correlation (4.17) become dominant due to the dependence on c s2 (see fig.1). This interesting subset of the parameter space, besides well reproducing PLANCK's data, has very distinctive features for secondary gravitational waves production [9] and for what concerns PNG as discussed in sections 5 and 7. We stress that the asymmetry in between the two c s1 and c s2 parameters starts once we imposed that only the PS of the ζ n field is the adiabatic/observed one, as predicted by the matching conditions during an instantaneous reheating phase. It is this last condition that makes the R π 0 power spectrum sensitive to small c s2 values.

Cubic Lagrangian (spatially flat gauge)
Primordial non-gaussianity is due to the interactions among the various dynamical fields during inflation and, at the leading order in perturbation theory, originates from cubic terms that can be divided into interactions among scalars (SSS), gravitons (TTT), and mixed ones: (TSS), (TTS). Thus the cubic Lagrangian can be split according to As we have discussed in section 3, there are up to ten parameters in the cubic Lagrangian; we focus on the region of parameters where c 2 s2 1 that is probably the most interesting from a phenomenological point of view. The enhancement of the size of PNG in the operators containing π 0 is due to the extra negative powers of c s2 in R π 0 and in the cross-correlation. It is convenient to associate at each cubic vertex a color code according to the power 10 of c s2 −1 of the corresponding contribution to the 3-point function of ζ n and then to f NL as shown in table 1. The color code goes from violet (no powers of c s2 −1 are present) to the red (highest negative power c s2 −8 ). Clearly, vertices labeled by a large wave-length color potentially give a large contribution to the 3-point function.
By expanding the action at the cubic order, we find in the scalar sector at leading order in 9 We used a notation typical for the 2-point function in solid (and also in some single-field) inflationary models [10] where it is present an extra factor c −5 L in adiabatic power spectrum, in our case c −5 L →ĉ −5 L . 10 Such classification is strictly valid in the case c 2 b = −1. Table 1: Operators in the cubic Lagrangian (5.1) classified according the powers c −n s2 , n = 0, ..., 8 of the corresponding contribution to ζ 3 n , in the limit c s2 1 and c 2 b = −1. When c 2 b = 0 the c s2 divergences are not as stiff as in the According to our classification scheme, for instance the vertex O 2 (orange) gives a larger contribution to the f NL than G 1 (green) but smaller than R 1 (red). The above operators are not all independent but can be related by total derivatives. For the violet operators we have the following spatial total derivative so that the V i form factors will appear in the following combinations For blue operators we get the following total derivatives The normalization of the vertices is chosen to simplify the Lagrangian (A.3) written in terms of ζn and that implies the presence in the observables of the following independent combinations B 1 + B 2 , B 4 − B 1 /2 and B 5 + B 1 /2. Cubic interactions among gravitons arise both from the Lagrangian U and the Einstein-Hilbert (EH) action. We do not consider the TTT, TTS, and TSS interactions coming from EH action (see for instance [3]), being subdominant with respect to the ones coming from U during inflation. Finally, the cubic TSS Lagrangian with a single graviton has the form The TTS cubic Lagrangian with two gravitons reads The cubic interacting Lagrangian for gravitons is given by In Appendix A, all vertices X = V, B, G, O, R are given as functions of the derivatives of the U Lagrangian. For the benefit of the reader interested in the details of the computation of PNG, the cubic Lagrangian is also written in Fourier space in terms of ζ n and R π 0 . By using (3.7) one has thatφ = a 1−3 c 2 b and the vertices generically denoted by X can be split into a time-independent part X 0 and a part proportional toσ whose time dependence is dictated by c 2 b according to where X σ is an X-dependent constant. The case c 2 b = −1 is special, being the full vertex time-independent. Actually, all the terms proportional toσ in the cubic action reconstruct total spatial and time derivative terms, when the last relation in (4.5) is taken into account. Thus, effectively we can forget aboutσ as soon as we are interested in the bispectrum of our adiabatic-entropic modes, see Appendix B.

Effective Theory Description
To compute PNG, we need to expand the action at least at the cubic order in the perturbations. Because of the action's shift symmetry, each π A field contains a least one derivative. Focusing on the scalar sector, schematically an interaction term with n scalar fields 12 will be of the form where the n-th derivative of U with respect of its arguments U n is assumed to be proportional to U ∼ H 2 times a dimensionless parameter λ n ; typically λ n will contain a combination of slow-roll parameters. The kinetic matrix of the π fields has the general structure given in 0 for the π 0 field and κ i=L ∝ 1 + c 2 1 for π L . By introducing the canonical Consequently, we get for Λ n Depending on the value of λ n there are at least two options: democratic and special. In the democratic case λ n ∼ λ 1 and all the derivatives of U are of order while in the special case λ n ∼ λ ∼ O(1). In both cases, the lowest cutoff Λ is obtained for n = 3 As expected Λ is proportional to the cutoff scale Λ 2 = M pl m of Lorentz breaking massive gravity in flat spacetime [5,6,12,13,31,32] with the graviton mass m ∼ H. Note that in the democratic approach, the cutoff scale is higher. In solid inflation [10] it was used the special approach and only a unique combination of derivatives of U is required to be small by consistency with slow-roll, namely In this paper we assume the democratic approach. By using (4.18) and (4.20), in order to avoid strong coupling during the inflationary period, we require that Λ H which gives Take the more stringent case . Thus, the most dangerous operator is a cubic interaction with three π 0 ; the absence of strong coupling implies that leaves enough space to explore the region of interest in the parameter space with a sufficiently small c s2 sound speed.

Scalar Bispectrum
Once the cubic Lagrangian is provided, primordial non-gaussianity can be computed evaluating the various 3-point functions by using the in-in formalism. The computation is rather straightforward, so we include some of the details in Appendix C.
We will focus on the c 2 b = −1 case, probably the most interesting case from a phenomenological point of view; the reason for that is two-folded. The relevant scalar fields ζ n and R π 0 have an almost flat power spectrum when −1 ≤ c 2 b ≤ 0 and for this particular case one can work out an analytic expression for the modes defined in (4.10) that are essential for the computation of the time integrals entering in the 3-point functions.The Fourier transform of a generic 3-point function depends on the 3-momenta k 1 , k 2 and k 3 satisfying In the computation of B we have retained the leading order term in the slow-roll expansion as in [10]. However, there are cases [33] where slow-roll corrections (π L , π 0 , lapse and shift next to leading terms) of the order can give sizeable corrections, where N k i is the number of e-folds at the time t i = (c si k i ) −1 of the sound horizon crossing of the i-th momentum, while t e is the end time of inflation. In this work, for each momentum, we consider sufficiently small to neglect those corrections.

Squeezed Configurations
To avoid excessively long expressions, the main focus will be on the squeezed configurations, the most constrained one by observations, see for instance [34]. In a squeezed configuration we have | k 3 | = k L | k 1 | = k S , and θ is angle between k 1 and k 3 then | k 2 | = k S + k L cos θ + · · · . The bispectrum of ζ n in the squeezed limit has the following form . 13 The traditional choice of the factor 6/5 comes from the following parametrization (local ansatz): where ζg is the gaussian part of the ζ-field. As well known other shapes and related fNL can be defined, as the equilateral, orthogonal and folded ones.
In our case the monopole term f M is not proportional to n s − 1 ∼ O( ) as in single clock inflation, and actually f M is in general of order 0 , if we consider the democratic approach. Also, the quadrupole term f Q is generally present and of order 0 . As expected, the Maldacena consistency relation is explicitly violated being the symmetry breaking pattern different from the one of single-field inflation. The structure of different 3-point functions in a squeezed configuration is similar as discussed in Appendix D.
The structure of a 3-point function in the squeezed limit of fields almost conserved on superhorizon scales and with a scale-free spectrum can be studied in general by extending the procedure exposed in [35]. Such fields at large scales can be considered as classical stochastic variables. The long "squeezed" modes can be incorporated in the background modifying linear dynamics of the short modes. In other words, introducing a long-short (k L − k S ) modes splitting of the relevant fields, the effect of the squeezed components is to generate sort of background fields whose effect is to give rise to a non-linear and local correction of the 2-point correlation function. In our case, in the scalar sector, the cubic Lagrangian is written in terms of ζ n and R π 0 fields. The squeezing in the j-th vertex will be equivalent to ks −(ks+k L ) to the effective quadratic Lagrangian (a convolution over k L is understood). In the case of single-field inflation, there is only one scalar mode given by the comoving curvature R and the effective quadratic action coincides with the quadratic Lagrangian up to a rescaling of the spatial coordinates [19] which is a residual symmetry of single-clock inflation in the limit k L → 0 where R L (x) is the long part of the comoving curvature. The equivalence between L (2) eff. and the rescaled quadratic action is the core of the Maldacena consistency relations [3], any deviation leads to a violation of such relations. By looking at the following ratio we can predict the nature of the consistency violations. Namely, an eventual time-dependence of Q j can be related to the scale dependence of f NL in the squeezed limit producing possibly dangerous interactions. This is never the case for c 2 b = −1. Still, the residual dependence on momenta directions, also related to the vertices' derivative structure, can give primordial anisotropies that can be very interesting from a phenomenological point of view. In general, f NL will be not suppressed by slow roll parameters and, as discussed above, with several features rather different from single-clock inflation. Let us now focus on the region of the parameter space whereĉ L , c s1 ≤ 1, while c s2 1. As shown in [9], in such a region, the secondary production of gravitational waves can significantly boost the tensor power spectrum enough to enter the sensitivity region of LISA. From (4.19) and (4.6), we have the following expansion we have set c 2 b = −1, andĉ 5 L has been defined in (4.20). By using the above expressions, one can compute f NL for the bispectrum of ζ n ; the result is given in table 2.
We have introduced For what concerns the wave functions of ζ n , R π 0 and their time derivative, we can deduce the following leading order behaviour in the small parameter c s2 Note that for theR π 0 field (as defined in (4.14)) there are two competing contributions.
Given that the cross-correlation r L/0 ∼ c −3 s2 (7.3), the dominant contribution to f NL is given by the vertices with the greatest number of R π 0 . A number of features follows.
1. The contribution to f NL of violet vertices, typical of solid inflation, is little sensitive to c s2 ; indeed, in those R π 0 is absent.

2.
Blue vertices have a single R π 0 and contribute to f NL via the cross correlation term r 3. Green vertices and orange vertices, with at least two R π 0 , contribute with a cross correlation of the form r 4. A vertexR π 0 3 , even if it has threeR π 0 contributes as r 5. Similarly, the red vertex can be rewritten asR π 0 R 2 π 0 and contributes as r Clearly, in the deep c s2 1 region, the contributions from red, orange and green vertices tend to be very large compared with the experimental bounds. We come back to that in section 9. In Appendix D the reader can find total squeezed expressions valid for any value of c s2 , where the contribution of each vertex O is factorized in an angular part and in an angle-independent part M (O) . Plots of the size of the above contributions are given in Appendix H.

Equilateral Configurations
Proceeding as for the squeezed configuration, the expressions for f NL in the equilateral configuration read Note that the green interaction ζ n R 2 π 0 , in the equilateral configuration assumes an orange behavior being the f NL proportional to c −6 s2 . The same holds for the folded shape whose results are very similar and we do not report for brevity. In section 9 we check that experimental constraints on both equilateral and folded shapes can be satisfied.

Tensor Bispectrum
Let us now study 3-point function containing at least a tensor. We do not consider vertices coming from the Einstein-Hilbert of the action as discussed in section 5. Before proceeding let us recall the expression for the two-point tensor correlator with P h the gaussian scale invariant tensor PS whose expression is As shown in [9], when c s2 c s1 , the secondary GWs production tends to overwhelm the primary one in the small c s2 limit. The same mechanism gives rise to the TTT and TTS one-loop correction (work in progress). This last one will be naturally dominant in the same limit where the secondary PS dominates. In this work, we limit our analysis to the three-level contribution, whose expressions in the TTS and TSS cases can leave sizeable and very interesting signatures. In the following sections we give explicit expressions for three tensor correlators two tensors and one scalars concluding with one tensor and two scalars

T T T
Let us start with the 3-point function with only tensors. From the cubic Lagrangian (5.6), we have that the relevant vertex coupling reads Thus for h p k 1 h q k 2 h r k 3 , by using (C. 15) and (C.16), we have where and δ p,q,r = 1 3 (δ pq + δ pr + δ rq ) .

T T S
Consider the 3-point function of the form ζ n k 1 h r k 2 h s k 3 . In the democratic approach, two cubic vertices are dominant: h h π L and h h π 0 ; see (5.5). The computation of the bispectrum is done by using (C.12) and (C.13).
• For the h 2 π L vertex we get where we have rescaled the momentum of the scalar field and defined the total rescaled momentum k T κ 1 = c s a k 1 , k T = κ 1 + k 2 + k 3 , a = 1, 2 ; (8.11) and .
• In a similar way, one can compute the contribution from the vertex h 2 π 0 : where In the squeezed limit, when the scalar momentum k 1 (see eq. (8.3)) is squeezed, the dominant contribution is the one from the h 2 π 0 and has a characteristic blue operator behavior with the presence of a single scalar cross-correlation, while the h 2 π L contribution is proportional to P ζn ; namely Let us define f NL for the TTS bispectrum as By squeezing along the scalar direction and keeping only the dominant blue vertex, we get the following expression In figure 10b, we show the logarithmic plot of the dominant Abs [f T T S ]. In this extreme case, setting B (tts) to one, the f T T S can get as large as 10 4÷5 for c s2 ∼ 0.1.

T SS
Finally let us analyse the 3-point function ζ n k 1 ζ n k 2 h k 3 .
As for the TTS case we can define a sort of f NL -like parameter in order to understand what is the typical strength of the TSS Bispectrum.
Thus, with the f T SS symbol we mean where k 3 is the momentum related to the tensor field h, while k 2 is one of the momentum related with the scalar field ζ n as presented in eq. (8.4). Consider a typical interacting TSS term where the generic scalar fields ξ i , ξ j stand for any among ζ n , R π 0 and their time derivatives. Thus, the TSS i − j interaction also selects the possible color related to the specific vertex 14 . Manipulating eq. (C.10) one can demonstrate that such a parameter can be described by a very compact formula where r (l) The scale dependent factor S describes the structure in polarization and spatial derivatives provided by the particular interaction For all our interactions, the S function does not depend on the polarization s, thus the f NL -parameter defined in eq. (8.19) does not depend on s and coincides with the average on the two polarizations. Finally, a shape integral needs to be defined: where the shape f i functions are given by .

19) is valid if and only if
S is symmetric under exchange of the scalar momenta, and this is always true in the squeezed, equilateral and folded shapes. In particular for the squeezed limit k 1 → k T ≡ k L , k 2 → k ζ ≡ k S and k 3 → k S + k L cos(θ) We now need to compute five TSS vertices contributions: h π 2 L , h π L π 0 , h π 2 L and h π 2 0 . Looking at vertices spatial derivative structure and contractions with polarization tensor, one can easily get that the folded limit always vanishes, i. e. S(k, k 2 , k 2 ) ≡ 0. Defining the following 2×2 symmetric matrices we get the shapes defined in table 4, while the equilateral formula are very similar and reported in Appendix E. Note that each vertex is classified with our color-code, starting from violet (asymptotic f T SS of order c 0 s2 ) arriving to green vertices (asymptotic f T SS of order c −5 s2 ). In tables 5 and 6 the reader can find the asymptotic c s2 1 of squeezed and equilateral f T SS .    Table 6: Equilateral f T SS for c s2 1. We set A Log = (3 γ e + 3 log(−k t) − 4) .

Solid vs Supersolid
Before giving an overview of primordial non-Gaussianity for a supersolid, it is interesting to briefly recap the main differences with solid inflation introduced in [10]. At the linear level, in the instantaneous reheating scenario, the adiabatic power spectrum transmitted to the radiation phase is pretty similar 15 . Indeed, we get [9] Solid : P Φ | Rad ∼ P ζn ≡ P 1 (4.20). Isocurvature perturbations are very suppressed (of order ζ n ). The crucial point is the fact that instantaneous reheating for a supersolid effectively almost filters out non-adiabatic modes. Isocurvature modes have a much bigger amplitude but are mostly dissipated during the transition. Such a modes, particularly the ones related to π 0 , are the main sources of the non linearities during inflation, affecting the scalar and tensor power spectra and bispectra. In solid inflation there is a single scalar mode π L , no isocurvature perturbation is present and all non linearities are due to the single mode whose amplitude is small as dictated by observed scalar power spectrum. The main difference between the two models resides in the presence in a supersolid of two independent scalar modes π L and π 0 with intrinsic mixing in the quadratic Lagrangian. Upon a non-trivial Hamiltonian diagonalization procedure, one gets the power spectra of ζ n and R π 0 (trivially related to π L and π 0 in the flat gauge). While, given the reheating adopted, the power spectrum of ζ n is strongly constrained by observations, this is not the case for the one of R π 0 that can be significantly enhanced when c s2 is small triggering sizeable differences with solid inflation at the level of non-Gaussianity and power spectra non linear corrections. The different non linear structure of solid and supersolids implies that not only the size of non-Gaussianity is different but also its form.
• The linear PS for tensor (gravitational waves) gets important corrections from the secondary production that can become of the same order or even dominant. As a consequence of the presence of π 0 in a supersolid not only can enhance the amplitude of the GW produced during inflation, but induces a blue tilt of the PS. A blue tilt makes the PS to grow as k n T −1 for large k values and potentially can enter in the LISA, DECIGO and ET sensitivity region. The secondary production originates from the pure supersolid interaction G (tss) h ∂R π 0 ∂R π 0 of green type. In [9], the correction to the linear tensor PS was estimated to be While for solids the tilt tends to be suppressed by the slow-roll parameter , for the supersolid we get n T − 1 = 12 (1 + c 2 b ) + 12 c 2 b + 2 η . (1 + c 2 b ) we naturally get a conspicuous deviation from scale invariance toward the blue.
• For bispectra in the scalar sector there are important differences between solid and supersolid inflation. Consider a squeezed configuration. While in solid inflation, if the special approach used in [10] is considered, a pure quadrupole term is present, in a supersolid both a monopole and quadrupole term are found. The above can be interpreted as the indissoluble presence of a fluid component (responsible for the monopole term) and a solid component (that induces a quadrupole term). In addition, a monopole term is induced by the supersolid operators w X/Z and a quadrupole term by superfluid operator χ. It is worth to stress that the size of non linearities in supersolids can be enhanced up a factor 1/c 8 s2 for small c s2 .
• For a squeezed shape in the bispectra of the TSS sector, we get the same angular structure for both solids and supersolids For a supersolid, the amplitude can be enhanced by a factor 1/c 5 s2 for small c s2 .
• Finally, for the bispectra (squeezed shape) in the sector TTS both a monopole and a quadrupole are present; namely While in solid inflation their size is similar, for a supersolid the monopole gets enhanced by a factor 1/c 3 s2 for small c s2 . This leads a characteristic signature of supersolid inflation: the non-linear and local corrections to the tensor PS that becomes rather different from solid inflation, as discussed in [36] where TTS predictions for the standard solid inflation [10] scenario are compared with another "isotropic" (in the sense that the f T T S gives a monopole) model [37].

PNG scenarios in supersolid inflation
As it is evident from the previous analyses, the amount of non-gaussianity generated by supersolid inflation can be important. As already discussed, the presence of entropic perturbations with a sizeable amplitude (during inflation) generates not only noticeable non-gaussianity but also important non linear effects as secondary gravitational waves production. Clearly such enhancements can become leading for observables that start tiny at tree level, as it is the case for the gravitational sector. Related to this problem, as already discussed previously, it will be interesting to check the amount of non linear corrections that could in principle enhance or leave a peculiar angular and scale dependence to the primordial tensor PS. In particular the presence of important corrections to f T T S are expected to be related to sizeable non linear and local corrections to the tensor PS. Similarly, important corrections to f SSS and f T SS are expected to non linear correct the scalar PS. Such a program is under study. So, let us come back to our "tree level" evaluation of the bispectra. Recalling that stability of the model gives the constraint c s1 <ĉ L < c s2 (we fixed for convenienceĉ L = 1/2), as one can infer from the figures in section H, we can discriminate two main regions in the (c s1 , c s2 ) space parameters. An intermediate region where c s2 is not too small if compared to c s1 (c s2 c s1 ), and an other region where one of the two diagonal sound speeds is very small, c s2 c s1 . In the first region (c s2 c s1 ), almost all the cubic interactions classified in this work equally participate in defining the total f NL -parameters for SSS, TTS and TSS interactions. In this case, some of the main sections' results need to be substituted with their total expressions given in Appendixes D and E. In particular, for the purely scalar sector, we have a summation of several terms that have to satisfy stringent constraints in the squeezed limit [34], and this can be easily achieved by adjusting our free parameters. At the same time, setting the coloured X-couplings to one as shown in figures (10b, 11,12,13) the TTS and TSS f NL can easily achieve values of order 10 2 for c s2 around 0.3 making supersolid inflation an interesting model. A more quantitative and accurate analysis of this region is left for detailed future work. The main motivation to consider the region c s2 c s1 is the possibility to enhance via secondary production the gravitational waves boost during inflation driven by a supersolid [9]. Following the same line of [9] taking c s2 1 we will argue soon how to maximize the mixed scalar-tensor f NL parameters. Let us start in the scalar sector for which f has a minimum in the region c s1 → 1 and c s2 → 1/2, while the maximum is attained for c s1 → 1/2 and c s2 → 1/2.
• The contribution of the remaining operators has only a minimum in the region c s1 → 1 and c s2 → 1/2 and diverges in the limit c s2 1.
Once we impose that |f (SQ) NL | < 10 for each vertex contribution, the results are the following absolute bounds for the couplings: When c s2 ≤ 0.1 the constraints are much stronger: where we normalized c s2 asĉ s2 ≡ c s2 0.
An interesting way to classify the vertices is to divide them into two main classes: operators compatible with a fluid or superfluid medium and operators compatible with a solid or supersolid medium, as shown in table 7. Table 7: Operators in the cubic SSS Lagrangian classified medium-wise.
The operators b, y, χ pertain to a fluid or a superfluid, τ Y /Z characterize solids while w Y /Z supersolids, see [9]. The couplings in L 3 are derivatives of the Lagrangian U with respect to the operators (2.5) and can also be divided according to their presence in the cubic Lagrangian of a fluid/superfluid and or of a solid/supersolid. Last but not least, only the coupling λ 6 depends on second derivatives with respect to operators characterizing both fluid/superfluids and solid/supersolids. As a result, a generic vertex can be ascribed to a specific class of medium. In the table 7 it is shown the contribution of a vertex to f is set to be ≤ 10 and c s2 ∼ 0.1. The result is that R, O, and G operators are rather constrained and very sensitive to the "portion" of fluid/superfluid of the medium. Thus the result is clear: if we want to enter a region with small enough c s 2 values, we have to mitigate the squeezed red, orange vertices typical of the fluid/superfluid component of the medium to avoid too large contribution to f (SQ) NL 16 . In this way, we can have the non-gaussianity in the scalar sector completely under control and entirely constrained by PLANCK results [34], and at the same time to be able to get a sizeable secondary tensor production (of the same order of the primary one or even dominant) as described in [9] and to maximize the TSS/ TTS primordial non-gaussianity.
Let us give one among several descriptive examples. Consider λ 2,6,10 free parameters as slightly c s2 -dependent. The remaining λ 9 and λ 7 as purely constant parameters. These last ones are defined by w-operators, playing a crucial role in the tensor sector. By setting the λ 2,6,10 parameters as in figure 2 (explicit expressions are given in Appendix G), and expanding the total f NL squeezed in c s2 we get f (SQ) For instance, in this particular example (see Appendix G) we reproduce a solid-like SSS scenario where f Q is set to an optimal value around 10 for small c s2 values as it can be inferred from figure 3. At the same time, with the same parameters, we verified that the total equilateral and folded f NL assume values between 0 to ±50 in a large range of the c s2 parameter as shown in figure 4, and the situation is almost unchanged varying c s1 17 . In order to not have an excessively high equilateral/folded shape, we are forced to take λ 9 between 0.1 and 0.01.  With this particular choice, we have the possibility to access to an interestingly small-c s2 region whose effect is to enhance the f T T S and the f T SS parameters as it can be deduced by figures 5.  Figure 5: On the left the dominant h R 2 π 0 f NL -like parameter with sin(θ) 2 set to one and λ 9 = 0.05. On the right the dominant h 2 R π 0 f NL -like parameter with set to 1, and k S t ∼ 10 −5 . For instance, considering ∼ 10 −2 we get f T T S of order 10 3 /10 4 in the region c s2 ∼ 0.1. 17 We used the λ parameters appearing in red and orange monopole to compensate the squeezed green vertex. Consequently, the equilateral/folded shape is of green type (O(c −6 s2 )) at maximum.
Note how, taking λ 9 relatively small, the dominant f T SS contribution is always between 10 ÷ 10 3 in the region (c s2 ∼ 0.1 − 0.2, c s1 ∼ 1) in agreement with the current WMAP constraints [38], while the f T T S parameter is even more enhanced thanks to the λ 6 parameters (∼ 10 5 if λ 9 ∼ 10 −2 ÷ −1 ). Finally, the TTT f NL is well within the PLANCK constraints, being the V T parameter almost untouched by this parameterization and defined by λ 7 and c L only (see (8.5)).

Conclusions
A considerable effort is devoted to search for signatures of inflation by studying the CMB and the large scale structure of the Universe. The progress in detecting GWs and the forthcoming space interferometer LISA give new opportunities to test the physics of inflation, and more in general, the early Universe on scales very different from the ones of the CMB [39][40][41].
In this context, primordial non-gaussianity can be used to distinguish single clock inflation from multi-field models. Specifically, in single-clock inflation the 4D-diffeomorphism of GR are broken down to the 3D-diffeomorphisms of the constant inflaton field hypersurface and non-gaussianity is tiny [3], the production of gravitational waves tends to be suppressed for modes that re-enter the horizon during radiation domination 18 , making almost impossible the direct detection by LISA [42][43][44][45]. We have studied the physical impact of a change of the symmetry breaking pattern during inflation. By using effective field theory based on four Stückelberg-like scalar fields, the complete breaking of diffeomorphisms down to a global group can be described in terms of a self-gravitating supersolid computing f NL for squeezed and equilateral configurations extending the results of [9]. Let us schematically highlight the main differences among single-field inflation, solid inflation (democratic approach) and our supersolid inflation model.
1. Symmetry breaking pattern: broken time reparametrization in single field inflation (one scalar DoF), broken spatial diffs. in solid inflation (one scalar and two vector DoFs), completely broken time and space diffs. in supersolid inflation (two scalar and two vector DoFs).

Perturbations
• Single field inflation: ζ and R equivalent at superhorizon scales and reheating independent.
• Solid inflation: ζ and R not equivalent at superhorizon scales, weakly time dependent a − 4 3 c 2 2 and reheating dependent. 18 This is a consequence of the red-tilted spectral index for tensors in single-field inflation. 19 The presence of the term a 4 M 2 pl H 2 c 2 2 h 2 ij in the gravitational action (4.2), becomes a "mass term" once we reduce to Minkowski space (this is the origin of the denomination "mass"). In dS space its presence introduces a mild time dependence (a − 8 3 c 2 2 ) for the superhorizon power spectrum of the gravitational waves.
• Supersolid inflation: two independent scalar perturbations R π 0 and ζ n . The former is the analog of the comoving curvature perturbation in single field inflation, while ζ n is the analog of ζ in solid inflation. In the instantaneous reheating approximation, ζ n determines the primordial seed of adiabatic perturbations, while R π 0 is mostly related to entropic perturbations.
4. Bispectrum of ζ n in a squeezed configuration.
• Single field inflation: a tiny monopole term, f SQ NL ∼ O( ). • Solid inflation: a sizeable quadrupole term, f SQ NL ∼ O(1) Y 0 2 (θ) • Supersolid inflation: sizeable monopole and quadrupole terms, It is evident that the violation of the consistency relations for solid and supersolid inflation is related both to the magnitude and the angular structure. While in solid inflation, an angular average restores the consistency relation [46], this does not happen in supersolid inflation. As it is evident, the main difference between supersolid and the other two alternative models (single filed and solid inflation) is the presence of two independent modes with phonon-like dispersion relations and a kinetic mixing present throughout the inflationary period. By studying the relevant 3-point functions, we have shown that the form and size of primordial non-gaussianity are very different from single-clock inflationary models. The most general supersolid produces a considerable enhancement for all sectors (SSS, TTS, and TSS) 21 in the small c s2 limit. Thus we have essentially two main possibilities.
• Do not tune the free parameters. In this case, we cannot consider values of c s2 too small, effectively reducing the allowed region of the parameters c s1 and c s2 ). Contrary to solid inflation, monopole and quadrupole f SQ NL terms arise, violating the Maldacena consistency relation in the SSS,TTS and TSS sectors. Getting the enhancement of GWs during inflation is more difficult.
• We have enough free-parameters, to mitigate the too-large contribution to f NL of SSS green ∝ O(c −5 s2 ), orange ∝ O(c −6 s2 ) and red ∝ O(c −8 s2 ) vertices in the small c s2 limit. In this case, we can have an interesting enhancement of the TTS and TSS sector; indeed the presence of the w-operators, specific of a supersolid still gives The former could leave a signature in LSS next-generation experiments, giving a specific imprint 22 on the distribution of galaxies [33,48,49], while the latter could play a crucial role in the generation of anisotropies of the gravitational waves background [36,37,50,51].
We also mention that in this paper we suppose a suitable number of e-folds needed to solve the horizon problem in order to avoid the issue with statistical anisotropy discussed in [17].
Regardless of the choice made, supersolid inflation is a very interesting model whose features could be tested in next generation experiments probing the early Universe.

Note Added
After the completion of the present work the preprint [52] was posted with some overlap with our results.
over momenta and the delta functions enforcing momentum conservation in each vertex is understood. Thus, the cubic Lagrangian in the scalar sector can be written as where k = | k| andk = k k . The values of dimensionless couplings of the vertices are given in the following table.

Vertex
Value

A.2 TSS Sector
In the TSS sector we have The expressions for the vertices are given in the following table

A.3 TTS Sector
For the TTS sector we have the vertices are only of violet and green type given following table.

Vertex
Value

B Entropy Term and Total Derivatives
In several points of the paper we stated thatσ appears only as a coefficient multiplying a total derivative term. Let us show how it comes about. Let us consider the operators of interest as a background (3.2) plus a perturbation (all orders) b =b + δb, Y =Ȳ + δY, χ =χ + δχ ; (B.1) and let us expand the Lagrangian for the scalar fields in (2.4) in the following way Now, from this expansion we can identify the derivatives of U that containsσ. To do that, we have to identify where theσ appears. For Thus, we see that U Y =bσ + · · · , U b =Ȳσ + · · · and finally U bχ =σ + · · · . So the coefficient δU |σ in the expansion of U that multipliesσ is given by δU |σ = (Ȳ δb +b δY + δb δχ) .

(B.4)
Note that when we implement the constraints in (4.5), we have to impose U Y =bσ + ..., U b =Ȳσ + ..., U bχ = 0 + ... and U bY =σ + ... that implies δU |σ = (Ȳ δb +b δY + δb δY ) =σ (Ȳ δb +b δY + δb δY ) = δ(b Y ) (B.5) Last step is to show that the term b Y form a total derivative. With the field content of our model it is easy to build the following total derivative term [53]: Thus, at cubic order, we have that all the terms in the expansion of the Lagrangian for the scalar fields in (2.4) proportional toσ can be collected in L σ given by It is well known that spatial derivative terms never contribute to the primordial bispectrum. Surprisingly, in general this is not the case of a time derivative term [54][55][56][57]. We can write the cubic total derivative expression proportional toσ as Using perturbation theory, such a local term gives a contribution local in time of the form Such terms are commonly estimated considering a quadratic field redefinition and absorbing the time derivative terms in the cubic Lagrangian [3,47,58]. In our particular case, Q does depend by π 0 and π L conjugate momenta and then it commutes with ζ n ∼ ∇ 2 π L and the contribution to the bispectrum of ζ n (B.9) is zero. Notice that, it is the choice of an instantaneous reheating that selects ζ n as the reference scalar perturbation that sets the initial conditions for the radiation domination phase of the Universe. Quite different scenarios can be considered computing ζ, or R bispectra where the conjugate momenta appear [9]. As final remark we point out that the vertices V 4/5 and G 1/2 do not contribute to (B.6); indeedσ cancels due to the condition M 1 ≈ 0.

C PNG and the In-In Formalism
In this section a general description of the in-in formalism and the main formula to get our results are given. The in-in formalism is based on the interacting picture, each quantum field 23ξ evolves following dξ dt = −i ξ ,Ĥ 2 , (C.1) whereĤ 2 is the quadratic(free) Hamiltonian. The leading contribution to the 3-point function of fξ 3 operator, by expanding the evolution operator, will be given by By definition of the interacting Hamiltonian, the reader can easily verify that supersolid inflation can be included in the class of theories such that whereL (3) is the cubic Lagrangian density. Thus, in Fourier space the bispectrum reads (C.4) 23 To make clear that we are dealing with quantum fields anˆis used.
where L q 1 q 2 q 3 is the Fourier transform of L (3) . Equation (C.4) is completely general and can be specialized to the case scalar, tensor or mixed case.

SSS:
vertex of the form L k, k k = X M 2 pl H 2 a n D(k, k , k ) ξ ωk ξ δ k ξ ρk , (C.5) where ξ Greek letter can be any scalar field among ζ n , ζ n , H −1 R π 0 , H −1 R π 0 . For the bispectrum of ζ n we have where δ l,i,...,j m,n,...,0 = δ l m δ i n ... δ j o . The dimensionless function D encodes the structure of spatial derivatives in the vertex structure, while the n and X are vertex dependent constants.

TSS:
vertex of the form L k, k k = X M 2 pl H 2 a n D ij (k, k , k )h ij k ξ ωk ξ δ k , (C.8) which leads to B ζ 2 n h p = dk 3 dk 3 dk 3 I , (C.10)

TTS:
vertex of the form L k, k k = X M 2 pl H 2 a n D lm (k, k , k ) h im k h il k ξ ωk , (C.11) which leads to: (C.13)
(C. 16) The results presented can be obtained by evaluating the above integrals. In the cases c 2 b = −1, 0, the integral J can be computed analytically, while when c 2 b ∈ (0, 1) a numerical approach is needed.

D Squeezed SSS f NL
In this appendix we give the full expressions for f NL in a squeezed configuration. Among ten vertices, only four give an independent contribution to f NL . The contribution from each independent operator O of (A.3) to f NL can be split in an overall amplitude M times a momentum-independent part wich gives the monopole and the quadrupole structure according with It is useful to define the following recurrent combinations in f NL of the diagonal sound speeds In the case of c 2 b = −1, at the leading order in SR, ζ n and R π 0 fields have analytic solutions of the form (−k t c sl ) , l = 1, 2 ; (D. 5) used to evaluate the in-in time integrals. We choose the following four independent vertices for their simplicity (D.9) • ζ n ζ n R π 0 vertex (D.11) ζn + r (D.17) L/0 + r   For the independent parameters in the cubic Lagrangian we get λ 2 = UΛ , X X X + 3 UΛ , X YY + 3 UΛ , X X Y + UΛ , YYY , λ 6 = −∂ X c 2 2 − ∂ Y c 2 2 , λ 7 = 1 9 UΛ , τ Z + UΛ , w Z , λ 9 = 1 9 UΛ , w Y + UΛ , w Z , Thus, (F.1) indeed satisfies all the requirements for a consistent slow-roll dynamics.

G Boosting the tensor sector
In this appendix we discuss the parameterization used in section 9 to keep the SSS sector under control. Figure 2 was obtained by the following choice for the five independent parameters {λ 2 , λ 6 , λ 7 , λ 9 , λ 10 } in the cubic Lagrangian NL | ≤ 60 in order to get an almost violet squeezed f NL without exceeding the equilateral/folded experimental constraints, see figure 6. This is only one possibility to get the squeezed f NL of order c 0 s2 and within the experimental constraints obtaining f (SQ) NL ≈ 2.3 + 9.5 cos(2 θ) − (6.9 + 5.5 cos(2 θ))ĉ s2 + −11.4ĉ 2 s2 + O(ĉ 3 s2 ) , (G.4) whose total plot is shown in figure 3, and allowing us to enhance the tensor PNG having access to a sufficiently small c s2 values.

H Bispectrum figures
In Thus, for each f NL we plot M functions by setting (X 1 + X 2 cos 2 θ) to one, and c 2 b = −1.  Figure 12: Squeezed TSS plots. Note that the shape h π 2 0 is the dominant one. We set V-B constants and sin(θ) 2 set to one.