Boosting GWs in Supersolid Inflation

Inflation driven by a generic self-gravitating medium is an interesting alternative to study the impact of spontaneous spacetime symmetry breaking during a quasi de-Sitter phase, in particular the 4-dimensional diffeomorphism invariance of GR is spontaneously broken down to $ISO(3)$. The effective description is based on four scalar fields which describes the excitation of a supersolid. There are two phonon-like propagating scalar degrees of freedom that mix non-trivially both at early and late times after exiting the horizon giving rise to non-trivial correlations among the different scalar power spectra. The non-linear structure of the theory allows a secondary gravitational waves production during inflation, efficient enough to saturate the present experimental bound.


Introduction
Inflation is the most compelling way to solve the drawbacks of the hot big bang and simultaneously generate the seed of the primordial perturbations to be used as initial conditions for latter stages of post inflationary evolution. The simplest class of models is single clock inflation, where time diffeomorphisms is non-linearly realized, whose predictions are largely independent on how the Universe is reheated. Indeed, according to the Weinberg theorem on cosmological perturbations [1,2], at large scales and under mild assumptions, the curvature perturbations of the constant density hypersurfaces ζ, or equivalently the comoving curvature R, are conserved and can be used to set the primordial initial conditions for the scalar sector at the beginning of radiation domination. The situation is different for models characterized by different symmetry breaking patterns, featuring more degrees of freedom and Weinberg theorem does not hold anymore: neither R nor ζ are conserved and equivalent at superhorizon scales and the details of reheating have to be taken into account [3,4,5,6,7]. That is exactly what happens when a fluid [8] or solid [9] drives inflation. In this work we present an effective field theory (EFT) description suitable to describe the complete breaking of spacetime diffeomorphisms during inflation by using the minimal set of four scalar fields {ϕ A , A = 0, 1, 2, 3} sporting a suitable set of internal symmetries. As matter of fact, ϕ A can also be interpreted as the coordinates of a self-gravitating non-dissipative medium [10,11,12,13,14,15] that in our case is a supersolid. Besides a complete analysis of linear evolution of scalar and tensor modes with the computation of the corresponding power spectra, we consider the secondary production of gravitational waves during inflation, exploiting the cubic tensor-scalar-scalar vertex of the theory which allows to saturate the experimental bound set by Planck without upsetting the scalar 3-point function. A detaled analysis of non-Gaussinity can be found in a companion paper [16].
The outline of the paper is the following. In section 2 we briefly review the effective field theory description of a supersolid at the leading order in derivates. In section 3, the dynamics of the two independent scalar perturbations is careful analyzed both at classical and quantum level, computing the relevant scalar power spectra and constraining the parameter region by using Plank data. Section 4 is devoted to study, in the instantaneous reheating approximation, how the seed of primordial perturbations are transmitted to the radiation dominated phase in a ΛCDM universe. In section 6 primary and secondary gravitational waves production during inflation is considered. Our conclusions are drawn in section 7.

Super Solids and Inflation
Several features of inflationary models can be traced back to the spontaneous symmetry breaking pattern: in single field inflation, the non-trivial time dependent configuration of the inflaton breaks time reparametrization leaving unbroken the space diffeomorphisms of the t =const. hypersurface. However, there are other possibilities, for instance an inflationary model where spatial diffeomorphisms are non-linearly realized was studied in [9] by working with three scalar fields. In a similar fashion, one can consider a more general case in which all diffeomorphisms are broken by the background configuration of four scalar fields which will be the background configuration for the inflationary phase. The existence of a spatially homogeneous background is allowed by the presence of global symmetries of the scalar field action. Consider a special multi-field model of inflation based on four scalar fields ϕ A , A = 0 , 1 , 2 , 3 with shift symmetry [11,15] and SO(3) intermal symmetry ϕ a → ϕ a = R a b ϕ a , a = 1, 2, 3 R t R = 1 . (2. 3) The "vacuum" configuration (2.1) has a residual global "diagonal" ISO(3) symmetry. Indeed a global spatial rotation R i j x j can be absorbed by a corresponding inverse internal transformation of ϕ a and the same is true for a global translation x i → x i + c i thanks to the shift symmetry (2.2). Among the shift symmetric operators with a single derivative of ϕ A (2.5) where u µ plays the role of a normalized four velocity such that u µ ∂ µ ϕ a = 0 and B ab = C ab , W ab = B ab − C 0a C 0b C 00 , a = 1, 2, 3 . (2.7) By using the relation ς = χ −1 b y and the Cayley-Hamilton theorem, only 7 among those operators are independent. Thus, we arrive to the action which can be interpreted as the relativistic generalization of the low-energy effective Lagrangian describing homogeneous and isotropic supersolids at zero-temperature [17,18]. Such an action is the most general at leading order in a derivative expansion compatible with (2.2) and (2.3) and is rather useful to study systematically the symmetry breaking pattern of spacetime symmetry during inflation. By inspection of the energy momentum tensor (EMT), the energy density and the pressure are given by According to the Noether theorem there are four conserved current three associated to solid configurations that spontaneously break translation invariance and one associated to the frictionless flow of the superfluid. In particular, the particle number density n sf of the superfluid component can be expressed in terms of the Noether current J 0 µ as This allows us to define the superfluid density per lattice site σ as (2.14) As we will see, at cosmological level, the perturbations δσ generate non-adiabatic contributions (for this reason, we will regard δσ in the following as isocurvature or entropic perturbations). Similarly, for the rest of the paper we identify n (2.13) with the usual particle density n . While u µ represents the 4-velocity of the normal component of the supersolid, is the 4-velocity (irrotational) of the superfluid component. One of the key features is that two independent phonon-like excitation are present; in general, the supersolid perturbation can be written around flat space At the quadratic level we have [19] whereρ andp are the background values of the energy density and pressure (constant in space and time) while ∆ = δ ij ∂ i ∂ j ; finally the conformal time derivative of a function f is denoted by f . The parameterM a = M 2 Pl M a are proportional to first and second derivatives of U and are given in appendix A. Notice that the space shift symmetry is crucial to have a homogeneous EMT even if the scalar fields have non-trivial background values. The properties of the EMT are largely determined by the symmetries of U as discussed 2 in [14,15]. It is useful to summarize the main features associated with the presence or absence of some of the operators (and related mass parameters) in the Lagrangian which corresponds to specific internal symmetries: • Perfect Fluids: -U (b) : only {ϕ a , a = 1, 2, 3} are present; the Lagrangian is invariant under internal volume preserving diffeomorphisms V s Diff: ϕ a → Ψ a (ϕ b ) , det |∂Ψ a /∂ϕ b | = 1, a, b = 1, 2, 3. -U (χ) : it is the most general Lagrangian for a perfect irrotational fluid with ϕ 0 only. • Superfluids U (b, y, χ) : invariant under transformations of ϕ a corresponding to V s Diff and ϕ 0 → A more detailed analysis about thermodynamical properties for general supersolids is planned for a future work.
Stability of (2.17) imposes the following conditions [19] As we will see such conditions are necessary for existence of the Bunch-Davis (BD) vacuum in a inflation phase driven by a supersolid. During a quasi deSitter period the most convenient parametrization of the mass term is through some c 2 i parameters such that (where w is defined in (3.5)) with the assumption that c 2 i are slowly varying in time ( 2 In [14,15] the set chosen independent operators is different from our choice without changing the physics.

Slow Roll
Cosmological perturbations in the flat-slice gauge are described by Perturbations in a generic gauge are discussed in Appendix B. The background EMT tensor has the form of the one of a perfect fluid with energy density and pressure given by (2.9) and (2.10) evaluated on FLRW; EMT tensor conservation, at the background level, is equivalent to 3 where Our benchmark values for c 2 b will be c 2 b = 0 which givesφ = aφ 0 and c 2 b = −1 leading toφ =φ 0 a −4 . Inflation takes place when We will be mainly interested in SR inflation 4 for which the following dynamical parameters are small Note that in a quasi dS phase, the adiabatic speed of sound is given by Both Ψ and F are non-dynamical fields and their linear equations of motion can be solved in terms of π 0 and π L , at the leading order in SR and working in Fourier basis, one finds The leading order in linearized dynamics is described by the following action In [20] it was 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. Moreover the correct implementation of the Stuckelberg trick at the background level requires a non-trivial background for ϕ 0 satisfying (3.2). 4 As discussed in [21] super SR is also possible; actually when M2 = 0, as for fluids, this is the only viable regime with small . (3.12) Up to boundary terms one can always take D antisymmetric. The peculiarity of (3.9) is the mixing in between the two propagating DoF (degrees of freedom) present at kinetic level due to the matrix D and at mass level being the matrix M non-diagonal. Such a mixing is unavoidable unless the parameters c b and c 2 2 , c 2 1 are unnaturally tuned and is a key property of a superfluid component in the solid which will be the origin of cross-correlations in the two and three points function of any scalar perturbation. As a result, the study of scalar linear perturbations is a bit involved and to get rid of the mixing by a suitable field redefinition few steps are needed. A similar system of coupled modes, described by (3.9), was encountered when studying the non-thermal production of gravitinos [22], multi-field inflation [23], chromo-natural inflation [24,25] and in effective theories of inflation [20]. As far as we know, our analysis is the first complete one that does not rely on special choices of parameters. The strategy to quantize the quadratic action will be the following. We start from the original fields π 0,L that describe physically two Nambu-Goldstone modes around a non-trivial Lorentz breaking background solution. The quadratic action controlling the dynamics of such modes (3.9) exhibits both kinetic (the presence of D) and mass mixing effects (non-diagonal M). A similar kinetic mixing is also encountered in mechanical systems with gyroscopic forces like the Coriolis force or in presence of magnetic fields; it is worth to stress that the the D mixing can take place when at least two fields are present. The first step is to make the fields canonical by a time-dependent field redefinition Π = K 1/2 π (3.17). At this level the corresponding Lagrangian L(Π, Π ) (3.18) is characterized by non trivial D − mixing and a time dependent non diagonal mass matrix. The classical equations of motion correspond to a coupled system of second order equations or, alternatively, to a two decoupled fourth order differential equations. The quantization of the system goes through the choice of the Bunch Davis (BD) vacuum i.e. the study of the Lagrangian in the UV (k → ∞) regime (C.1) where k dominates over all other scales. In this regime the mass term is diagonal and time independent. Thanks to this feature we can write a decoupled system of quantum oscillators and quantize it with the usual canonical rules in

Quantization and Power Spectra
In this section we compute the Power Spectra (PS) amplitudes and later we analyze the next to leading order correction in the SR expansion. Before proceeding let us fix the notation in order to define the PS of a general quantum scalar field ξ(x) with two DoF. In Fourier space we set where the (cl) subscript stands for classic solution, and with the Latin index i = 1, 2 we discriminate between two different DoF whose related annihilation and creation operators obey the standard canonical commutation relations (3.14) Thus, the ξ two point-function reads 15) and the common scale-invariant PS is defined as usual (3.16) The first step to compute quantum correlators during inflation is to introduce the canonical field Π defined as The conditions (2.17) guarantee that the matrix K is positive definite. Given that the matrix elements of K are time-dependent, besides turning the kinetic term into a canonical one, the transformation (3.17) also alters the form of D and M; thus, the quadratic Lagrangian in (3.9), when written in function of the new canonical fields becomes and we have defined We have kept only the leading contribution in the SR parameters, obtaining the following equations of motion In order to quantize (3.18) we need to remove the kinetic mixing introduced by D. Our strategy is the following: in the UV, at very large k, M becomes diagonal and time-independent. Thus, at very large k, the original Lagrangian (3.18) is equivalent to L (UV) 2 , in accordance with the equivalence principle. In that regime, by using a canonical transformation, one can reduce the Hamiltonian associated to L (UV) 2 to a system of two canonical free fieldsΠ = (Π L ,Π 0 ) linearly related to Π Thus, the unique Fock space vacuum is the Bunch-Davies vacuum for the system. Details can be found in Appendix C. The existence of the Bunch-Davies vacuum requires the frequencies squared , ω 2 2 = k 2 c 2 s2 to be strictly positive or equivalently c 2 si > 0 , i = 1, 2. In addition we restrict ourselves to the case of subluminal "diagonal" sound speeds: 0 < c 2 si < 1. The condition (2.17) are sufficient conditions for c 2 si > 0 and when expressed in terms of (2.18) gives 5 We have checked that there is a large region of parmeters c i where the stability conditions hold together with c 2 s1 , c 2 s2 < 1; moreover in such a region, the first two conditions when rewritten in terms of c s1 and c s2 , become where we choose the convention c s2 < c s1 . The equations of motion (3.23) constitute a coupled linear system of two second order differential equations with time-dependent coefficients. Finding explicit solution is not an easy task; of course, one could solve the equations numerically. However, from a physical point of view, it is more transparent to quantize the system focusing on the following values of c b : c 2 b = −1 and c 2 b = 0 which have a special interest and for which an analytic solution can be found. Neglecting SR corrections, the coupled system of second order equations can be written as a fourth order equations for both Π 0,L fields. Remarkably c 2 b = 0 and c 2 b = −1 gives identical equations (that are valid as soon as c 2 in the evolution equations. The solutions can be written as a linear combination of Bessel functions of order 5/2 and 3/2; the integration constants are fixed by imposing that subhorizon, where x = −k t 1, the solution matches the one given in (C.13) and (C.14) which represents flat space modes; such a choice is equivalent to select the Bunch-Davies vacuum. Thus, as quantized free (Gaussian) fields we have that Π L and Π 0 are given by (3.28) the expression for C k are the creation operators for the fields defined in eq.(C.11). In single field SR inflation, naturally two gauge invariant scalar quantities can be considered when studying the dynamics of superhorizon modes: the curvature ζ of the ρ = const. hypersurface and the curvature R of hypersurface orthogonal to the scalar component of the fluid 3-velocity According to the Weinberg theorem [1], in standard single field inflation, both ζ and R are conserved and coincides at superhorizon scales; as a result the power spectrum of primordial perturbations during inflation is given in terms of the Fourier transform of the 2-point function of ζ or equivalently of R. In our case, the Weinberg theorem does not hold and, besides the above curvature perturbations, additional gauge invariant scalar perturbations can be considered. In particular the curvature ζ n of constant particle number n-hypersurface (remember that n = n , (2.13)) and curvature of the ϕ 0 =const. hypersurface; namely Note how R π0 is the comoving curvature related to the superfluid component (2.15). In Appendix B the reader can find the expression of the various curvature perturbations in terms π L and π 0 in a generic gauge. In the spatially-flat gauge (3.1) we have that As we will see soon, one can deduce that if −1 ≤ c 2 b ≤ 0 both R π0 and ζ n will have almost scale-free power spectra 6 , namely n s = 1. The uniform curvature perturbation ζ can be obtained from eq. (B.15) at leading order in SR and for superhorizon scales and similarly, from eq. (B.14), for the comoving curvature As anticipated, by using (3.28), being C (j) , the power spectra of ζ n , R π0 , ζ and R will be scale-free, up to SR corrections. Moreover, as shown in (B.14) and (B.15), ζ and R are linear combinations of ζ n and R π0 and their time derivatives, the same conclusion applies to their spectral indices. Thus, one can narrow the parameter region −1 ≤ c 2 b ≤ 0 where all the relevant curvature perturbations have an almost scale-free PS.
We start by computing amplitudes where the leading order SR expansion is sufficient. From (3.28), (B.14), (B.15) and (C.22) we get the following PS, valid for ;(3.34) ; (3.35) while for the cross-correlation we have As we will discuss in section 4, for the simplest reheating scenario the seed of primordial perturbations is given by the power spectrum of ζ n . Let us set where is the standard scalar PS in canonical single field inflation slow-roll prediction, whileP is a suitable dimensionless form factor depending on c 2 L , c 2 s1 , c 2 s2 and c 2 b that that can be read out from eq.(3.34). The value of the Hubble parameter during dS is denoted by H i . It is interesting to compare the above expressions with other existing models on the market. General single field models, in the effective field theory approach [26], when the sound speed is different from one, giveP = c −1 s ; while in adiabatic solid inflation modelP reduces to c −5 L (see (3.22)). ThusP can be interpreted as a sort of "effective sound speed" parameter in order to compare our predictions with different inflationary models. It should be stressed that the singular behavior when c 2 s1 or c 2 s2 is sent to zero or coincide, signals the simple fact that there is no way to change the number of propagating DoF in a controlled way. This, for instance, manifests when trying to re-obtain the adiabatic solid result from the supersolid one by imposing c s i → c L , resulting in a divergence proportional to 1/(c 2 s1 − c 2 s2 ), see (3.34). One can chooseP such that, in the stability region, the amplitude of the ζ n power spectrum is of order 10 −9 as required by observational constraints as shown in Figure 2 (this guess will be clearly demonstrated in section 4). We setP to a constant, extracting c 2 L as a function of c s2 , c s1 , c 2 b . When one of the two diagonal sound speeds tends to the longitudinal one, for instance c s1 → c L , then c L reduces to its maximal valueP − 1 5 . The maximal allowed area corresponds to a maximal longitudinal speed c L = 1 2 and P = 32. Thus, by taking 5 <P < 100, there is a sufficient large region in the parameters space spanned by c L and c s1 , c s2 to get a good agreement with data.
It is useful to study the behavior of the various power spectra when one of the sound speeds is much smaller than the other: c s2 c s1 withP fixed. From our findings (3.34), we get This gives, for the different values of c b that we are using, the approximate relations ; thus if we takeP = 32, the two speeds of sound are in the region: It is precisely the constraint onP that introduces a dramatic asymmetry, boosting P Rπ 0 . For c 2 b = −1, P Rπ 0 is naively enhanced by a factor 1/c s2 however taking into account 3.40, an extra enhancing factor 1/c 6 s2 is introduced; namely Similarly, for c 2 b = 0, an enhancement from 1/c s2 to 1/c 2 s2 is obtained So, the constraint imposed by the observed P ζn increases the sensitivity to small sound speeds of P Rπ 0 7 . Thus, at fixed P ζn , the P Rπ 0 becomes dominant being proportional to negative powers of c s2 , see 3 . In order to generate flat power spectra in a slow-roll regime requires that −1 ≤ c 2 b ≤ 0. Moreover, we were able to find an analytic solution for the modes at any time t only for c 2 b = −1 and c 2 b = 0; such values will be considered in the rest of the paper. It is convenient to trade the three independent parameters c 2 0,1,2 to c 2 L , c 2 s1 and c 2 s2 . While c 2 L = 4 3 c 2 2 − 1 can be interpreted as the longitudinal speed typical of Solid Inflation, the other two are the sound speeds corresponding to the two DoFs without any mixing (basically harmonic oscillators) described by (3.24). Finally, by matching the amplitude of the scalar power spectrum to the observed value of 10 −9 one can fix c 2 L as a function of the remaining two free parameters c 2 s1 and c 2 s2 . Let us note that our results, when comparable, do not agree with the one in [20]. The reason is the missing parameter c b (see footnote 3) and the treatment of the extra scalar degree of freedom in addition to the one present in solid inflation [9]. 7 For the cross-correlation Pζ n Rπ 0 we find the following expansion

Slowroll Corrections at Superhorizon Scales
In this section, we give the slowroll corrections of the primordial PS, focusing on the two exact solutions obtained for c 2 b = −1 and c 2 b = 0. Being the superhorizon behavior determined by the c 2 b value even at the leading order, the analysis of the large scales Π L and Π 0 fields in dS approximation is fundamental in order to get the parameter space where power spectra are scale free. By manipulating the system of second order coupled equations (3.27), in the large scale limit x = −k t 1 , we can get the following two independent fourth-order equations: The solutions can be expressed in the following form where the C coefficients are unspecified constant at this stage. The C L/0,1/2 refer to homogeneous solutions while C L/0 3/4 to particular solutions of the original system (3.27). Thus, we can understand the effect of c 2 b on the superhorizon evolution by obtaining particular solutions for Π L , sourced by Π 0 and vice versa. Table 2: Leading t and k behaviour for x = −k t → 0 of Π L,0 and ζ n , R π 0 for different values of c 2 b . Table 3: Inside the region of flat PS (−1 ≤ c 2 b ≤ 0) we give the late time structure of ζ n and R π 0 .
The boundary values for c 2 b = − 5 3 , −1, 0, 2 3 are shown in Table 3. The relation between ζ n and R π0 to the canonical fields asymptotically implies the following k and time dependence Taking into account the various transformations, the almost scale-free PS of ζ n and R π0 is obtained when the coefficient C L,1 and C 0,3 dominate. Thus, at the leading order, one realizes that ζ n is also always almost scale-free in the region − 5 3 ≤ c 2 b ≤ 2 3 , while R π0 selects a smaller region −1 ≤ c 2 b ≤ 0. We will focus on this last region. Recall that R and ζ are simple functions of (ζ n , R π0 ) and their time derivative, see (3.32,3.33).
Finally, let us outline the form of the relevant scalar fields at super-horizon scales, obtained by computing the next to leading slowroll corrections of the canonical normalized fields: the explicit form of the constants A a are not relevant here.
Let us explain what superscripts (ad) and (en) mean. The (ad) part of a field stands for its adiabatic part, in the sense that its adiabatic-tilt n (ad) s will not be affected by the presence of the c 2 b parameter strictly related to the presence of propagating superfluid density per lattice site (i.e. non-barotropic) perturbations. On the contrary, the superscript (en) stands for the entropic part of a field and the related tilts n (en) s will be c 2 b dependent. A crucial feature is that the behaviour of ζ n on super-horizon scales is determined by a single purely adiabatic power law given in terms of α (ad) and n (ad) s − 1. As we will show in the next section, in the case of an instantaneous reheating, it is precisely ζ n that determines the transition to the radiation era, setting the adiabatic part of the initial conditions; moreover, the non-adiabatic part will be determined by the difference ζ − ζ n . If c 2 b is far from the interval [−1, 0], the time dependence from α (en) will overwhelm the homogeneous R π0 solutions, leaving only its particular adiabatic solution. In practice, when we are far enough from the boundary values of c 2 b , also the other fields will be single-tilted; however in this region, the link between super and subhorizon amplitudes needs to be computed numerically.
On the contrary, at the leading order in slow-roll, when c 2 b = 0, −1, we get analytical solutions in the form of almost scale-free power laws for all relevant scalars (3.51) The above form was used in the previous section to compute the leading order amplitudes of the primordial PS for R π0 . In practice, when c 2 b = 0, −1 we cannot discriminate between the adiabatic and entropic R π0 parts. Such a degeneracy is removed by next to leading slowroll corrections. In the case of an almost instantaneous reheating, the slow-roll leading order computation of the Γ ΛCDM , primordial Non-Gaussianities and GWs back-reaction will be sufficient. For completeness we give the result of next to leading slowroll corrections n (ad) (3.52) The (ad) tilt is formally the same as the one found in [9], being obtained by solving the same superhorizon equations of motion. However, starting from a supersolid, the solid inflation limit simply does not exist: the extra degree of freedom cannot be smoothly switched off. In principle, two (en) tilts are present in the other fields and exists only on the edges of the region. Finally, let us sketch the main steps to get the above slow-roll corrections on superhorizon scales: 1. Trade the system of coupled equations for (ζ n , R π0 ) for an equivalent but simpler to analyze involving (R, δσ); 2. Find the canonical fields (R c , δσ c ) and find the leading superhorizon behaviour at the leading order in slow-roll as done for (Π L , Π 0 ). Define the (ad) part of R as its "homogeneous" (e.g. c b independent) component (coherent with the fact that ζ n is a purely "adiabatic" field) and the (en) part of R as the δσ-sourced solutions.
3. Compute the R and δσ slow-roll corrections on superhorizon scales.

Reheating
Once the seed of primordial perturbations is produced, it is important to study how the Universe reheats and gets to the radiation domination era. In single clock inflation, the hypothesis of the Weinberg theorem are satisfied [1] and the inflationary predictions are largely independent from reheating, however this not the case when more then one field is present, as for solid and supersolid inflation, where nor R neither ζ are conserved on super horizon scales and moreover R = ζ. As a consequence of the presence of ϕ 0 , the pressure perturbation is not proportional to δρ thus the Γ signals the presence of non-adiabatic perturbations. Dealing with more than one component like in ΛCDM, non adiabaticity can also be present when the relative energy density perturbations of two components are different: δ i = δ j . The total non-adiabaticity Γ tot contains both the intrinsic contribution for each component of the form (4.1) and the "relative" part Γ rel that takes into account that δ i is not simply caused by the "universal" temperature perturbation. In the case of ΛCDM with a barotropic equation of state for all the components only Γ rel is present and then Γ ΛCDM ≡ Γ rel ; at superhorizon scales one gets where ζ 0 is the adiabatic constant contribution. For a recent discussion see [7]. A pragmatic approach is to assume that reheating takes place instantaneously on a time-like hypersurface Q given in terms of a 4-scalar q as q =constant, or expanding at the linear order in perturbation theoryq + δq = constant .

(4.3)
A generic physical quantity ξ will be denoted by the subscript ξ − when evaluated at the end of inflation, and with ξ + when evaluated at the end of reheating phase. Thus, the change of ξ across Q will be simply written as and the transition will be dictated by the Israel junction conditions [27]. By generalizing the results in [28], in the general gauge (B.1), such conditions read at the linear level At the background level the junction conditions imply that both a and H are continuous on Q. The quantity ζ q represents the gauge invariant curvature perturbation of a constant q-hypersurface and thus it is continuous across Q. From the transformation properties (B.4), one can easily show that the junction conditions are gauge invariant. As a reasonable assumption, we will take Q to be the particle number density n. Intuitively, in the approximation of an instantaneous reheating, the rate for any channel for the decay of inflatons into a particle A becomes very large and the decay itself is democratic, in the sense if n A is the number density of the particle A and n is the total number density, then from the above relation and the particle number conservationn + 3 Hn = 0 we have that δn n ± = 0 . (4.9) In the flat gauge, where Φ = 0, such a condition is precisely (4.5) with q = n [ζ n ] ± = 0 , (4.10) When the field ϕ 0 is absent, namely M 0 = M 1 = 0 (solid inflation limit), one is back to standard case where reheating takes place at a constant energy density ρ hypersurface like in [9,21]. The continuity of ζ n can also be shown following the same lines of [21] by a generalization of the procedure given in [29]. By using the definition of ζ n and δσ we have that (4.11) By integrating by parts the relation which gives R and by using the time-time component of the Einstein equations, see [21], one gets [R] ± = 1 3 (4.12) where the effective intrinsic entropic perturbation Γ eff before/after reheating is defined as follows: Then (4.12), is equivalent to (4.14) demonstrating our intuition (4.10). By using (3.6) and (4.10) the second junction condition (4.6) reads Let us consider the most important case where, after the Universe reheats, a vanilla ΛCDM radiation dominated era is reached, for which at super-horizon scales From the above relation and by using (4.15), the jump of ζ across the reheating hypersurface is where, being ζ n continuous, ζ n+ = ζ n− and has been denoted simply by ζ n . Finally, one can calculate the total amount of non-adiabaticity Γ ΛCDM present at the beginning of the radiation era. Indeed, by comparison with (4.14) The jump of ζ is given by (4.17), thus finally, taking into account that the above relation refers superhorizon scale, from the results of Appendix B and C we arrive at There is still a point to address, take a generic field X that satisfies a second order evolution equation with two independent solutions: one X cg , growing or constant with scale factor a and a second one X d decreasing with a. Clearly, the physical relevant solution is X cg ; however, even if the junction conditions prescribe that [X] = 0, the constant/growing mode alone can be discontinuous. A classic example is given by the gauge invariant Bardeen potential Φ gi , which according to (4.7) is continuous in the transition at constant ρ with a sudden change of equation of state in ΛCDM; however, from the continuity of ζ constant mode, one gets Things are different in our non-adiabatic case and a clear understanding of the behaviour of Φ gi constant mode is crucial to predict the correct back reaction of tensor modes during radiation domination. Indeed, the validity of (4.20) crucially implies that the Φ gi gains a factor −1 entering radiation domination. For simplicity, in the rest of this section we will work in Newtonian gauge, where Φ gi coincides with Φ. For each classic scalar field it is convenient distinguish among a constant, decaying (absent during inflation) and entropic (particular solution proportional to the non-adiabatic source term proportional to Γ) modes. Once the decaying modes are under control, in principle, it is still possible a reshuffling of constant and entropic modes in the junction conditions. Focusing on the entropic source Γ ΛCDM relative to ΛCDM where dark energy is just a cosmological constant; neglecting baryons during radiation domination, we have two fluids: dark matter and photons as discussed in [7] and Γ ΛCDM assumes the form with s 0 (k) a scale dependent constant that is determined by using (4.19) at t + = −t − 8 . At superhorizon scales, the non-adiabatic contribution to ζ reads ζ| en = s 0 (k) a Ω m 3 a Ω m + 4 Ω r , (4.22) while the contribution to ζ n is (4.23) Thus, during inflation ζ n acts always as a source term for ζ, R and δσ when −1 ≤ c 2 b ≤ 0. The same exactly happens, during the radiation domination where any entropic contribution from ζ n is compensated by an opposite contribution from ζ or R leading to [ζ n ] | constant modes ≡ 0 . (4.24) Following [29], expressing ζ n in terms of the Bardeen potentials in the Newtonian gauge, we get eq. (4.20) is non longer valid. Imposing that we get (4.27) Considering the early stages of the radiation domination, where dark energy is negligible, we have Φ rad → − s 0 (k) 5 a 3 a 3 − 2 a 2 a e + 8 a a 2 e + 16 a 3 e , (4.28) where a eq = Ωr Ωm is the scale factor at the matter radiation equality and we normalized the today's scale factor as a 0 = 1. The same results could have been obtained by directly solving the Φ equation of motion or equivalently by expressing ζ in terms of Φ and Φ and enforcing that ζ n is not affected by non-adiabatic perturbations. Thus, eliminating Γ eff from (4.25), we can extract the constant Φ mode during the radiation phase which is similar to the standard result with ζ replaced by ζ n . The result (4.29) is not compatible with eq. (4.20) that would imply the transmission of the R π0 in the constant mode of Φ. Thus, Φ gets an enhancement of order −1 when the Universe transits into the radiation era without a further enhancement due to the presence of R π0 during inflation. Summarizing, in the case of an instantaneous reheating, ζ n determines initial conditions at superhorizon scales for the standard evolution for the ΛCDM scenario with small deviations from a perfectly adiabatic spectrum of primordial perturbations.

Primordial Non Gaussianity: a preview
Primordial Non-Gaussianity (NG) is very important to distinguish among different models of inflation. Single field inflation with its characteristic symmetry breaking patter gives a small amount NG both in the scalar and tensor sector, with the scalar part peaked in the local shape. A complete analysis of NG in supersolid inflation will be given in a companion paper [16], here we will outline some of the results needed to study the secondary production of GWs. Given the presence of two scalars and tensor fields, the full cubic action for a supersolid is quite complicated. Cubic terms can involve three scalars (SSS), one scalar and two tensors (TTS), two scalars and one tensor (TSS) and three tensors (TTT); each contribution to the cubic Lagrangian L (3) in Fourier representation has the following general structure L (3) ∼ ω a m D k k k ξ 1, k ξ 2, k ξ 3, k , where ω is a constant that sets the overall size of the vertex and m determines its time evolution in terms of the scale factor a; finally D is a function of the momenta and is determined by the structure of spatial derivatives acting on the fields entering the vertex denoted by ξ i,k which can be any combination of ζ n , ζ n , R π0 , R π0 , and h s ; h s is the spin two tensor field (indices are omitted). In general, one can show that The value of ω is determined by the relative size of the derivatives of the Lagrangian density U of the scalar sector with respect to the rotational invariant independent operators. In [9] it was assumed that there is partial cancellation among the derivatives of U such that, even in slow-roll, ω ∼ 1. Such extreme choice maximizes the deviation from single field inflation, pumping up local NG to f NL ∼ −1 which is in trouble with recent Plank constraints [30]. Here we take a more conservative approach, considering that each derivative of U is of order in slow-roll, leading to with α an order one quantity. As a result, we get that f NL ∼ O( 0 ) and, in addition, the cutoff of the effective field theory describing a supersolid is higher. Compared with NG in solid inflation, the presence of an additional scalar, related to the superfluid component, introduces non-adiabatic perturbations controlled by the parameter c 2 b which has an important effect on any 3-point function involving R π0 and, as we have seen, on the PS of R π0 itself as discussed in section 3.1.
In particular, when c 2 b → 0, we can show that the local f NL tends to be unacceptably big and strongly scale dependent, unless some rather unnatural tuning is made. As a result, when primordial NG is considered the best choice is to take c 2 b ≈ −1. As it will be shown in the next section, in this case supersolid inflation features a rather interesting boost of the secondary gravitation waves production during inflation thanks to the cubic mixed TSS that is promising for future experiments.

Gravitational Waves
Given the current experimental upper bound on the tensor to scalar ratio r ≤ 0.5, it is important to discriminate among different inflationary models by telling how much far or close to this upper limit a given model can be, considering that reasonably, in the next few years, we will be able to probe the range of r ∈ 10 −1 ÷ 10 −2 . Our analysis is similar to the one in [31], where secondary gravitational waves generated by a spectator scalar field was studied. However, in that specific case, taking into account the related secondary scalar PS, considerably reduces the ratio r [32]. On the contrary, in our supersolid model of inflation, the dominant cubic scalar vertex (SSS) is essentially unrelated to the dominant tensor-scalar-scalar (TSS) cubic one. That gives us room to effectively enhance r to get close to its experimental upper limit with only the secondary tensor production. That feature singles out supersolid from single field inflationary models where the dominant GW production is not very sensitive to NG and gravitational waves back-reaction is much smaller than the one generated during the radiation phase as it was observed originally in [33,34,35] and later extended in [36,37].
Spin two tensor perturbations in the metric tensor are defined by and the corresponding quadratic Lagrangian is where S ij is a transverse-traceless source term. The evolution equation for GWs is where we neglect the mass M 2 being proportional to , see (2.18). The leading contribution of the source terms is given by cubic interaction terms containing one spin two field h ij and it is composed by universal non linear interaction terms from the Einstein Hilbert Lagrangian and by interactions in the "matter" sector S ij = S . The leading structure of the EH interactions comes from derivatives of scalar perturbations and has the following structure The matter contribution S (Matter) ij , depending on the period, changes effectively during the universe evolution. In our case, during the inflationary period (where Matter = Inflaton medium), among all the possible TSS vertices, the dominant one is given by with α an order one constant. Given the presence of R π0 , the size of the source is very sensitive to the value of c s2 (typically ∝ c −3 s2 ). During radiation domination (Matter=Radiation Fluid), with photons represented as a perfect fluid, the source term becomes In our specific case, during inflation, we get that the Einstein Hilbert term is always suppressed S (Infl) S (EH) , while during the radiation phase nothing more than what is described in [36,37] happens; the only difference is that the Bardeen potential is proportional to ζ n instead of ζ, see (4.29). The tensor PS has two contributions: one (primary PS) P (1) h from the quantum fluctuations during the dS period and calculated with the homogeneous quadratic action of h ij , and an another classical contribution (secondary PS) P (2) h coming from the interactions of h ij with the other scalar fluctuations present and which can be calculated by finding the particular solution of (6.3) proportional to S (M atter) . The computation of the the primary tensor PS is standard and denoting with H i the Hubble parameter during the dS phase, we have Notice that ζ n fluctuations represent the primordial seed for scalar perturbations during the radiation phase. For the secondary PS, P h , the particular solution can be obtained by using the Green method, one has In the c 2 b = −1 case, the secondary PS is then where S k (t ) S p (t ) is a Gaussian 4-point correlator.
During the inflationary period the above 4-point correlator is proportional to R 4 π0 and in the limit of a small c s2 one gets the following estimate for the secondary scale-invariant PS P (2) h ≈ α 2 2 2 π 2 c s2 P 2 Rπ 0 = α 2 2 2 π 2 c s2 P 2 Rπ 0 = α 2 2 γ 2 2 π 2 c 13 s2P 2 P 2 ζn , (6.10) where we have defined γ such that The final expression for the total tensor PS is given by The presence of the coupling constant (6.5) which controls the TSS vertex gives rise to the question whether the cubic scalar interactions can give a sizable contributions to the scalar PS. The SSS dominant vertex has following structure in supersolid inflation where β is an order one dimensionless constant. As discussed in detail in [16], one can either choose β so that is small enough to safely neglect the secondary scalar production or take β ∼ 1. The second option makes the detailed secondary scalar production much more complicated; however, a reasonably estimate is given by ζn ∝ 2 β 2 2 π 2 c s2 P 2 Rπ 0 = 2 β 2 γ 2 2 π 2 c 13 s2P 2 P 2 ζn , (6.14) with the total scalar PS given by P ζn + P (2) ζn . The possibility to have a regime where the secondary tensor production is dominant while the secondary scalar contribution is negligible, namely Taking γ ∼ 1, P ζn = 10 −9 andP = 32 we get 0.15 β 2 2 1/13 c s2 0.2 α 2 1/13 . (6.17) Obviously, a small value β enlarge significantly the parameter region fur such possibility. Even the case where both the secondary scalar and tensor production are leading can be of interest. Clearly, a complete computation is more involved, a rough estimate gives which is of course very sensitive to the detail of the non-linear structure of the theory. Thus, there is a regime where the secondary production of gravitational waves during inflation can be comparable or greater than the primary one.

Conclusion
The pattern of symmetry breaking during inflation is very important. We have systematically explored the physical consequences of the breaking of full set of diffeomorphism of general relativity down to ISO(3). The breaking pattern is triggered by the background configuration of four scalar fields and, in order to allow dS spacetime as solution, we have considered an additional set of internal symmetries comprising SO(3) internal rotations and four shift symmetries. The four scalars ϕ A can by interpreted as the coordinates of a supersolid embedded in spacetime and the corresponding effective Lagrangian we have studied is the most general one consistent with the given symmetries at the leading order in a derivative expansion. As a comparison, in the effective description of single clock inflation [26] the residual symmetry comprises three dimensional diffeomorphism with one scalar and two tensor propagating modes, while in our supersolid inflation we have two scalars, two transverse vectors and two tensors. Interestingly, as a benefit of the supersolid interpretation, the scalar field fluctuations can be interpreted as phonons modes and non-adiabatic perturbations. Given the symmetry breaking pattern and the number of propagating modes, the difference with single clock inflation are important both a the linear and non-linear level. At the linear level, the symmetry breaking pattern gives rise to a peculiar kinetic mixing between the two scalars that makes the quantization and the computation of the linear power spectra non-trivial. A similar mixing is found in chromo-natural inflationary models [23,24,25], non-thermal production of gravitinos [22], multi-field inflation [23] and in effective theories of inflation [20]. Our analysis and results differ from the previous ones: we do not use perturbations theory to resolve the kinetic mixing but rely on Hamiltonian analysis and a set of canonical transformations to reduced the dynamical system to two uncoupled harmonic oscillators in the limit of large momentum k. As a consequence, crosscorrelations in the scalar power spectra are unavoidable. The presence of the scalar ϕ 0 associated to the superfluid component introduces the important parameter c 2 b for the superhorizon evolution of the scalars. The hypothesis of the Weinberg theorem are violated both because of the presence of the anisotropic stress which is not negligible in the k → 0 limit and also perturbations are in general nonadiabatic. Thus, neither the comoving curvature perturbation R nor the curvature perturbation ζ are conserved or equal on superhorizon scales, though their superhorizon evolution is only due to slow-roll corrections. In the range c 2 b ∈ [−1, 0] all thee relevant scalar power spectra are scale-free, modulo small slow-roll corrections, in agreement with experimental constraints. Because of the presence of two scalar propagating degrees of freedom there is no smooth limit that leads to solid inflation [9] and thus the predictions at the level of linear power spectra are rather different. The system of coupled second order differential equations for the linear evolution of the two independent scalar perturbations are complicated enough due to the non-trivial kinetic mixing to elude an analytical solution for a generic time t unless c 2 b = 0, −1. Luckily enough such values the boundary values for c 2 b where the relevant power spectra are almost scale-free. Among the various scalar perturbations, we select the power spectra of the curvature perturbation ζ n of the constant particle number n hypersurface and and curvature perturbation R π0 of the constant ϕ 0 hypersurface and the relative cross correlations, studying in detail their properties as a function of c 2 b and the speed of sounds of the two independent diagonal scalar modes. In the instantaneous reheating approximation, by extending the analysis in [28], we give a detailed analysis on how seed of primordial perturbations are transmitted to the standard hot radiation dominated era of ΛCDM, predicting, besides the standard adiabatic, a small isocurvature part which can be written as a linear combination of ζ n and R π0 evaluated at the end of inflation. Also the prediction for primordial non-Gaussianity are rather interesting; we leave a detailed account for a companion paper, focusing on the secondary production of gravitational waves during inflation. The structure of the tensor-scalar-scalar cubic vertex is such that is possible to enhance the secondary production, saturating the experimental bound, still keeping the scalar bispectrum within the limits set by Planck. In conclusion supersolid inflation is an interesting alternative to single clock inflation to explore different symmetry breaking patterns with a clear experimental signature.
Consider an infinitesimal coordinates transformation; in the scalar sector we have that By using this parameterization, R, ζ and σ can be easily written in terms of π L , π 0 and their derivative w.r.t. conformal time. Namely, with the suffix gi understood, we have The above expressions are valid at the linear order in the slow-roll parameters.

C Canonical Transformation
In the UV (large k) the Lagrangian (3.18) becomes . (C.7) The transformed Hamiltonian reads both K new and M new are diagonal and thus, for fixed k, (C.8) describes two uncoupled harmonic oscillators with frequencies k 2 c 2 s1 and k 2 c 2 s2 , where (C.10) Quantization of (C.8) is straightforward: the Bunch-Davies vacuum is vacuum state |0 of the Fock space corresponding to the following field operators k a k (2) e −i k cs2 t + A (C.16) By using (C.13-C.14), one can compute the free-field (Gaussian) average of any operator that can be expressed in terms of π L and π 0 . It is rather useful to simplify as much as possible the expression of power spectra to rewrite all the parameters of interest in terms of the two "diagonal" sound speeds c s1/2 , c b and c L defined in (3.22). One gets that , As a consequence of the fact that under the exchange of c s1 ↔ c s2 we have that C (1) L/0 , all power spectra will have the same property. Such a symmetry simply reflects the conventional choice c s2 < c s1 or c s1 < c s2 . Finally, note that in the two analytic cases c b = 0, −1, C L/0 sub-horizon coefficients are easily transmitted on superhorizon scales, and eq. (3.34) can be obtained considering that ζ n and R π0 classical solutions reduce to ζ n (j) → − C (j) (C. 22)