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 I SO(3). The effective description is based on four scalar fields that describe the excitations of a supersolid. There are two phonon-like propagating scalar degrees of freedom that mix non-trivially both at early and late times and, after exiting the horizon, give 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 and with a blue-tilted spectral index.


Introduction
Inflation is the most compelling way to solve the drawbacks of the hot big bang model and simultaneously generate the seed of the primordial perturbations to be used as initial conditions for the latter stages of Universe's evolution. The simplest class of models is single clock inflation, where time diffeomorphisms are 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 for which the Weinberg theorem does not hold. In this case, R and ζ are not conserved and R = ζ at superhorizon scales; thus, 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 -1 -

JHEP01(2021)185
a matter of fact, ϕ A can also be interpreted as the coordinates of a self-gravitating nondissipative medium [10][11][12][13][14][15] that in our case is a supersolid. A complete analysis of the linear evolution of scalar and tensor modes, together with the computation of the corresponding power spectra is given. In addition, we consider the secondary production of gravitational waves (GWs) during inflation, exploiting the cubic vertex of the theory involving a tensor and two scalars, saturating the experimental bound set by Planck without upsetting the scalar 3-point function. The secondary production can give a blue tilt to the spectral index, an important feature for the direct detection of the stochastic GWs background produced during inflation. A detailed analysis of non-Gaussianity 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 are 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.

Supersolids 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 ϕ 0 =φ(t) , ϕ i = x i (2.1) 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) internal symmetry 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 -2 -

JHEP01(2021)185
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 spacetime scalars shift symmetric operators with a single derivative of ϕ A (2.5) where plays the role of the medium four-velocity such that u µ ∂ µ ϕ a = 0 and By using the relation ς = χ −1 b y and the Cayley-Hamilton theorem, only 7 among those operators are independent. Thus, we arrive at the action that 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 it 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 JHEP01(2021)185 density n sf of the superfluid component can be expressed in terms of the Noether current J 0 µ as n sf = −u α J 0 α ; (2.12) while the density of lattice sites n is identified as the projection of the off shell conserved current 1 J α = b u α along the four-velocity u µ , namely [17] n = −u α J α = b . (2.13) 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 an isocurvature or entropic perturbation). 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 excitations are present. In general, the supersolid perturbation can be written around a flat space-time as 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 derivative of a function f with respect to conformal time is denoted by f . The parametersM 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 the related mass parameters) in the Lagrangian depending on a specific set of internal symmetries.
-U (b, y): the most general isentropic perfect fluid; the Lagrangian is invariant under V s Diff and ϕ 0 → ϕ 0 + f (ϕ a ).
As a remark, in the literature solids are typically associated with the presence of only three scalar fields {ϕ a , a = 1, 2, 3} and a Lagrangian of the form U (τ X , τ Y , τ Z ), see for instance [9]. However, introducing a fourth scalar ϕ 0 and enforcing the following field dependent shift symmetry the allowed operators are b, τ Y , τ Z and y and the resulting theory U (b, τ Y , τ Z , y) describes an adiabatic solid in the sense that the entropic perturbation δσ is a conserved quantity as discussed in [15,20]. The term supersolid is reserved to the case where, in addition to the phonons of the solid, also the entropic perturbation δσ propagates. A more detailed analysis of thermodynamical properties for general supersolids is planned for a future work. Finally, stability of (2) imposes the following conditions [19] As we will see, such conditions are necessary for the existence of the Bunch-Davies (BD) vacuum in an inflating 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 ( 3 The operator τX can be written as a function of (b, τY , τZ ), so only three operators are independent.

JHEP01(2021)185 3 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; the conservation of the background EMT is equivalent to 4 where Our benchmark values for c 2 b will be c 2 b = 0 which givesφ = aφ 0 , and c 2 b = −1 leading tō ϕ =φ 0 a −4 . Inflation takes place when We will be mainly interested in slow-roll (SR) inflation 5 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 space, one finds In [21] it was setφ = a that 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). 5 As discussed in [20] ultra SR is also possible; actually when M2 = 0, as for fluids, this is the only viable regime with small .

JHEP01(2021)185
The following action describes the linear dynamics at the leading order of a slow-roll expansion in : Up to boundary terms, one can always take D antisymmetric. The peculiarity of (3.9) is the mixing of 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 mixing is unavoidable unless the parameters c b and c 2 2 , c 2 1 are unnaturally tuned, and it is a key property of a superfluid component in the solid at 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 [21]. 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 (nondiagonal M). A similar kinetic mixing is also encountered in mechanical systems with gyroscopic forces like the Coriolis force or in the presence of magnetic fields; it is worth to stress that 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 two decoupled fourth-order differential equations.
The quantization of the system goes through the choice of the BD by studying 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 The quantum oscillator dynamics is recovered by a canonical transformation at Hamiltonian level (involving also the conjugate momenta) Introducing the Hamiltonian formalism allows us to decouple the two DoF with a canonical transformation and to select the BD vacuum. The main steps are summarized in figure 1.

Quantization and power spectra
In order to define the Power Spectrum (PS) of a general quantum scalar field ξ(x), in Fourier space we set where the (cl) subscript stands for classic solution, and the latin index i = 1, 2 indicates the two indepedent scalar modes whose annihilation and creation operators obey the standard canonical commutation relations (3.14) Thus, the ξ 2-point function reads 15) and the ξ scale-invariant PS is defined by (3.16) The first step to compute quantum correlators during inflation is to introduce the canonical field Π defined by The conditions (2.18) 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 (3.20) , (3.21) and we have defined At the leading order in slow-roll, the equations of motion are the following 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 timeindependent. 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 Π L (UV) 2 Thus, the unique Fock space vacuum is the BD vacuum for the system. Details can be found in appendix C. The existence of the BD vacuum requires the frequencies squared ω 2 1 = k 2 c 2 s1 , ω 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 conditions (2.18) are sufficient to have c 2 si > 0 and when expressed in terms of (2.19) gives 6 (3.25) 6 We assume the null energy condition 1 + w > 0.

JHEP01(2021)185
We have checked that there is a large region of parameters 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 system of second order linear 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 for which an analytic solution can be found. Neglecting SR corrections, the coupled system of second order equations can be written as two indepedent fourth order equations for Π 0,L . Remarkably, the case c 2 b = 0 and c 2 b = −1 gives the following identical equations valid if Analytic solutions are possible due to the absence of the terms Π (3) i 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 1, the solution that represents flat space modes matches the ones given in (C.13) and (C.14); such a choice is equivalent to select the BD vacuum. Thus, Π L and Π 0 are quantum free (Gaussian) fields given by (3.28) the expression for C 0 ; j = 1, 2 can be found in appendix C and a (j) 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

JHEP01(2021)185
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 ( keep in mind that n = n ), (2.13)) and curvature of the ϕ 0 =const. hypersurface; namely The comoving curvature R π 0 is related to the superfluid component (2.15). The expression of the various curvature perturbations in terms π L and π 0 in a generic gauge can be found in appendix B. In the spatially-flat gauge (3.1) we have that The uniform curvature perturbation ζ can be obtained from eq. (B.15) at leading order in SR 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, in the region −1 ≤ c 2 b ≤ 0, all the relevant curvature perturbations have an almost scale-free PS. From (3.28), (B.14), (B.15) and (C.22), at leading order SR expansion, we get for ; (3.34) where we have introduced the scalar PS P SF in canonical single field inflation where H i denotes the value of the Hubble parameter during dS. 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 whereP is a dimensionless form factor depending on c 2 L , c 2 s1 , c 2 s2 and c 2 b that can be read out from eq. (3.34). 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 of the PSs 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 trying to re-obtain the PS for an adiabatic solid result from the supersolid one by imposing c s i → c L , leading to divergence proportional to 1/(c 2 s1 − c 2 s2 ) as one can see from (3.34). In the stability region, one can chooseP such that the amplitude of the ζ n power spectrum is of order 10 −9 as required by observational constraints as shown in figure 2. 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 andP = 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), 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 figure 3. Let us briefly recap the relevant used parameters. The quadratic Lagrangian contains 5 mass parameters (2), (2.19) or equivalently 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 . 7 For the cross-correlation Pζ n Rπ 0 we find the following expansion  While c 2 L = −1+ 4 3 c 2 2 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 [21]. 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].
Taking c s1 and c s2 as independent parameters, the behavior of c 1 and c 0 is important in the study of the amount of isocurvature perturbations and the secondary gravitational waves production. The complete expression is given in appendix C, we quote here the leading contribution for small c s2 (3.46) Note that we have used that P ζn = const. ∼ 10 −9 as a fixed constant value. Thus, both c 0 and c 1 are suppressed for small c s 2 as shown in figures 4 and 5.

Slow-roll corrections at superhorizon scales
In this section, we give the slow-roll 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 powerspectra are scale-free. By manipulating the system of second order coupled equations (3.27), in the large scale limit x 1 , we can get the following two independent -15 -JHEP01(2021)185  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.
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 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 slow-roll corrections of the canonical normalized fields: (3.51) 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. This parameter is 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 tilt n (en) s will be c 2 b dependent. A crucial feature is that the behavior of ζ n on superhorizon 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 singletilted; 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.52) 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 slow-roll 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 slow-roll corrections

JHEP01(2021)185
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 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. 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.

Degeneracy breaking:
• for c 2 b = −1, δσ is a dominant decoupled (en) source on superhorizon scales, which means that δσ = δσ (ad) + δσ (en) ≡ δσ (en) ⇒ δσ (ad) (ζ n , R π 0 ) = 0 ; (3.54) • When c 2 b = 0, R is a dominant decoupled (ad) source on superhorizon scales, which means that R = R (ad) + R (en) ≡ R (ad) ⇒ R (en) (ζ n , R π 0 ) = 0 . (3.55) Following the above steps, one arrives at eq. (3.51) and, in addition, the degeneracy is resolved by (3.56) A degeneracy in the amplitude persists for c 2 b = −1 . In that case, by using eq. (3.56), we get As will be shown in the next section, ζ n provides the seed for adiabatic perturbations at the beginning of radiation domination; as a result CMB data [27] imply that its spectral index n (ad) s has to be red tilted. Thus, when c 2 b = 0, R π 0 will be red-tilted too. However, which can be blue-tilted. The consequences for the secondary production of gravitational waves is rather interesting and studied in section 6.

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 of reheating, however this is not the case when more then one field is present, as for solid and supersolid inflation, where neither R nor ζ 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 timelike hypersurface given in terms of a 4-scalar q as q =constant, or expanding at the linear order in perturbation theoryq + δq = constant . A generic physical quantity F will be denoted by the subscript F − when evaluated at the end of inflation, and with F + when evaluated at the end of reheating. Thus, the change of F across q will be simply written as

JHEP01(2021)185
and the transition will be dictated by the Israel junction conditions [28]. By generalizing the results in [29], 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 qhypersurface, 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 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 the standard case where reheating takes place at a constant energy density ρ hypersurface like in [9,20]. The continuity of ζ n can also be shown following the same lines of [20] by a generalization of the procedure given in [30]. 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 [20], one gets [R] ± = 1 3 (4.12) where the effective intrinsic entropic perturbation Γ eff before/after reheating is defined as follows: -20 -

JHEP01(2021)185
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 superhorizon 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 to superhorizon scales, from the results of appendix B and C we arrive at Note that the contribution to the transmitted Γ stays small for small c s2 . Indeed, from eqs. (3.46), we get that 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 physically relevant solution is -21 -

JHEP01(2021)185
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. A clear understanding of the behavior of Φ gi constant mode is crucial to predict the correct back reaction of tensor modes during radiation domination. Indeed, the validity of (4.21) 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 to distinguish among 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, a reshuffling of constant and entropic modes in the junction conditions is still possible. 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.23) while the contribution to ζ n is (4.24) 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 to ζ n is compensated by an opposite contribution from ζ or R leading to (4.25) Following [30], expressing ζ n in terms of the Bardeen potentials in the Newtonian gauge, we get we get (4.28) Considering the early stages of radiation domination, where dark energy is negligible, we have 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.26), we can extract the constant Φ mode during the radiation phase which is similar to the standard result with ζ replaced by ζ n . The result (4.30) is not compatible with eq. (4.21) 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 an essential tool to distinguish among different models of inflation. Single field inflation with its characteristic symmetry breaking pattern gives a small amount NG 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 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 dimensionless function of the momenta and is determined by the structure of spatial derivatives acting on the fields entering the

JHEP01(2021)185
vertex denoted by ξ i,k which can be any combination of ζ n , H −1 ζ n , R π 0 , H −1 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 the presence of a 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 [31]. Here we take a more conservative approach, considering that each derivative of U is of order in slow-roll expansion, 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 introduces non-adiabatic perturbations controlled by the parameter c 2 b . This parameter 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 exciting 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 close to the limit the prediction for r can be. Indeed, in the next few years, we will be able to probe the region ∈ 10 −1 ÷ 10 −2 . Our analysis is similar to the one in [32], 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 [33]. 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 [34][35][36] and later extended in [37,38].
Spin two tensor perturbations are defined by

JHEP01(2021)185
where h ij is the transverse and traceless part of the metric tensor. During Inflation, the corresponding quadratic/cubic Lagrangian can be written as where S ij is a transverse-traceless quadratic-source term. The evolution equation for GWs is where we neglect the mass M 2 being proportional to , see (2.19). The leading contribution to S ij comes from the cubic interaction terms containing one spin two field h ij and two scalars. There is a "universal" contribution from cubic terms in the Einstein-Hilbert Lagrangian and a graviton scalar interactions in the "matter" sector; namely . The leading structure of the EH interactions comes from derivatives of scalar perturbations and has the following structure The matter contribution S (Matter) ij changes effectively during the universe evolution. 9 In our case, during the inflationary period (where Matter = Inflaton), among all the possible TSS vertices, the dominant one is given by the following cubic lagrangian (see the structure in (5.1)) where c 2 2 ≡ c 2 2 τ +c 2 2 w ; c 2 2 τ = − a 2 9 H 2 U τy + U τz ; c 2 2 w = − a 2 9 H 2 U wy + U wz ; (6.8) with c 2 2 w the part of M 2 10 proportional to the derivatives of U with respect to the operators The analysis of the role of the operators w i is interesting. At the zero and first order in the perturbation theory they are degenerate with the operators {τ i }, so are sensitive only to the solid structure of the medium. It is only at second order that the {w i } start to discriminate a solid from a supersolid. From the structure of the {w i }, we see that the ϕ 0 scalar field, related to the superfluid part, is intrinsically coupled to the ϕ a fields, 9 During matter/radiation domination (Matter=Matter/Radiation Fluid), with DM/photons represented as a perfect fluid, the source term becomes . 10 See appendix A. and c 2 s 2 < c 2 L < c 2 s1 . Given the presence of R π 0 , the size of the source is very sensitive to the value of c s2 (typically ∝ c −3 s2 ). 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 [37,38] happens; the only difference is that the Bardeen potential is proportional to ζ n instead of ζ, see (4.30). The tensor PS has two contributions: one (primary PS) P Remember that ζ n fluctuations represent the primordial seed for scalar perturbations during the radiation phase. The particular solution of (6.3) can be obtained by using the Green function (g k (t, t )) method (6.10) The above solution can be used to extract the PS P (2) h for the secondary production of GWs as where S k (t ) S p (t ) represents a Gaussian 4-point scalar correlator. During the inflationary period, the above correlator is proportional to R 4 π 0 ; in the limit of a small c s2 one gets the following estimate for the secondary scale-invariant PS where we have defined γ such that

JHEP01(2021)185
The final expression for the total tensor PS is given by 14) The presence of the coupling constant (6.6) which controls the TSS vertex gives rise to the question of whether the cubic scalar interactions can give sizeable contributions to the scalar PS. The SSS dominant interaction Lagrangian for the scalar one loop corrections to P ζn has the following structure (5.1) and when c 2 b ≈ −1 As usual, U by and σ will be taken to be of order to get β i order one, as discussed in the previous section. Furthermore, while β 1 vanishes in the absence of w-operators, β 2 is generically different from zero. For generic β i ∼ 1, the computation of the non-linear correction to the scalar PS is complicated. A reasonable estimate is given by 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 is valid only if is very small, and then is much more tuned. The presence of the parameter c 2 2 w makes the gravitational waves secondary production dominant for a suitable region of the parameters space, even if the R π 0 is not efficiently transmitted in the scalar sector -27 -

JHEP01(2021)185
after inflation. Even the case where the secondary scalar and tensor production are both relevant is interesting. A rough estimate gives which is very sensitive to the detail of the non-linear structure of the theory. Let us mention that, even if a cubic coupling between tensors and transverse vectors of the schematic form h ij (π i π j + ∂π i ∂π j ) exists, we do not expect an enhancement similar to the one found due to the scalars. The transverse vector sector is very similar to solid inflation, and in the limit of small c s 2 with π i defined by the linear theory. Thus, the secondary production of GWs from the vector sector is much smaller than the one from the scalar sector which is of order P Rπ 0 2 . All the above expressions are given at the leading order in a slow-roll expansion and the PS are scale-free, modulo small slow-roll corrections. Finally, let us estimate the tilt of secondary GWs production. The next to leading corrections to the primordial tensor PS can be obtained by simply substituting P Rπ 0 in eq. (6.12) with the complete expression k 3 |R π 0 | 2 /(2 π 2 ). In the c 2 b = −1 case, we have three contributions with three different tilts (see eq.
which is blue tilte when The presence of a blue n T parameter will be an interesting tool to test inflationary models in future high sensitivity experiments of Gravitational waves detection [39]. Indeed, for modes that re-enter the horizon during radiation domination, the GWs energy density spectrum [40,41] goes as Ω GW ∝ k n T −1 , (6.28) 11 The standard CMB-like pivot scale is k * = 0.  Table 4. Relations among power spectra (leading order) and their fate in a instantaneous reheating, in the limit c s2 1. ad = adiabatic (P if n T − 1 > 0, it could grow up to the frequencies f = c k/2π in the milli-Hertz band, where we expect the maximum of LISA sensitivity [42,43].

Conclusion
An extreme synthesis of the production of the seeds for cosmological perturbations in supersolid inflation is given in table 4 where the magnitude and the fate of the leading order power spectra of scalar and tensorial perturbations are shown. Adiabatic perturbations are related to the solid part of the medium (ζ n ∝ π L ∝ ϕ i ). The presence of large entropic perturbations, related to its superfluid component (R π 0 ∝ π 0 ∝ ϕ 0 ), potentially can enhance the PS of the other fields by next to leading corrections. The secondary production is generically suppressed for inflaton-like fields, 12 tensor perturbations play the role of spectator fields (with a PS proportional to H 2 i /(M 2 pl )) and the interaction with entropic scalar perturbations enhances considerably their secondary PS. Indeed, the following scenario is possible. 13 We have systematically explored the physical consequences of the breaking of the 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 a solution, we have considered an additional set of internal symmetries comprising SO(3) internal rotations and four shift symmetries. The four scalars ϕ A can be interpreted as the coordinates of a supersolid embedded in 12 Inflaton perturbations in a quasi-deSitter background has its modes proportional to H 2 i /( M 2 pl ). In the supersolid case we have two scalars (ζn and Rπ 0 ) and one transverse vector field π i T . Transverse vectors decay subhorizon, we also expect a suppressed contribution to the secondary PS and will not be discussed in this paper. 13 Actually, in section 6 we get an extra 1/cs2 enhancing factor from phase space integration.
ad and P (2) h refer to the next to leading contribution to the adiabatic scalar and tensor power spectra. The above consideration are valid in the case of a small cs2. From the above results we can extract some general conclusions related to the pattern of symmetry breaking during inflation.

JHEP01(2021)185
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 significant both at the linear and non-linear levels. 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 (but different) 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 [21]. 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, cross-correlations in the scalar power spectra are unavoidable. The presence of the scalar ϕ 0 associated with the superfluid component introduces the important parameter c 2 b for the superhorizon evolution of the scalars. The hypothesis of the Weinberg theorem are explicitly violated. Indeed, we get both the presence of the anisotropic stress which is not negligible in the k → 0 limit, and perturbations are non-adiabatic in general. Thus, neither the comoving curvature perturbation R nor the curvature perturbation ζ are conserved, moreover they differ on superhorizon scales, though their superhorizon evolution is only due to small slow-roll corrections. In the range c 2 b ∈ [−1, 0], all the 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, these boundary values for c 2 b are such that 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 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 [29], we analyze how the seed of primordial perturbations are transmitted to the standard hot radiation dominated era of ΛCDM. Besides the standard adiabatic component, a small isocurvature part 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 is rather interesting; we leave a detailed account for a companion paper, focusing on the secondary production of gravitational waves during -30 -

JHEP01(2021)185
inflation. The structure of the tensor-scalar-scalar cubic vertex is such that it is possible to enhance the secondary production, saturating the experimental bound, still keeping the scalar bispectrum within the limits set by Planck. Finally, the spectral index of GWs PS can be blue-tilted, enhancing the chance of a direct detection of the primordial stochastic background.
In conclusion, supersolid inflation is an interesting alternative to single clock inflation to explore different symmetry breaking patterns with a clear experimental signature.
The first step is to find the Hamiltonian density corresponding to (3.18), which reads where P is the conjugate momentum of Π P = Π + D UV Π . (C.4) The decoupled system can be obtained by applying a canonical transformation of the form: .
(C.7) The transformed Hamiltonian reads (C.9) 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 Quantization of (C.8) is straightforward: the Bunch-Davies vacuum is the state |0 of the Fock space corresponding to the following field operators where C (j) (C.16) By using (C.13)-(C.14), one can compute the free-field (Gaussian) average of any operator expressed in terms of π L and π 0 . In order to simplify as much as possible the expression -34 -

JHEP01(2021)185
of power spectra, it is rather useful 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). We get Note that under the exchange of c s1 ↔ c s2 we have that C (1) L/0 . As a consequence, all the 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

(C.22)
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.