Higgsed Gauge-flation

We study a variant of Gauge-flation where the gauge symmetry is spontaneously broken by a Higgs sector. We work in the Stueckelberg limit and demonstrate that the dynamics remain (catastrophically) unstable for cases where the gauge field masses satisfy γ < 2, where γ = g2ψ2/H2, g is the gauge coupling, ψ is the gauge field vacuum expectation value, and H is the Hubble rate. We compute the spectrum of density fluctuations and gravitational waves, and show that the model can produce observationally viable spectra. The background gauge field texture violates parity, resulting in a chiral gravitational wave spectrum. This arises due to an exponential enhancement of one polarization of the spin-2 fluctuation of the gauge field. Higgsed Gauge-flation can produce observable gravitational waves at inflationary energy scales well below the GUT scale.

With the increasingly exquisite measurements of the spectrum of temperature and polarization fluctuations in the Cosmic Microwave Background, the basic inflationary paradigm is in good shape [11]. The measured fluctuations are adiabatic, Gaussian, and there is evidence at the 5-σ level of a red tilt from the CMB alone. While there is currently no evidence for gravitational waves, upcoming experiments such as CMB Stage 4 [12] will probe tensor-to-scalar ratios as low as r ∼ 10 −3 .
In this work we study a massive or Higgsed variation of the model of inflation called Gauge-flation, first proposed in refs. [13,14]. 1 The remarkable aspect of Gauge-flation is that it does not contain scalar fields. Instead the theory utilizes non-Abelian gauge fields in a classical configuration to generate an epoch of accelerated expansion. The Gaugeflation model can be obtained from a related model, Chromo-Natural inflation [16][17][18], by integrating out an axion about the minimum of its potential [17,19,20]. Unfortunately, Gauge-flation and Chromo-Natural inflation are ruled out at the level of the fluctuations [21][22][23][24]. In regions of parameter space that yield acceptable scalar density fluctuations, the tensor-to-scalar ratio is too large; conversely, in the regions where the tensor-to-scalar ratio is acceptable, the scalar spectrum is too red-tilted. In this work, we augment the Gaugeflation model by introducing a Higgs sector which spontaneously breaks the gauge symmetry [25,26]. Recently we demonstrated that breaking the gauge symmetry in Chromo-Natural inflation allows that model to generate spectra that are consistent with current data [26], and in this paper we demonstrate that Gauge-flation too can generate acceptable spectra in a broken phase.
While there is certainly no shortage of inflationary models on the market [27], and many that fit the data well [28], most rely on (one or multiple) slowly rolling scalar fields to generate an extended period of nearly exponential expansion. In these models of inflation, the amplitude of the gravitational wave spectrum is set only by the energy scale during inflation. Obtaining a large amplitude gravitational wave spectrum generically requires that inflation occurs at an energy scale near the energy associated with grand unification, and a large tensor-to-scalar ratio requires that the inflaton roll a distance in field space that is comparable to the Planck scale [29]. As we demonstrate, the remarkable feature of Higgsed Gauge-flation is the generation of observable gravitational waves at much lower energy scales -in this model, gravitational waves mix with exponentially enhanced gauge field fluctuations, resulting in their subsequent amplification. This phenomena is also observed in Higgsed Chromo-Natural inflation [26] and in models of inflation that have an accompanying spectator Chromo-Natural inflation-like sector [30][31][32] Classical non-Abelian gauge fields lead to striking phenomenology in cosmological settings [33,34], most notably chiral gravitational waves [23,30,31,[35][36][37][38][39][40][41][42][43] and the facilitation of gravitational leptogenesis [44]. Classical non-Abelian gauge fields have also recently been employed in generalized multi-Proca theories [45] to build stable inflationary models that do not require gauge invariance [46], and to generate inflation models with Horndeski couplings [47,48].
Throughout this work, we use natural units where the speed of light and the reduced Planck constant are set to unity, c = = 1.

Gauge-flation: Inflation from non-Abelian gauge fields
We consider the theory of Gauge-flation [13,14], which is described by the action and consider a general SU(2) gauge field, A µ , adopting the conventions of Peskin and Schroeder [49] for its action. In particular, the field-strength tensor and covariant derivative are defined as 2 where g is the gauge field coupling, not to be confused with the determinant of the spacetime metric. We normalize the trace over the SU(2) matrices, which we denote J a , so that where abc are the structure functions. The dual field strength is definedF µν = µναβ F αβ /2, and our convention for the antisymmetric tensor is 0123 = 1/ √ −g, while our spacetime metric signature is (−, +, +, +). Here and throughout, Greek letters denote spacetime indices, Roman letters from the start of the alphabet denote gauge indices and Roman letters from the middle of the alphabet denote spatial indices. Appendix A outlines our remaining conventions and notations.
In addition to the field content of Gauge-flation, we consider the addition of a symmetry breaking sector proposed in ref. [26], 3 which we write in Stueckelberg form [50,51] The fields ξ a are the Goldstone modes corresponding to fluctuations of the Higgs along its vacuum manifold.

Background solutions
The background evolution in Gauge-flation is found by considering the gauge fields in the classical flavor-locked configuration where J a is a generator of SU(2) satisfying the commutation relations and f abc are the structure functions of SU (2). Note that for SU (2), f ijk = ijk , where ijk is the completely antisymmetric tensor in three-dimensions. On the background field configuration in eq. (2.5), the field strength tensor components are, (2.7) Here and throughout a prime, , denotes a derivative with respect to conformal time. This field configuration results in a stress tensor that is consistent with the symmetries of Friedmann-Robertson-Walker spacetime. For these degrees of freedom, the mini-superspace action takes the form (see also ref. [25]) where an overdot here and throughout represents a derivative with respect to cosmic time, and N = a on the background. Sheikh-Jabbari and Maleknejad [13,14] demonstrated the existence of inflationary solutions for this system in the absence of the symmetry breaking terms (Z 0 = 0). They pointed out that, while the terms in the action arising from the Yang-Mills field have the equation of state of radiation, p = ρ/3, where the term proportional to κ has the equation of state of a cosmological constant. That is (2. 10) This implies that if ρ κ ρ YM , then the background spacetime undergoes a phase of accelerated expansion.
The addition of the symmetry breaking sector generates additional contributions to both the energy density and the pressure, note the equation of state is w = −1/3, and thus the presence of the symmetry breaking sector does not affect the conditions for accelerated expansion, which remains ρ YM ρ κ [25]. However, since ρ Z 0 + p Z 0 = 0, successful slow roll inflation requires ρ Z 0 ∼ ρ YM ρ κ in order to ensure H 1 (see eq. (2.15) below). The equation of motion for the gauge field vacuum expectation value (vev) that follows from the action at eq. (2.8) is Hφ a + g 2 Z 2 0 φ a = 0. (2.12) For the remainder of this work we instead use the variable ψ = φ/a, in terms of which, eq. (2.12) isψ The equations of motion for the metric are the Friedmann constraint and which can be combined to read We introduce the standard Hubble slow roll parameters, (2.20) Alternatively, using eq. (2.16), can be expressed as Equations (2.21) and (2.20) can be used to express κ and ψ in terms of , γ, δ and M 2 The parameter γ is identical to the parameter m 2 ψ defined in ref. [26]. In this work we use γ to be consistent with the nomenclature frequently used for gauge-flation, as for example in ref. [13,14,21]. We use M for the contribution to the mass due to the Higgs VEV, rather than ω as used in ref. [25] Differentiating eq. (2.21) with respect to cosmic time and using eq. (2.20) gives the exact expression Performing the same differentiation on eq. (2.20), we arrive at the equivalent exact relation Using the exact relations obtained above, the slow-roll parameters can be shown to satisfy the relations [25] to lowest non-trivial order. Note that η = O( ) and δ = O( 2 ), which implies that the relative change in ψ during inflation is much smaller compared to the corresponding change in H. Therefore, to a very good approximation, ψ ≈ constant throughout inflation. The Hubble parameter can be re-written using the definition of γ and the slow-roll approximation of as Finally, the total number of e-folds of inflation can be conveniently expressed only in terms of initial values as We end this section with a comment on the parameters required to fully characterize the background evolution of the system. As in Gauge-flation, the parameter κ can be eliminated by rescaling time t → t √ κ and the gauge coupling g → g/ √ κ [21]. This rescales the value of the Hubble parameter as H → H/ √ κ, and thus the value of κ is determined by fixing the overall amplitude of the scalar spectrum to A s 2 × 10 −9 . The remaining parameters of the theory are g and Z 0 .
In what follows, we use e-folding number, N = − ln(a/a 0 ), as our time parameter. For convenience, we choose a 0 = 1 at N = 60 e-folds before the end of inflation where necessary. In order to specify a background trajectory, (ψ(N ),ψ(N )), we need to specify the set {H in , ψ in ,ψ in , γ in , M in , g, Z 0 } for some initial N in . However, note that these seven quantities are not all independent. We can use the definition of δ to expressψ in aṡ Figure 1. The total number of e-folds of inflation using the full numerical evolution of the system (red) and using eq. Together with the Friedmann constraint eq. (2.14), eq. (2.28) and the definitions at eq. (2.19) provide four relations among the seven variables. This reduces the required parameters to three. If we specify γ in , ψ in , and M in , eq. (2.27) determines the length of the inflationary phase. Figure 1 shows the resulting number of e-folds for various parameter combinations. Increasing all three parameters (ψ in , γ in and M in ) leads to a decrease in the total number of e-folds of inflation. However, there remains a large region of parameter space where sufficient inflation is easily achieved. In the evolution of the perturbations we present below, we choose to specify γ in and M in at N in = 60 e-folds before inflation ends. In this case, eq. (2.27) specifies ψ in . Since eq. (2.27) is a very good approximation for all γ, M , and ψ, with the further (excellent) approximation that ψ ≈ ψ in , eq. (2.27) and eq. (2.14) can be solved for the subsequent values of M and γ as a function of N , the e-folding number measured with respect to the end of inflation.

Linear perturbations
In order to find the spectra of density and gravitational wave fluctuations in Higgsed Gauge-flation, we need to understand how the field and metric fluctuations evolve. In this section we derive the action to quadratic order in small fluctuations about the solutions described above in section 2.1. We begin by deriving the action for a the fluctuations of a general SU(2) gauge field about the background field trajectory before we specialize to a two-dimensional representation and introduce a scalar-vector-tensor decomposition of the fluctuations in section 3.1. Sections 3.2, 3.3, and 3.4, study the scalar, vector, and tensor fluctuations, respectively.
To proceed, we write the metric in ADM form [52], (3.1) where N is the lapse, N i is the shift vector, andh ij is the metric on the spatial hypersurface. At zeroth order in fluctuations, the FRW metric in conformal time corresponds to N = a and N i = 0 in our conventions. The metric on the hypersurface,h, can be decomposed into scalars, vectors and tensors by writing where γ ii = ∂ i γ ij = 0, and ∂ i C i = 0. The coordinate invariance of general relativity allows to impose four conditions on the fields in eq. (3.2). For this work, we choose spatially flat gauge, where the time threading and spatial coordinates are chosen so that A = B = 0. The remaining spatial reparametrizations can then be used to set C i = 0, which completely fixes the coordinates. We further write 5 so that det[h ij ] = a 8 to all orders in perturbation theory. Inserting the ADM metric at eq. (3.1) into the action in eq. (2.1) we find, In this expression, E ij is related to the extrinsic curvature of the spatial slices and ∇ i is the covariant derivative constructed fromh. Note that spatial indices are raised and lowered usingh ij and its reciprocalh ij . We proceed by expanding the lapse and shift about a Friedmann-Robertson-Walker spacetime where α (1) and N i (1) , and α (2) and N i (2) , are first and second order in fluctuations, respectively. As is well known, in order to obtain the quadratic action we require only the constraints at linear order, and thus we drop the subscripts in what follows.
In spatially flat gauge, neglecting gravitational waves for a moment, the curvature of the spatial slices vanishes, 3 R = 0, and the connection for ∇ i (the covariant derivative compatible with the metrich ij on the hypersurface) vanishes and thus ∇ i → ∂ i . The Einstein-Hilbert action to quadratic order in scalar and vector fluctuations is given by We denote the fluctuations in the gauge field by in terms of which the gauge field and Stueckelberg action to quadratic order in field fluctuations, and scalar and vector metric fluctuations is given by (3.8) The terms δ 2 L YM and δ 2 L κ are given by and we have defined and The Goldstone modes contribute at quadratic order in fluctuations via The addition of a Higgs sector thus yields an additional mass term for the gauge field fluctuations. Note, however, that retaining gauge-invariance requires us to also add the Goldstone modes ξ a . Finally, the quadratic Lagrangian density for the transverse-traceless components of the metric, and their interactions with the gauge field fluctuations is given by where γ 2 = γ ij γ ij . In order to proceed we need to choose a specific representation for the gauge field. We focus on a two-dimensional representation in what follows for simplicity.

Two dimensional representation and scalar-vector-tensor decomposition
Specializing to the case of a N = 2 dimensional representation of SU(2), the representation matrices are the Pauli matrices, J a = σ a /2, and we can decompose the gauge field fluctuations into scalar, vector, and tensor fluctuations. In order to make contact with the existing literature, we decompose the gauge field, Goldstone, and metric fluctuations as where Y , θ, α, δψ, U , ξ, and M are scalars; Finally, t ia is a transverse and traceless tensor t ii = ∂ i t ia = ∂ a t ia = 0. We fix the gauge for the gauge field fluctuations by setting U = U j = 0, (3.19) which is equivalent to choosing Ψ a i to be symmetric under exchange of i ↔ a. At quadratic order, the Lagrangian separates into separate scalar, vector, and tensor pieces as usual, and in what follows we consider each type of fluctuation separately.

Scalar fluctuations
After gauge fixing, there are five scalar fluctuation degrees of freedom in this theory, δψ, M, and Y , which arise from the gauge sector, ξ arising from the Higgs sector, and α and θ that arise from the metric perturbations. The quadratic action for these degrees of freedom reads where we have integrated by parts, discarded a boundary term, and made use of the background equations of motion. Note that in the limit Z 0 → 0, the above action does not quite agree with the corresponding expression in ref. [21]. The difference arises due to a slightly different choice of parametrization of the gravitational constraints. In this work, we have chosen to split the metric using ADM variables. Note that the fields Y , α, and θ appear in the action without time derivatives. As described in detail in ref. [21], these fields are algebraic constraints, and can be integrated out by solving their linear equations of motion and substituting the solutions back into the action. While this is a straightforward procedure, the result is extremely messy and we do not reproduce it here.
Denoting by X = (δψ, M, ξ), we redefine the fields using the transformation After integration by parts and discarding boundary terms, this field redefinition puts the action in the form In principle it is possible to choose a redefinition, U, that sets the kinetic matrix, T, to the identity matrix, however, this requires a much more complicated transformation which makes the algebra much more involved. As discussed in detail in ref. [21], this is not necessary to evolve the fluctuations. All that is required for our purposes is that T approaches the identity in the limit k aH in order to impose the initial conditions. While the matrices T, K, and Ω 2 are obtained in a fairly straightforward manner as we have described above, they are extremely long, and their full form is not particularly illuminating. In appendix B, we present slow-roll expansions of the matrices.
At early times, k aH, the symmetric kinetic matrix, T, is while the anti-symmetric K matrix has non-zero entries 25) and the symmetric Ω 2 matrix has entries While, for superhorizon modes k aH, (3.28) and and , , (3.30) We note that, while these expressions are extremely accurate in the asymptotic regimes, they are not accurate near horizon crossing, −kτ ∼ 1. Therefore, in order to solve the equations numerically, we are required to use the full expressions presented in appendix B.

Initial conditions and quantization
We set the initial conditions for the fields by canonically quantizing them, and using Bunch-Davies conditions in the asymptotic past. We expand the fields into modes [21] ∆ where we impose the canonical commutation relation between ∆ i and its canonically conjugate momentum .

(3.32)
We decompose the canonical momenta, π i , into the same set of creation/annihilation operators as above, The relations in eqs. (3.31) and (3.32) can only be simultaneously imposed if the condition is obeyed. As pointed out by ref. [21], eq. (3.34) can be imposed as an initial condition, which then holds at all times if the initial conditions satisfy which is equivalent to imposing that the products ππ † and QQ † are real.
In the limit x → ∞, the fields remain coupled, and the separate fields cannot be simply quantized independently. That is, Q ij is not simply proportional to δ ij . Instead, we expand the solutions into normal modes and impose the initial conditions on these solutions to fix the constants. To identify these modes, we use a Wentzel-Kramers-Brillouin (WKB) method.
Working in the limit k aH, and using the expressions for the matrices above, the equations of motion for the fluctuations become Q + αQ + βQ = 0, (3.36) where in the limit k aH, the matrices α and β are given by Adopting a WKB ansatz for the mode functions and substituting into the system of equations, neglecting terms of order O(ω /ω) and O(ω /ω), we find six solutions for the frequencies In order for the system to be stable, all of these instantaneous WKB frequencies must be real. Thus there is an instability in the system for parameters such that γ = g 2 ψ 2 /H 2 < 2, as was found for the original model in ref. [21]. The corresponding mode solutions are, up to an irrelevant phase where the c ij are constants and the a i are the vectors Demanding the solutions approach the positive frequency solutions as x = −kτ → ∞ means we can set c 1j = c 3j = c 5j = 0. The remaining constants now need to be set by imposing the quantization conditions above. Working in the limit −kτ → ∞, it is then straightforward to see that a solution that satisfies the initial conditions is Slow mode: We show the solutions to all three independent modes in figure 2. Notice that the effect of the Higgs vev and accompanying Goldstone fluctuations boosts the final amplitude of the fluctuations. These dynamics are what allows the model to become consistent with the data -the scalar curvature fluctuations are boosted, thus lowering the tensor-to-scalar ratio. In numerically solving the system, we initialize the system including 1/kτ corrections to the solutions described above. This allows more efficient and accurate evaluation starting at later times.

Superhorizon solutions
We can solve the system to a very good approximation in the superhorizon regime, k aH, by expanding the T, K, and Ω 2 matrices in both k/aH and , keeping only the lowest-order non-trivial terms. We also use the relation τ = −(1 + )/aH, which is accurate to first order in .
The matrices at zeroth order in a series in k/aH and lowest non-trivial order in become and Finally, Ω 2 becomes (3.44)

Gauge-flation
We study the M = 0 case-regular Gauge-flation-separately, in order to separate the superhorizon behavior that arises due to the Higgs field from the generic superhorizon evolution. For M = 0 the 3 × 3 system of equations becomes 2 × 2, since the Higgs degree of freedom becomes trivial and decoupled at all times. We look for solutions of the form where ∆ (2) = {∆ 1 , ∆ 2 }, and ∆ (2,0) is a constant two-vector. At late times, the eigenvalues, n, are , leading to the eigenvector In figure 3 we show the late-time (k aH) evolution of the scalar fluctuations. Shown is the numerical solution, as well as the power-law solution from eqs. (3.46) and (3.47) -our analysis accurately captures both the growth rate and the eigenvectors.

Higgsed Gauge-flation
We now move to the full Higgsed Gauge-flation case of M = 0. In this case the Higgs fluctuations are coupled to the other two modes. However, in the superhorizon limit k/aH → 0 the Higgs fluctuation decouples regardless of the mass M , splitting the 3 × 3 system into a 2 × 2 and 1 × 1 system. Setting ξ = ξ 0 (−kτ ) n ξ at late times results in (3.48) The remaining 2 × 2 system for {∆ 1 , ∆ 2 } is solved as before, using an ansatz for the solution that scales as (−kτ ) n . There are four solutions, the three non-growing modes are while the late-time growing mode has the exponent with the corresponding eigenvector (3.51) The excellent agreement of our analytical results with the numerical solution of the corresponding equations is shown in figure 4. It is worth commenting on the superhorizon evolution of the scalar power spectrum P ζ , as defined in appendix D. For the purposes of this section, the exact expression of P ζ is unimportant and we only consider its parametric dependence for superhorizon modes. To first order in slow-roll the Hubble parameter evolves as where we take H * to be the Hubble scale at horizon crossing of a particular mode with comoving wavenumber, k * , i.e. −k * τ * = 1. Furthermore the effective mass parameters γ and M also flow with time. In particular where M * is the value at horizon crossing. The case of γ is in principle more complicated, since both ψ and H evolve in time. However, ψ = ψ * |τ /τ * | δ and since δ = O 2 we can regard it as a constant, leading to Furthermore, η = O( ), which implies that the flow of must also be taken into account.
To lowest order Numerical evaluation of γ and show very good agreement with eq. (3.54) and (3.55), respectively. The late-time power spectrum is dominated by the contribution of M(τ ), therefore we consider only the late-time behavior of its pre-factor. Following the notation of appendix D, the asymptotic behavior is where the second-to-last equality holds for γ + M 2 /2 1. This approximation shows good to excellent agreement with numerical data for all tested values of γ and M .
The parametric time-dependence of the superhorizon scalar power spectrum is Since n scal > 0 for all values of γ > 2 and M > 0, the scalar power spectrum (slowly) decays outside the horizon, |kτ | → 0. This is different to the case of single-field slow-roll inflation, where the time evolution of the prefactor exactly cancels the time evolution of thee superhorizon field fluctuation , so that the scalar power spectrum is exactly constant at late times. Alternatively, the evolution of the curvature perturbation in Gauge-flation indicates the presence of an isocurvature mode, which is absent in the standard single-field scenario. The late-time decay rate of the scalar power spectrum is plotted in figure 5.
The consequence of this decay is that the tensor-to-scalar ratio increases during inflation. Hence, we find that the disagreement of Gauge-flation with Planck data is worse than originally computed in ref. [21]. It is worth noting that the form of the scalar power spectrum given in eq. (3.58) is missing a wavenumber-dependent prefactor. This prefactor arises from the different amount of enhancement that each mode undergoes due to the evolution of the background; each mode sees a slightly different background. Therefore, an estimate of the spectrum tilt n s cannot be read-off immediately from this expression, as it can be in the case of simple single-field inflation, where n s = 1 − 2 − η.

Vector fluctuations
We turn now to the vector fluctuations. Working in fourier space we expand the transverse vector modes into helicity states as After introducing these modes, we find that the action splits into two non-interacting pieces corresponding to the positive and negative helicity states. Explicitly, the vector action reads Figure 5. The late-time decay rate of the scalar power spectrum for different values of γ and M .
where the quadratic Lagrangian density for the vector modes is Note that Y ± and N ± appear in the action without time derivatives, and are thus algebraic constraints. We can thus solve their equations of motion, and insert the solutions back into the action. To leading order in slow roll, the shift vector is given by After substituting back, the resulting action is complicated and not particularly enlightening; we do not reproduce it explicitly here.
Denoting by V ± = (M ± , ξ ± ), we redefine the field using the transformation After making this transformation, the action takes the form Again, it is possible to choose a transformation that sets the matrix T ± to the identity for all times, however, this requires a more complicated transformation which significantly complicates the algebra. As discussed above, for our purposes this is not required. All that is required is that T ± approaches the identity for k aH, so that we can set the initial conditions, and that T ± is invertible for all times, so that we may smoothly evolve the equations of motion.
The full matrices T ± , K ± , and Ω 2 ± are obtained in a straightforward manner, their exact forms are messy and not particularly illuminating. In appendix C we present slowroll expansions of these matrices that are valid at all scales. For use in the subsequent sections, we present expansions of the matrices in the limits k aH and k aH. On subhorizon scales, k aH, the symmetric T ± , and antisymmetric K ± matrices are and the symmetric Ω 2 ± matrix is Analogously to the scalar case, this transformation-which diagonalizes the kinetic matrix at early times-is found by neglecting the gravitational constraints.
On superhorizon scales, these matrices are and (3.70)

Initial conditions and quantization
Following the procedure outlined above in section 3.2.1 for the scalar modes, we quantize the vector modes. We begin by expanding the fluctuations into modes Working in the limit k aH, and using the expressions for the matrices above, the equations of motion for the fluctuations become W ± + β ± W ± = 0, (3.72) where in the limit k aH, the matrix β ± is given by Adopting a WKB ansatz for the mode functions and substituting into the system of equations, neglecting terms of order O(ω /ω) and O(ω /ω), we find four solutions for the frequencies In order for the system to be stable, all of these instantaneous WKB frequencies must be real. Therefore, there is an instability in the system for parameters such that γ = g 2 ψ 2 /H 2 < 1. However, since the scalar sector requires γ > 2, this does not present a further restriction on the model. The corresponding mode solutions are, up to an irrelevant phase, where the c ij are constants and the a i are the vectors In figure 6, we show the evolution of the Higgs vectors and gauge field vectors, as well the evolution of the shift vector.

Superhorizon solutions
We can solve the equations of motion for the vector modes in the superhorizon limit by using the ansatz Using the asymptotic matrices in eqs. (3.68)-(3.70), to first order in , we find (3.80) The growing mode has n = −1 − , and corresponding eigenvector (3.81) In the limit k aH and inserting the late-time solution of eq. (3.81) into eq. (3.82) while keeping the lowest order term in 1/a we arrive at the solution where we used the fact that the scale-factor grows like a ∼ |τ /τ * | −1− during inflation. 8

Tensor fluctuations
We now turn to the tensor degrees of freedom. The addition of the Goldstone modes in the Stueckelberg limit does not add any new degrees of freedom to the tensor sector. At linear order in perturbation theory these fluctuations are gauge invariant under SU (2) transformations and coordinate transformations for the gauge field and metric fluctuations respectively. Furthermore, neither are subject to the Einstein, or Gauss law constraints at this order. We expand the tensor modes into a helicity basis in Fourier space where the polarization tensors satisfy (see, e.g. ref. [53]) 86) 8 We note here that in this model, we find that the vector part of g0i metric perturbation grows exponentially during the inflationary phase. While vector modes typically decay following inflation, this exponentially growing shift potentially invalidates the FRW background during inflation. One can show that the scalar part of the shift vector also generates an exponentially growing g0i perturbation, even in the limit that the Higgs is absent, M → 0. and (±) (k) are the helicity vectors from above. Inserting these decompositions into the action, it splits into two decoupled pieces corresponding to the left-helicity and right-helicity modes. Neglecting boundary terms, the action takes the form where the antisymmetric matrix K ± has entries and the symmetric Ω 2 ± matrix has entries ) , .
In writing these expressions we have made no slow-roll approximation, however, we have made use of the background equations of motion. Note that we recover the results of ref. [21] in the limit that Z 0 → 0, as expected. We now eliminate as many variables as possible.
Making use of the background equations of motion, and expanding to leading order in slow roll, we obtain ) (3.93) Note that in the limit that a → 0, K ±,12 → 0 and Ω 2 ± → k 2 1, which indicates that γ ± and t ± are the canonically normalized fields in the far past. 9 The original model of Gauge-flation is ultimately ruled out, because for a range of k/aH one of the helicitiest ± experiences tachyonic growth due to its mass becoming 9 Note that the recent work of ref. [43] quantized a different combination of γij and tij. In that work the kinetic term is not diagonal in the limit k aH due to couplings proportional to ψ. negative. This leads to exponentially large gauge tensor modes, which in turn source large gravitational waves. Indeed, we see that the effective mass oft − is negative during the range and thus, rather than making the tensors more stable, the effect of the mass in fact makes the instability worse. Examples of the amplification of the tensor modes are shown in figure  7. However, the additional scalar dynamics behave in such a way as to also boost the scalar spectrum, thus reducing the tensor-to-scalar ratio for a fixed value of the Hubble rate.

Sub-horizon solutions and Born approximation
The rather complicated looking system of equations that results can actually be solved analytically to an excellent approximation, as first pointed out in refs. [23,24]. Note that the off-diagonal terms in the K and Ω 2 matrices that couple the gauge field perturbations to the metric fluctuations are slow-roll suppressed by √ . This slow-roll suppression allows us to develop a series solution (in powers of √ ) to the equations of motion using the Born approximation. Further, as in refs. [23,24], the evolution of the gauge modes is dominated by its mass term. The negative helicity mode of the gauge field tensor remains heavy throughout the evolution, and it can be neglected to a good approximation, leaving the negative helicity gravitational wave mode undeflected. We thus focus on the positive helicity modesγ + andt + .
To leading order in slow roll, neglecting interactions with the gravitational wave sector, the equation of motion for the gauge field reads Introducing the variable z = 2ikτ , and the parameters where M α,β (2ix) and W α,β (2ix) are the Whittaker M and W functions. We set the values of the constants A k and B k , by imposing the Bunch-Davies vacuum conditions in the asymptotic past. That is, we demand that the solutions approach canonically normalized positive frequency free plane waves as x = −kτ → ∞, (3.99) In this large x limit, the Whittaker functions have asymptotic expansions, Using eq. (3.100), we find that the constants are given by (3.102) 10 Note that the basis of solutions to the Whittaker equation given bŷ leads to a much simpler set of coefficients A k , B k . However, we choose to work with the basis in (3.98) to make contact with the work in ref. [23,24].
To leading order in slow-roll, the positive-helicity tensor modes of the metric obey the equation of motion where x = −kτ . Equation (3.103) can be solved as a series in √ using the Born approximation. This solution consists of a homogeneous piece that solves the free equation of motion, and an inhomogeneous piece that is sourced by the gauge field fluctuationŝ Here,γ + 0 is the homogeneous solution that matches onto the Bunch-Davies vacuum as x = −kτ → ∞. To leading order in √ , the inhomogeneous part of the solution iŝ is the Green's function, and W [. . .] is the Wronskian.
We can proceed in exact accordance to the analysis found in refs. [23,26] for the physically and mathematically related models of Chromo-Natural and Higgsed Chromo-Natural Inflation. The positive-helicity gravitational wave undergoes exponential amplification and the late-time solution is well-approximated by where the solution is written as the sum of a free and a sourced part, arising from the interaction with the gauge field fluctuations. The functions I 1 , I 2 , I 3 arise from the three distinct integrals of eq. (3.105) The evolution for the negative helicity gravitational wave is largely independent of the potential parameters and follows closely the free expanding-universe solution The total power spectrum at late times is given by the sum of the power in the positiveand negative-helicity modes where we define the power in the each mode as Since the negative-helicity modes follow the free-field result, their power is (3.112) while the power in the positive-helicity modes can be written as because the free and sourced parts are uncorrelated. It is now straightforward to compute the chirality parameter ∆χ = which is plotted in figure 8. We see that the gravitational wave spectrum is strongly polarized for most of the parameter-space plotted, especially for larger values of M or γ.

Superhorizon solutions
On superhorizon scales corresponding to k aH, we can also solve for the evolution of the system very accurately. In this limit, the parity violating terms become irrelevant, and thus we can ignore the difference between the positive-and negative-helicity modes. In this limit, We then look for a solution of the form where ∆ 0 is a constant vector. Substituting into the asymptotic mode equation yields the four solutions for n (two for each independent equation of motion) (3.118) Note that the only growing mode is the solution with n = −1 − , which corresponds to the amplitude ratio Therefore, with the scale factor given by eq. (A.2), on superhorizon scales the gravitational wave spectrum becomes constant. Despite the appearance of an apparent mass for the graviton, the superhorizon evolution is such that the resulting gravitational waves become constant. 11 While it appears that the tensor tilt can be read off eq. (3.119) to be given by the standard n t = −2 , in the case at hand, the prefactor A ± has significant scale dependence, in contrast to standard single-field slow-roll inflation. We demonstrate below that the tensor spectral tilt, n T , can take both negative as well as positive values, corresponding to a red or blue spectrum, respectively.

Phenomenology
In this section we explore the observational consequences of Higgsed Gauge-flation. We numerically compute both the positive-and negative-helicity gravitational wave amplitudes γ ± as well as the density fluctuation ζ. 12 As described above in section 2.1, in the slow roll approximation the evolution of the background can be completely specified by the values of {γ in , M in } N in e-folds before the end of inflation. In what follows, we parameterize the resulting spectra in terms of these parameters. The value of the parameter κ is then determined by matching the amplitude of the scalar spectrum to the observed value.
In order to solve the equations of motion for the fluctuations numerically, we work with e-folding number N as our time variable. For the background, treating ψ ≈ constant, we use the analytical expressions for {γ(N ), M (N ), (N )} as described above in section 2.1. This method allows for a straightforward computation of the observables closer to the end of inflation, 60 e-folds after horizon-crossing. 11 In the parameterization of ref. [43] this superhorizon solution corresponds to the vanishing of the 'genuine tensor perturbation of the gauge field.' In that work, the slow-roll suppressed backreaction of the gravitational wave modes onto the gauge fields was dropped. Generically, this backreaction sources a growing mode of the gauge field. However, we demonstrate here that the solution corresponding to the adiabatic mode corresponds to a solution where these gauge modes are absent. 12 We verified through explicit computation that our results are unchanged if we instead work with the comoving curvature perturbation in spatially-flat gauge R = Hδu, and on superhorizon scales k aH, |ζ| = |R|.
As a check for accuracy, we performed the same computation in conformal time τ , simultaneously solving the background and fluctuation equations numerically. We initialize the system close to the slow-roll attractor using the following procedure: , for chosen values of γ in and M in at N in = 60 e-folds before the end of inflation.
3. Finally, the initial velocity of the ψ field isψ in These initial conditions start the calculation close to the exact numerical solution for all values of γ and M . We found excellent agreement between these methods (evolution in conformal time vs e-folding number) in both the evolution of each mode-functions, as well as the resulting values of r, n s , and n T . As shown in sections 3.2.2 and 3.4.2 respectively, the scalar modes evolve outside the horizon, while the tensor modes become constant. This means that the tensor-to-scalar ratio is not constant, but rather evolves from horizon exit until the end of inflation. Even though the evolution is rather weak, since it involves the slow-roll parameter = O(0.01), a simple order-of-magnitude calculation shows that the tensor-to-scalar ratio evolves as (τ comp /τ end ) O( ) = O(10), where τ comp is the conformal time where we compute r and τ end is the end of inflation. Thus in the 60 e-folds between horizon-crossing and the end of inflation the tensor-to-scalar ratio can vary by an O(10) factor. For this reason, we compute the observables at 5 e-folds before the end of inflation. We did not choose to evolve further, in order to ensure the accuracy of the slow-roll formulas. At N = 5 e-folds before the end of inflation the numerically computed (exact) solution for and the slow-roll expression start to deviate. Additionally, at this point 0.1, and the -expansion of the equations of motion begins to break down after this point.
Beyond the tensor-to-scalar ratio, the superhorizon evolution of the scalar power complicates the numerical computation of the spectral index n s . For the standard single-field slow-roll inflationary models, the power spectrum freezes outside the horizon and the spectral index can be evaluated by comparing P ζ (k) for neighboring values of the wavenumber k ± ∆k at some constant value of conformal time τ , e-folding number N , or x = −kτ , provided the mode has left the horizon. 13 In the present case the first two methods (constant τ or N ) are equivalent sufficiently far outside the horizon, but the case of taking x = −kτ constant is not. The value of n s computed using the two (inequivalent) methods can differ by O( ) = O(0.01). We choose to evaluate the spectral index by comparing the power-spectra at a fixed number of e-folds before the end of inflation. Even though the spectral index settles to a constant value after around 10 e-folds after horizon-crossing, even for the largest values of M used, we evaluate it at 5 e-folds before the end of inflation, due to the evolution of the tensor-to-scalar ratio as explained above. Figure 9. The spectral index n s (left) and the tensor-to-scalar ratio r (right) as a function of the mass M for γ = 4, 5, 6, 7, 8, 9, color-coded from red to green. The horizontal black-dotted lines correspond to the Planck 2-σ bounds.
The evolution of the scalar power in the last 5 e-folds is a source of uncertainty in our results. One the one hand, the continued decay of the scalar power likely increases the tensor-to-scalar ratio above the values we quote. On the other hand, since all modes evolve identically outside the horizon, the scalar spectral index is frozen and does not evolve during this period.

Spectral tilt and tensor-to-scalar ratio
For the case Gauge-flation, M = 0, we recover the results of ref. [21] by computing the spectral index at constant x = −kτ , or 3 e-folds after each mode has left the horizon. By computing n s at constant N we get the same dependence of n s on γ, however, the spectral index is larger by about ∆n s = 0.015 compared to the results of ref. [21]. This is a change of O( ) as expected. However, the tensor-to-scalar ratio is larger by a factor of 10 when calculated close to the end of inflation (N = 5) instead of close to horizon crossing (3 e-folds after horizon-crossing as in ref. [21]), due to the superhorizon decay of the scalar power spectrum, as explained in section 3.2. This pushes the observables of Gauge-flation even further from the Planck constraints.
As the Higgs vev, and thus M , is increased, both the tensor-to-scalar ratio, r, and the scalar spectral index, n s , develop a distinct spike for any value of γ, which is however significantly more pronounce for smaller values of γ. This feature is shown in figure 9. In order to understand the feature in r and n s , we revisit the behavior of the evolution of the scalar and tensor power spectra for M = 0. In figure 10 we show the evolution of the scalar power spectrum and the amplitude of mode that contributes dominantly to the scalar spectrum, M(−kτ ), for constant γ and various values of M bracketing the "spike". Figure 10 shows that at fixed γ and increasing M , the power in the scalar modes initially decreases before reaching a minima, and then ultimately increases quickly. In contrast, the power in tensor modes increases monotonically for constant γ and increasing M . The "spike" in the spectral parameters is due to the non-monotonic behavior of the scalar spectrum as M is varied; the tensor to scalar ratio reaches its maximum value where the scalar spectrum is minimized. Further, across this minima the spectral tilt goes from blue to red. This is counter-intuitive, since M increases during inflation, one might naively expect the opposite. However, as demonstrated below, γ and its evolution have a much stronger effect on the tilt of the spectrum (see figure 11). We note that this behavior was not seen in the related model of Higgsed Chromo-Natural Inflation [26]. However, in that work the region of parameter space corresponding to the (Higgsed) Gauge-flation model was not explored. For large-enough values of M , the observables of Higgsed Gauge-flation pass through the Planck-allowed region in the n s -r plane, as shown in figure 11. By varying both γ and M we can fill-up the whole of the allowed Planck region 14 in the n s -r plane for r 10 −5 .

Running of the scalar spectral index and inflationary energy scale
We compute the running of the scalar spectral tilt, α ≡ d n s /d log k, by locally fitting log P ζ as a function of log k using a second order polynomial around the mode-function k that leaves the horizon 60 e-folds before the end of inflation. The results for the parameters that correspond to the black dots of figure 11 are shown in figure 12. We see that for all but one of the cases shown the running is positive and also for all but one of the cases shown, the running is within the Planck limits.
14 It is worth re-iterating here that we do not compute the evolution of the fluctuations during the last 5 e-folds of inflation. Since we do not expect ns to vary, these last 5 e-folds may shift the curves of figure 11 vertically, hence filing a slightly different part of the Planck-allowed region. Figure 11.
The tensor-to-scalar ratio, r, (at k = 0.002h Mpc −1 ) as a function of the scalar spectral tilt, n s (at k = 0.05h Mpc −1 ) and the tensor spectral tilt, n T (at k = 0.002h Mpc −1 ) for models drawn from a grid of values for the parameters γ and M measured at N = 60 e-folds before the end of inflation. The mode k = 0.05h Mpc −1 is assumed to leave the horizon 60 e-folds before the end of inflation. In both panels each rainbow-colored line corresponds to a definite value of γ ranging from γ = 5 (red) to γ = 9 (green). Each purple-blue colored line corresponds to a definite value of M , ranging from M = 4 (far left purple-dotted line) to M = 7 (far right cyan-dotted line). The dots correspond to the values used in Table 1. The dotted purple-to-blue lines shown on the left panel do not appear on the right one, since they are either largely outside of the Planck-allowed regime (purple-dotted) or lead to a very blue tensor spectrum (cyan-dotted). In both panels, the shaded light red region correspond to the 1% limit of the linear regime, as discussed in section 4.1.
Until now, we have rescaled the gauge coupling, g, and cosmic time, t, by κ, eliminating it from the equations of motion as described above in section 2.1. However, the value of κ sets the Hubble scale, and the overall energy scale of inflation. We fix κ by matching the the amplitude of the observed scalar power spectrum P ζ 2.2 × 10 −9 [28]. Table 1 shows the required value of κ and the corresponding value of the Hubble scale. We note here that the standard inflationary result (see, for example, [12]) is violated.

Validity of the linear theory
We end this section with a brief discussion of the validity of our analysis. As in the case of Higgsed Chromo-Natural Inflation [26], the introduction of a Higgs sector enhances the tensor mode amplification, rather than suppressing it. Furthermore, the Goldstone mode dynamics contribute constructively to the amplification of the scalar and vector modes. It is thus of paramount importance to examine whether the use of the linear order perturbation theory and of the scalar-vector-tensor decomposition of the fluctuations is justified. Following the analysis of ref. [26], we provide an estimate of the non-linearity, rather than attempt a detailed analysis.
The gauge field mode amplification begins around −kτ ∼ M and ceases around the time when the mode exits the horizon at −kτ = 1, as shown for scalar, vector and tensor modes in figures 2, 6, and 7 respectively. In order to keep the linearized analysis under control, it is sufficient for the backround field fluctuations δA µ to be significantly smaller than the classical (background) gauge field value, which is given in eq. (2.5) asĀ µ = The gauge (solid) and metric (dotted) tensor modes, using the same color-coding. The gauge modes are multiplied by 0.1, in order to make the gauge and metric amplitude easier to compare visually. We show only three out of the six curves for clarity.
We consider this ratio as a measure of the validity of the linear analysis. We focus on the tensor part of the gauge fluctuations. In order to compute the ratio of eq. (4.2), we estimate the gauge field fluctuations as where we cut the integral off at the peak of the amplification, which is observed to be near −kτ ∼ M . The linearity condition of eq. (4.2) can thus be rewritten as where we used the definition of γ to swap H for g.
In figure 13, we plot the left hand side of eq. (4.4). Note that most of the cases of interest satisfy the linearity criterion at the 1% level or better. As in the case of Higgsed Chromo-Natural Inflation, the theory remains within the linear regime for lower values of the tensor-to-scalar ratio. A simple estimate of this can be done in a more general way using dimensional arguments. The background amplitude of the gauge field depends only on the number of e-folds, which we take to be N = 60, and the values of where we used the fact that P T = rP ζ and P ζ ≈ 2.2 × 10 −9 [28]. The advantage of this form of the linearity condition is that it involves only the Hubble scale and the tensorto-scalar ratio. Finally, we emphasize that both eqs. (4.4) and (4.5) should be used as order-of-magnitude estimates of the non-linearity, rather than a sharp cut-off. Before concluding this section, we revisit and justify the claim made throughout this work regarding the dominance of the contribution of the 'slow' scalar mode compared to the 'regular' and 'Higgs' modes in the scalar spectrum (see eq. (3.41)). Figure 14 shows the ratio of the final scalar power spectrum computed using the three initial conditions, slow, regular and Higgs. In the region of parameter space where the spectral tilt lies within the Planck bounds, M 5 (see figure 9), the regular and Higgs modes contribute negligibly to the final power spectrum. We are therefore justified in only using the "slow" mode to initialize the numerical evolution of the fluctuations in this region. Away from this region of parameter space, M 3, the otherwise subdominant modes can dominate the spectrum, however, in this region of parameter space the model is ruled out by the large tensor-to-scalar ratio, as demonstrated in figure 9.

Discussion and conclusions
In this work we have demonstrated that Gauge-flation can be made compatible with existing limits from Planck data by introducing an additional mass term for the gauge field fluctuations. We assume that the symmetry is spontaneously broken by a Higgs sector and the resulting Higgs boson is much heavier than the Hubble scale, and is thus irrelevant. We thus work with the theory in the Stueckelberg form, restricting the Higgs sector to the Goldstone modes which fluctuate along the vacuum manifold.
The introduction of a Higgs sector to Gauge-flation significantly changes the phenomenology of the fluctuations. On the one hand, the additional terms in the tensor action have the effect of exacerbating the existing chiral gauge-tensor instability, and boost the overall amplitude of the tensor modes. The resulting chiral spectrum of gravitational waves becomes constant on superhorizon scales, and is qualitatively similar to the spectrum of gravitational waves from Gauge-flation. On the other hand, the mass terms and associated scalar Goldstone modes significantly alter the scalar dynamics and resulting spectrum of density fluctuations. The theory remains (catastrophically) unstable for parameters such that γ < 2. However, as the Higgs vev, Z 0 , and thus M , is increased from zero (with γ fixed) the amplitude of the scalar power is initially suppressed before reaching a minima. This behavior leads to a feature in both the tensor-to-scalar ratio, r, and the spectral index, n s , at the points in parameter space where the scalar power reaches its minimum. By increasing the Higgs vev, and thus M , beyond this point, the overall amplitude of the resulting scalar density fluctuations grows monotonically, decreasing the tensor-to-scalar ratio r, and increasing the scalar spectral index n s . This behavior allows the theory to produce spectra of gravitational waves and density fluctuations that satisfy the latest Planck bounds. The scalar spectrum decays on superhorizon scales in this model due to the presence of an isocurvature mode. We evaluate the tensor-to-scalar ratio and scalar spectral index 5 e-folds before the end of inflation to minimize errors, however, the evolution of the spectrum during the last 5 e-folds of inflation represents a source of uncertainty in our results.
The exponential enhancement of the tensor modes means that observable gravitational waves may be produced in this model, despite inflation occurring below the GUT scale, and all fields evolving over sub-Planckian distances in field space. The model therefore violates some formulations of the Lyth bound. The production mechanism of gravitational waves is exactly analogous to the cases of Chromo-Natural Inflation, Higgsed Chromo-Natural Inflation, and Gauge-flation. The gravitational waves in these models predominantly arise from linear mixing with the (exponentially amplified) gauge field fluctuations. The form of the gravitational wave spectra produced in this model is therefore significantly altered from the usual form assumed in formulations of the Lyth bound. In contrast to standard inflationary scenarios which uniformly predict red tilted gravitational wave spectra, these gravitational waves can have either red-or blue-tilted spectra on CMB scales, with a strong favoring of a blue tilt, especially for lower values of the tensor-to-scalar ratio r. Furthermore, these gravitational waves have the distinct characteristic that they are chirally polarized and, to a very good approximation, consist only of a single helicity. Unfortunately, it seems unlikely that future CMB experiments will be unable to distinguish between unpolarized and chirally polarized gravitational waves [54,55], which would significantly reduce the space of viable inflationary models.
The running of the scalar spectral index is predominately positive and within observational bounds for all but the lowest calculated values of the tensor-to-scalar ratio. Between the linearity constraints for the equations of motion of the fluctuations and the Planck constraints on the running of the spectral index, Higgsed Gauge-flation can fill the whole Planck-allowed region on the n s -r plane for 10 −4 r 10 −2 , making this model especially interesting in anticipation of planned Stage-4 CMB experiments. These experiments are aiming to probe tensor-to-scalar ratios as low as r ∼ 10 −3 . The exponential sensitivity of the amplification of both the scalar and tensor power spectra makes some level of finetuning necessary to fit observations. In contrast to Gauge-flation, the addition of the Higgs sector causes vector perturbations of the matter sector to freeze out on super-horizon scales, we leave the further study of the consequences of these modes to future work. While we have estimated that the linear theory is under control, we leave the study of non-Gaussian features of this model for future work.
respect to conformal x are kept explicit (∂ x ). Our symmetrization and antisymmetrization conventions throughout are

B Details of the scalar action
In this appendix we present the details of the scalar action. After eliminating the algebraic constraints from the action, redefining the fields according to eq. (3.22), performing an integration by parts and discarding the boundary terms, the scalar action is put into the form The exact forms of the matrices can be obtained in a straightforward manner, however, they are long and complicated, and not particularly enlightening. We do not present their gory details here.

C Details of the vector action
In this appendix, we present the details of the vector action. After making the transformation in eq. (3.66), the action takes the form While it is straightforward to obtain the matrices exactly, they are long and not particularly enlightening. Each entry in the matrices are of the form of eq. (B.2). To perform the numerical evaluation, we expand each coefficient to leading order in slow-roll in the same way as described in section B. We find Finally, the entries of the mass matrix take the form Ω 2