Tensor Non-Gaussianity from Axion-Gauge-Fields Dynamics : Parameter Search

We calculate the bispectrum of scale-invariant tensor modes sourced by spectator SU(2) gauge fields during inflation in a model containing a scalar inflaton, a pseudoscalar axion and SU(2) gauge fields. A large bispectrum is generated in this model at tree-level as the gauge fields contain a tensor degree of freedom, and its production is dominated by self-coupling of the gauge fields. This is a unique feature of non-Abelian gauge theory. The shape of the tensor bispectrum is approximately an equilateral shape for $3\lesssim m_Q\lesssim 4$, where $m_Q$ is an effective dimensionless mass of the SU(2) field normalised by the Hubble expansion rate during inflation. The amplitude of non-Gaussianity of the tensor modes, characterised by the ratio $B_h/P^2_h$, is inversely proportional to the energy density fraction of the gauge field. This ratio can be much greater than unity, whereas the ratio from the vacuum fluctuation of the metric is of order unity. The bispectrum is effective at constraining large $m_Q$ regions of the parameter space, whereas the power spectrum constrains small $m_Q$ regions.


Introduction
We do not yet know how to quantize gravity over an entire spacetime, but we can quantize its perturbations around a specified background. In this case, the degrees of freedom of the gravitational field, including the transverse traceless tensor modes of the metric, should have ground state vacuum fluctuations [1,2]. These can be found in B-mode polarisation of the cosmic microwave background (CMB) [3,4]; thus, a detection of non-zero (primordial) B-modes is evidence for tensor fluctuations of the metric.
So far, no such evidence has been found [5]. CMB experiments provide constraints on the tensor-to-scalar ratio r, which is defined as the ratio of the power in tensor modes (P h (k 0 )) to the power in scalar modes (P ζ (k 0 )), at some wavenumber k 0 , r ≡ P h (k 0 )/P ζ (k 0 ). Currently, this ratio is constrained to be r < 0.07 (95 % C.L.) [5] at k 0 = 0.05 Mpc −1 .
We stress here that a detection of B-mode polarisation in the CMB is evidence for primordial tensor perturbations, but is not necessarily evidence for the vacuum fluctuation in the tensor metric. For the standard scenario of single-field slow-roll inflation, the tensor fluctuations of the metric, h ij , obey the equation where is the d'Alembertian operator in 4-dimensions. Equation (1.1) shows that if we find evidence for tensor fluctuations, in the absence of anything that can source them, they have to be necessarily quantum. However, there is no a priori reason to ignore sources in the right hand side of Equation (1.1). It is reasonable to think that there are (many) more than one field during inflation. While their energy density may be much smaller than that of the dominant inflaton field, they can still act as sources of perturbations. In general, we write h ij (t, x) = Π ij (t, x) , (1.2) and tensor perturbations are sourced by the anisotropic stress-energy Π ij , which is provided by quantum fluctuations of a field other than the metric. These sourced tensor fluctuations can be much larger than the vacuum one and can generate observable B-modes, invalidating the claim that B-modes are evidence for vacuum fluctuations of the metric. Consequently, there have been intense efforts to build inflation models where a sizeable r can be generated from sources, without violating stringent observational constraints on the scalar perturbation.
The sources include scalars [6][7][8][9], U(1) gauge fields [10][11][12][13][14], and SU(2) gauge fields [15][16][17][18][19][20][21][22]. How then, do we differentiate between B-modes generated from vacuum fluctuations of the metric and those from sources? Vacuum fluctuations of the metric are usually almost scale invariant, with a slightly red tilt (see [23] for the latest review). On the other hand, B-modes from sources can have a red or blue tilt or completely non-power-law spectra such as bumps, depending upon model parameters. Moreover, the tensor fluctuations produced by sources can be chiral (see section 3), and so can be seen as a non-vanishing TB/EB correlation in the CMB [24][25][26][27][28][29], whereas vacuum fluctuations produce parity-even B-modes. Finally, these modes can be highly non-Gaussian [30]. Because tensor modes from vacuum fluctuations of the metric are almost Gaussian [31,32], non-Gaussianity provides strong evidence for sourced tensor modes. Therefore, we hope that in any future detections of primordial gravitational waves (GWs), one would not only check for their amplitude (r) and scale-invariance, but also their non-Gaussianity and parity-violating correlations. Only after all these tests support nearly scale-invariant, Gaussian, and non-chiral primordial B-modes, can we confidently claim to have discovered vacuum fluctuations of the metric.
We focused on a particular set of model parameters in ref. [30], using a model proposed by Dimastrogiovanni, Fasiello and Fujita [20]. In this paper we shall give more detailed derivations of the bispectrum and present the results for wider parameter space. The rest of the paper is organised as follows : in section 2 we present details of the model that we consider. In section 3 we present the second-order Lagrangian for the tensor perturbations in our model and their imprint on the B-mode power spectrum. The third-order Lagrangian is presented in section 4 and is used to calculate the bispectrum of metric fluctuations. A detailed discussion of the deviation of the bispectrum from the equilateral shape is given in section 5. In section 6 we explore parameter regions of the model, which can be potentially observed in upcoming CMB missions. We conclude in section 7.

Model Setup
In the model of ref. [20], inflation is driven by a scalar inflaton φ, which is only minimally coupled to a pseudoscalar axion χ and SU(2) gauge fields, A a µ . The SU(2) gauge fields and the axion have negligible energy densities compared to the inflaton, and thus are called "spectator fields". They are coupled to each other by a Chern-Simons like interaction χFF . The Lagrangian is then given as, L = L GR + L φ + L spec , (2.1) where Einstein gravity L GR = M 2 P R/2 is assumed, the Lagrangian of the inflaton L φ is not specified, and L spec denotes the Lagrangian of the spectator fields. V (χ) is the potential of the axion field with the canonical kinetic term −(∂χ) 2 /2, the gauge field strength tensor F a µν is written in terms of the gauge fields as F a µν = ∂ µ A a ν − ∂ ν A a µ − g abc A b µ A c ν with g being the self-coupling constant, a dimensionless parameter λ controls the strength of the Chern-Simons interaction, f is a decay constant of the axion field, andF aµν ≡ µνρσ F a ρσ /(2 In the rest of this section, we discuss the background dynamics, while perturbations will be studied in the following sections. In this paper, we do not solve for the inflaton φ(t), but consider dynamics of the spectator fields in a de Sitter universe, where the Hubble expansion rate is constant. We also leave the axion potential V (χ) unspecified by assuming that it supports slow-roll of the background axion χ 0 (t) with the aid of the coupling to the SU(2) fields. While these assumptions are far from generic, they still capture the essence of physics of generation of non-Gaussianity, and are observationally relevant because they produce scale-invariant GWs.
As shown in [33,34], while the background axion slowly evolves, the homogeneous background component of the gauge fields has an attractor configuration which respects isotropy of the universe, where a(t) is the scale factor. Then we decompose these spectator fields into the background and the perturbation components as There also exist non-dynamical components δA a 0 that we integrate out. The equations of motion (EoM) for the background fields are given bÿ where the dots denote cosmic time derivatives ∂ t and H ≡ȧ/a is the Hubble expansion rate. The terms on the right hand side of eq. (2.5) slow down the time evolution of χ 0 (t) in addition to the Hubble friction term 3Hχ 0 , because a non-zero background value Q(t) is sustained by energy transfer from the kinetic energy of χ 0 through the coupling. Here we introduce two dimensionless parameters; Here, m Q is the effective mass of the SU(2) field around its vacuum expectation value (vev) normalized by the Hubble scale, and Λ characterizes the coupling strength between χ 0 and Q.
Note that the right hand side of eqs. (2.5) and (2.6) are proportional to m Q Λ. We consider the slow-roll regime, m Q 1 and Λ 1, in which Q is stabilized by its effective mass and χ 0 is significantly slowed down by the coupling. We can then drop all the terms with time derivatives in the EoMs except for the r.h.s. of eq. (2.6) and find [33] The Einstein equations at the background yield where ρ φ is the energy density of the inflaton and the slow-roll parameters are defined as We shall assume that ρ φ dominates in eq. (2.10

Amplification of Gravitational Waves
In this section, we study the tensor perturbations at linear level. Only the tensor perturbations are amplified due to tachyonic instability, while scalar and vector perturbations are not amplified for m Q > √ 2 in this model [35]. Although the scalar perturbations of χ and A a i are generated from the vacuum fluctuations, their contribution to the curvature perturbation ζ is negligible, unless the energy density of χ becomes comparable to that of the inflaton after inflation [20]. 2 The vector perturbations decay on super-horizon scales in any case. Therefore in this paper, we assume that the observed curvature perturbation was produced from the inflaton fluctuation δφ and concentrate on the tensor perturbations from the gauge fields. 1 In ref. [20], the background dynamics is numerically solved, and the perturbations are also solved with the time varying background quantities, mQ(t), B (t) and H(t). They find that mQ, B , H ≈ const. is a very good approximation for a sufficiently strong coupling, especially when one is interested in the range of wavenumbers observable by CMB. 2 The scalar perturbations of χ and A a i directly contribute to ζ through the density perturbation (e.g. δρχ ∂χV δχ). This channel is negligible, for instance, if χ reaches its potential minimum (i.e. V (χ), ∂χV (χ) → 0) during inflation. If χ0 acquires a non-negligible energy fraction after inflation, however, the contribution to ζ from δχ may be relevant. This implies that the spectator sector can produce ζ in a way similar to the curvaton mechanism [36,37]. We leave this intriguing possibility for future work.
To calculate the power spectrum and bispectrum of GWs we need to expand the action, equation (2.1), up to second and third order in perturbations, respectively. We write the tensor perturbations of the metric and the SU(2) gauge field [17,35] as: where · · · represents the scalar and vector perturbations of the SU(2) gauge fields which we neglect in this paper. We have imposed the transverse and traceless conditions on h ij and t ij , δ ij T ij = ∂ i T ij = ∂ j T ij = 0 (T = h and t). The inverse metric is given by . For later convenience, we redefine h ij as Precisely speaking t ai is not a tensor, since the index a is not a spatial index but the label of SU(2) gauge. Nonetheless, under the background configuration of eq. (2.3), t ai transforms as a tensor in practice, because the gauge index a is identified with a spatial index. Substituting these in equation (2.1), and expanding up to third order, we obtain the Lagrangian of the tensor perturbations as where τ −1/aH is the conformal time, prime denotes the derivative with respect to τ and we neglect terms suppressed by slow-roll parameters. The cubic Lagrangian L 3 will be discussed in the next section.
From the quadratic Lagrangian above, we obtain the following EoMs for GWs ψ ij (τ, x) and tensor perturbations of the SU(2) gauge field t ij (τ, x), Although terms linear in ψ ij also source the gauge field t ij , ψ ij is not as substantially amplified as t ij [17][18][19][20], and so we ignore its contribution as a source for t ij . To solve the dynamics of ψ ij and t ij , it is useful to decompose them with the circular polarisation tensors, where X = ψ, t and the properties of the polarisation tensors are summarized in appendix. A. Note that we normalize e p ij such that e R ij (k)e R ij (−k) = e L ij (k)e L ij (−k) = 1.
To proceed further, we quantize ψ and t and expand them in a perturbative series as [39] X p k (τ ) =X p 1 (τ, k) +X p 2 (τ, k) + . . . . (3.8) The first order components are written aŝ with the creation/annihilation operators,â p k andâ p † k , satisfying [â p k ,â q † −k ] = (2π) 3 δ pq δ(k+k ). We only consider GWs sourced by the gauge field in this paper, and assignψ 1 the same quantum operator ast 1 . The mode functions ofX p 1 satisfy linearised EoMs and their solutions induce the second order fieldsX p 2 through non-linear terms in the EoMs. In Fourier space the EoMs for the linear mode functions can be written as, where x ≡ −kτ . The minus and plus signs are for right-(R) and left-handed (L) modes, respectively. Only T R 1 undergoes an instability and it sources only Ψ R 1 [15,18,19,33], provided that m Q is positive. Therefore, we only consider the right-handed polarisation in the rest of the paper. The homogeneous solution for T R 1 can be analytically calculated and is expressed in terms of the Whittaker function W β,α (z) as where α ≡ −i 2m 2 Q + 7/4 and β ≡ −i(2m Q + m −1 Q ) [15,17,20]. Ψ R 1 can then be calculated using Green's function method, with where Θ(x) is the unit Heaviside function of x. In figure 1, we plot T R 1 , the source term DT R 1 , Green's function G ψ in the super-horizon limit and the sourced gravitational wave Ψ R 1 . The time integral of the source term DT R 1 multiplied by G ψ yields Ψ R 1 . The source term DT R 1 peaks around the horizon crossing and Green's function stops oscillating there. As a result, Ψ R 1 is mainly produced around the horizon crossing as well, as seen in the right panel.
Equation (3.14) can also be analytically solved and in the super-horizon limit we obtain, from which we obtain the power spectrum of h in the super-horizon limit [20] where F is a function of m Q , whose exact expression is given in [20] (F(m Q ) here is The function F(m Q ) can be approximated to be an exponential function of m Q , |F(m Q )| ≈ e 1.8m Q . Note that eq. (3.18) is derived under the assumption of (m Q , B , H) = const. However, as long as the time variations of these background quantities are slow, eq. (3.18) with m Q (t), B (t), H(t) at the horizon crossing time k = a(t)H(t) gives a good approximation of P sourced h (k). It should be also noted that only the right-handed polarisation modes contribute to the above P sourced h (k).

Bispectrum of Gravitational Waves
In this section, we calculate the tensor bispectrum of the right-handed GWs, B RRR h , in the super-horizon limit: The three-point correlator of the right-handed GWsψ R =ψ R 1 +ψ R 2 can be written as becauseψ R 1 satisfies Gaussian statistics. We calculateψ 2 using the second order EoMs for the tensor perturbations which are derived from the cubic Lagrangian.
The cubic tensor Lagrangians introduced in eq. (3.3) are given by where we organize terms such that L = O(ψ 2 t) and we neglect the O(ψ 3 ) terms. The coefficients of the cubic Lagrangians are They satisfy a hierarchical relationship, The tree-level contributions from L to the tensor bispectrum are illustrated as Feynman diagrams in figure 2. As we see below, the contributions from the three diagrams to the gravitational wave bispectrum are also hierarchical, 3 In what follows, we calculate these three contributions in order. When we show plots in this section, we use the following parameters as an example, The viable parameter space will be explored in section. 6.

Diagram (i)
This diagram arises from the self-interaction of the SU(2) gauge field, and thus is absent in Abelian theory. Here, the second order gravitational waveψ R 2 is sourced linearly byt R 2 , but the second order gauge field perturbationt R 2 is produced byt R 1 via non-linearity. The where the source term S (i) ij is evaluated with the first order solutiont R 1 (τ, k). Although in second order it is no longer true that the right-handed polarisation is sourced only by the right-handed tensors, the exponential amplification of the right-handed modes ensures that terms containing the left-handed modes are exponentially smaller. Thus we only use the right-handed polarisation of the gauge field perturbation to evaluate the source term. Note that this source term contains g explicitly because of the non-Abelian nature of the vertex.
Expanding the tensor perturbation with the tensor polarisation as before, one finds the EoM in Fourier space aŝ Using the homogeneous solution, equation (3.13), Green's function for equation (4.11) can be written as [39], where Im(z) denotes an imaginary part of a complex number z. Dependence of Green's function on the homogeneous solution also ensures that the second order left-handed polarisation of the gauge field is sub-dominant, even if sourced by the first order right-handed polarisation.
Then, the second order gauge field is given aŝ which yields the second order sourced metric perturbation (c.f. equation (3.8)) aŝ Substitutingψ R 1 and the above expression forψ R 2 into eq. (4.2), we obtain As discussed in appendix. A, contraction of the polarisation tensors is calculated as where we have defined with r 2 ≡ k 2 /k 1 and r 3 = k 3 /k 1 . Using this, we obtain Since we are interested in the bispectrum in the super-horizon limit kτ → 0, Green's function G ψ (τ, η 3 , k 3 ) can be reduced. By changing the integration variables from η 3 and η to y ≡ −k 1 η 3 and z ≡ −k 1 η, we obtain starts undergoing a tachyonic instability, to avoid incorporating unphysical vacuum contributions. The integration result is not sensitive to the precise value of the cutoff, as we checked numerically for several different values of x max . Using the super horizon solution for Ψ R 1 , equation (3.17), and ψ ij = aM P h ij /2 we finally obtain The factor of δ D (k 1 + k 2 + k 3 ) ensures the triangle condition, namely that the three wave vectors k 1 , k 2 , and k 3 form a closed triangle. This is a consequence of homogeneity and isotropy of the Universe. In the same way, the other two terms in eq. (4.2) can be calculated. Combining them, we obtain the contribution from the diagram (i) as [30], where W 1 (z) = W β,α (−2ir 2 z)W β,α (−2ir 3 z), W 2 (z) = W * β,α (−2iz)W β,α (−2ir 3 z), and W 3 (z) = W * β,α (−2ir 2 z)W * β,α (−2iz). Figure 3 shows the tensor bispectrum from the diagram (i). We shall discuss the shape of the bispectrum in detail in section 5.

Diagram (ii)
The bispectrum from the diagram (ii) is generated by a non-linear sourcing of the gravitational wave via the O(ψtt) terms in the Lagrangian, L 3 . The EoM for the second order gravitational wave with the corresponding source terms in Fourier space is given by where the source term is written as the sum of two parts, 1ij has terms without time derivatives, The second order gravitational wave is given using Green's function, We can now compute the first term in eq. (4.2) produced via the diagram (ii) as The expectation values of the 4-point functions are given by . Note that these functions are invariant under interchange of q 1 ↔ q 2 . As a result, upon integrating the Dirac delta functions, the polarisation factors in equation (4.30) yield where Ξ ≡ Ξ(r 2 , r 3 ) has been defined in equation (4.19). Substituting this in equation (4.30), we obtain We can similarly evaluate the other two terms in eq. (4.2) contributed from L (ii) 3 . Taking the super-horizon limit, we obtain [30], whereW 1 (y) = ∂ y W β,α (−2ir 2 y)∂ y W β,α (−2ir 3 y),W 2 (y) = ∂ y W * β,α (−2iy)∂ y W β,α (−2ir 3 y), W 3 (y) = ∂ y W * β,α (−2ir 2 y)∂ y W * β,α (−2iy), and as before, we have introduced y ≡ −k 1 τ . Figure 4 shows the gravitational wave bispectrum from the diagram (ii). The magnitude of this bispectrum is smaller than that of the diagram (i) by a factor of a few for this choice of parameters (however, see section 4.4 for more details on the relative amplitudes and shape of the bispectrum), and they have the opposite signs. The bispectrum from the diagram (ii) peaks in the equilateral configuration, r 2 = r 3 = 1.

Diagram (iii)
For the diagram (iii), we consider the second order gravitational waveψ R 2 sourced by a first order gauge field perturbationt R 1 and a first order metric perturbationψ R to be, where the source term is again written as the sum of two parts, 1ij (τ, q 1 , q 2 ) + S The first part depends on the fields without time derivatives, whereas the second part includes time derivatives of the gauge field perturbation, However, we find that the first part multiplied by the polarisation factor vanishes: e L ij (k 3 )S (iii) 1ij (τ, −k 1 , −k 2 ) = 0, for all the permutations of k 1 , k 2 , and k 3 . Hence we consider only the second part. With Green's function, we find the second order gravitational wave aŝ The first term in eq. (4.2) from the diagram (iii) yields The expectation value is calculated as Then, equation (4.43) reads Combining with the other two terms in eq. (4.2), we obtain the bispectrum from the diagram (iii) as where we defineW 1 (y) = ∂ y W β,α (−2ir 2 y)Φ(−2ir 3 y) + ∂ y W β,α (−2ir 3 y)Φ(−2ir 2 y), (4.48) Figure 5 shows the momentum dependence of the bispectrum from the diagram (iii). We find that the contribution from the diagram (iii) is almost 7 orders of magnitude smaller than that of the diagram (i) and (ii), justifying that we neglected its contribution in our previous work [30]. This diagram is also zero in the folded limit and the bispectrum peaks in the equilateral configuration. The contribution from this diagram is so small that we do not compare the templates to it.

Total bispectrum
Combining the three contributions which we have calculated in the previous subsections, eqs. (4.23), (4.36) and (4.46), we obtain the total bispectrum of the sourced GWs in our model as The B dependence of Υ is weak, since the third line in eq.(4.50) from the diagram (iii) is negligible compared to the others. In figure 6, we show the shape of the total bispectrum (k 1 k 2 k 3 ) 2 B RRR h (k 1 , k 2 , k 3 ). We find a mild peak around r 2 r 3 ≈ 0.6 as a result of the combination of the plateau of B  figure 15. To quantitatively measure similarity, we calculate the "cosine" between the shape of our bispectrum and the equilateral shape, cos(B h · F eq ) [40]. Definition of the cosine is described in appendix B. In figure 7, the cosine is shown as a function of m Q . Note that the cosine depends only on m Q , because H and B change only the overall amplitude as long as the diagram (iii) is negligible. We find that the cosine rises from about 0.5 for m Q ∼ 2 to around 0.9 for m Q 2.5 and it varies up to 1% for the parameter range of interest (see section 6). This is because for values of m Q < 2.5 the total bispectrum receives significant negative contribution from the second diagram close to the equilateral limit, thus suppressing the total bispectrum relative to the peak in this region. As a result, the shape becomes quite different from the equilateral shape.
Around 90% similarity to the equilateral shape implies that our gravitational wave bispectrum is reasonably characterized by the amplitude at the equilateral limit, r 2 = r 3 = 1.
In the equilateral limit, the factors in the bispectrum become Ξ(r 2 = r 3 = 1) = 27 64 (4.51) and Υ eq (m Q ) ≡ Υ(m Q , 1, 1) where Re[z] denotes a real part of a complex number z, and the small contribution from the diagram (iii) is ignored. In figure 8, we plot Υ eq . For 3 m Q 5, Υ eq is well approximated by the following expression: Here, the relative error of this fitting formula is less than 1%. The ratio of the bispectrum to the squared power spectrum of GWs from the vacuum fluctuation of the metric, B vac h /(P vac h ) 2 , is of order unity [31,32]. The ratio for the sourced GWs can be much greater than unity. From eqs. (3.18) and (4.49), the ratio in the equilateral limit is given by : Υ eq (m Q ) defined in eq. (4.52) and its fitting formula (eq. (4.53)) are plotted as the red dots and blue line, respectively.
In figure 9, we plot B B sourced h /(P sourced h ) 2 in which steep exponential dependence of the bispectrum and power spectrum on m Q cancels out, though milder exponential dependence remains. We find a simple relation for a specific range of m Q ; for instance, where Ω A ≡ ( B + E )/2 (1 + m −2 Q ) B /2 (see eq. (2.10)) is the energy density fraction of the background SU (2) gauge field. Note also, that there is a "kink" in the bispectrum at m Q ∼ 2.25. This corresponds to the value of m Q for which the bispectrum from the second diagram is larger in magnitude than diagram one at the equilateral configuration, because of which the total bispectrum becomes negative. Since we plot the absolute value of the bispectrum, this appears as a "kink" in figure 9.
Dependence on the energy density fraction of the gauge fields in eq. (4.55) is analogous to the curvaton mechanism [36,37], where a similar relation holds for the scalar non-Gaussianity parameter, f NL ∼ Ω −1 σ . Ω σ is the energy density fraction of the curvaton field at its decay time. Therefore the origin of the dependence in eq. (4.55) may be understood in a similar way as the curvaton case [41]: suppose that the metric perturbation h is given by h = c 1 t where t is the mode function of the gauge field. At the same time, t is expanded as t = t (1) + t (2) + O(t (3) ) such that t (2) = c 2 (t (1) ) 2 . Then B h /P 2 h ∼ c 2 /c 1 . From equations (3.2) and (3.5), we see that c 1 ∝ √ B H/M P . We also see from equation (4.9) that c 2 ∝ g. Thus B h /P 2 h ∼ gM P /( It should be noted, however, that this relation only holds when the gauge field has the dominant contribution to both the tensor power spectrum and bispectrum, and thus, is not valid in the limit B → 0. Numerical Analytic Fit Figure 9: Ratio of the absolute value of the sourced tensor bispectrum to the sourced tensor power spectrum squared (normalised by B ) as a function of the parameter m Q , along with its fitting formula (eq. (4.55)), plotted as the red dots and blue line respectively. It rises exponentially as m Q increases. Since the y-axis shows the ratio normalised by the energy density fraction of the gauge field, B , the ratio of the correlation functions can be very large for allowed values of B .

Peak of Bispectrum
The total tensor bispectrum given in eq. (4.49) has a peak not at the equilateral limit r 2 = r 3 = 1 but at r 2 = r 3 ≈ 0.6 (see figure. 6). In this section, we study why this happens, by looking into the evolution of the tensor perturbations of the SU (2) gauge field t ij . The shape of the tensor bispectrum is determined by N i (m Q , r 2 , r 3 ) andÑ i (m Q , r 2 , r 3 ) (i = 1, 2, 3) as well as Ξ(r 2 , r 3 ) in eq. (4.49), while the contributions fromN i (m Q , r 2 , r 3 ) are negligible. For our current purpose, it suffices to focus on the case with r 2 = r 3 in which the three momenta, k 1 , k 2 , k 3 , form an obtuse-angled isosceles triangle. In other words, we concentrate on a cross-section surface of the 3D plot, figure 6. For r ≡ r 2 = r 3 , the r dependence of N i (m Q , r) andÑ i (m Q , r) is shown in figure 10. We do not plot |N 3 | and |Ñ 3 | which are the same as |N 2 | and |Ñ 2 | for r 2 = r 3 , respectively. We find that only N 1 grows significantly as r decreases, while the others have moderate dependences. To understand its behaviour, we look closely at the second line of N 1 in eq. (4.24) for r 2 = r 3 , In the integrand of I 1 , the first part, Im[W * β,α (−2iy)W β,α (−2iz)], represents Green's function for t R 2 given in eq. (4.14), and the second part, z(1 + 2r) − 5m Q − 2m −1 Q W 1 (z), represents the non-linear source term from the first order tensor perturbations, t R 1 × t R 1 . They are plotted in figure 10 for r = 1 and r = 0.5. Basically the second part is shifted by a factor of ≈ 2 along the z-axis, as r is reduced to the half. However, Green's function has a bigger amplitude at larger z without oscillations up to z 10. Note that the non-linear source term contains W 1 (z) ≡ W β,α (−2ir 2 z)W β,α (−2ir 3 z), indicating that the two sourcing modes t R 1 have momenta k 2 = k 3 = rk 1 in the case of r ≡ r 2 = r 3 , while the momentum of the sourced mode t R 2 is k 1 in the process of N 1 . The fact that Green's function for t R 2 is larger on sub-horizon scales implies that the source effect from t R 1 × t R 1 to t R 2 is more efficient when the sourcing modes t R 1 (rk 1 ) have lower momenta (i.e. a smaller r) and get amplified before t R 2 (k 1 ) crosses the horizon. In other words, N 1 becomes larger for a smaller r, because atypical Green's function G t allows the sourcing effect to be active deep inside the horizon.
This non-linear sourcing process of t 2 through the diagram (i) shows clear contrast from the linear sourcing process from t 1 to ψ 1 discussed in section 3. There, Green's function for ψ, G ψ , rapidly oscillates inside the horizon and does not allow t 1 to induce ψ 1 on sub-horizon scales, as shown in figure 1. In cases where only such normal Green's functions are involved, the shape of the bispectrum is typically equilateral, since all the modes are mainly produced around the horizon crossing. Nonetheless, in our case, Green's function for t R is peculiar due to tachyonic instability, and the peak of the bispectrum deviates from the equilateral limit.
The total contribution to the bispectrum from N i andÑ i is maximal in the folded limit r = 0.5. However, Ξ(r 2 , r 3 ) arising from the tensor polarisations is also an important factor determining the shape of the tensor bispectrum. Ξ is multiplied to the total bispectrum eq. (4.49) as an overall factor and it vanishes at r 2 = r 3 = 0.5. In figure 11, we illustrate how Ξ changes the shape of the bispectrum on the r 2 = r 3 plane. Ξ suppresses the bispectrum at lower r and vanishes at r = 0.5. In fact, Ξ vanishes not only at r 2 = r 3 = 0.5, but at all points on the line r 2 + r 3 = 1 (i.e. the folded limit), because of conservation of angular momentum. The Feynman diagrams in figure 2 can be seen as processes in which two spin-2 particles collide and one spin-2 particle comes out. In particular, in the case of a head-on collision, which corresponds to the folded limit, k 2 + k 3 = k 1 , the cross-section vanishes, because the angular momentum is contributed only by spins (i.e. no orbital angular momentum) and the h ) is plotted as solid green line as an function of r ≡ r 2 = r 3 . The peak is located at r ≈ 0.6 and the deviation from the equilateral shape (red dotted) is remarkable. The blue dashed line shows the case without Ξ, namely 10 −1 ( h )/Ξ which is multiplied by 0.1 for illustrative purpose. The parameters are given in eq. (4.8). spin of the system cannot be conserved as 2 ± 2 = ±2.
In summary, the peak of the tensor bispectrum B h is located at r ≡ r 2 = r 3 ≈ 0.6 for the following two reasons -(i) Among several contributions to the sourced tensor bispectrum B RRR h , the biggest one comes from ψ R 2 (k 1 )ψ R 1 (k 2 )ψ R 1 (k 3 ) ∝ N 1 in which the two linear perturbations of SU (2) gauge fieldt R 1 (k = rk 1 ) non-linearly induce the second order onê t R 2 (k 1 ) and subsequentlyt R 2 (k 1 ) sources the second order GWψ R 2 (k 1 ). In this process, the amplitude of the second order fluctuations is larger when the momentum of the first order perturbations, rk 1 , is smaller, because in this case t R 1 gets amplified when the second order t R 2 is still deep inside the horizon where the source effect is more efficient (i.e. Green's function G t has a bigger amplitude). Hence, N 1 , which dominates the tensor bispectrum, is a decreasing function of r (see figure 10). (ii) The polarisation tensors also yield a r dependence as an overall factor Ξ(r) to the total bispectrum. Ξ(r) is a growing function of r and vanishes at r = 0.5. Multiplying Ξ(r) changes the blue dashed line into the green line in figure 11. As the result of (i) and (ii), we obtain the bispectrum with a peak at r 0.6.

Parameter Search
In this section, we constrain the parameter regions from present observations and selfconsistency of the model. We also clarify the parameter regions where the power spectrum or bispectrum of the sourced GWs will be detectable by upcoming CMB observations. Note that there remain four parameters, H, m Q , B and g, and one relationship g 2 B M 2 P = m 4 Q H 2 . Eliminating g, we are left with three free parameters, H, m Q and B , in our model.

Tensor-to-scalar ratio
Currently the CMB observations put an upper bound on the tensor-to-scalar ratio r as where P ζ is the power spectrum of the curvature perturbation and k CMB = 0.05 Mpc −1 . In our model, not only the vacuum fluctuation of the metric but also the sourced GWs contribute to P h . Substituting eq. (3.18), the total tensor-to-scalar ratio is given as where the dimensionless scalar power spectrum, ∆ ζ ≡ k 3 P ζ /2π 2 ≈ 2.2 × 10 −9 [42], and that of the tensor metric vacuum fluctuation, ∆ h vac ≡ k 3 P vac h /2π 2 = 2H 2 /π 2 M 2 P , are introduced. Translating the upper bound on r into the constraint on our model parameters, we obtain where the conventional tensor-to-scalar ratio contributed only from the tensor metric vacuum fluctuation is defined by Since the upcoming CMB B-mode polarisation observation missions aim to achieve a sensitivity r ≈ 10 3 , the parameter region predicting r ≥ 10 −3 is particularly interesting. In our model, we find (6.5)
From the constraint on f tens NL we find, where the first constraint applies when Υ eq < 0 and the second when Υ eq > 0. 4 The factor of 2 √ 2 in the denominator comes from the difference of the normalisation of the polarisation tensors. In [42], e R ij (k)e R ij (−k) = 2 is adopted.

Consistency of the model
In addition to these observational constraints, we discuss the restriction imposed by selfconsistency of the model. Scalar perturbations of the spectator sector have a fatal instability on sub-horizon scale if m Q < √ 2 [35]. Hence we demand m Q > √ 2 in our model. Since B approximately indicates the energy density fraction of the background SU(2) gauge field, B is positive and small. As found in [43], if B is too large, its effect on the evolution of the inflaton perturbation significantly alters the spectral index n s , because B contributes toḢ through eq. (2.11). To keep this effect negligible, it is required where t * is the time at which CMB modes leave the horizon. On the other hand, since B can be rewritten as B = m 4 Q H 2 /(g 2 M 2 P ), if one lowers B by fixing m Q and H, one would confront a large self-coupling constant g of the SU(2) gauge fields which leads to a non-negligible backreaction from SU (2) tensor perturbations to the background dynamics. In order to avoid large backreaction, we need to have [20] where B andB are functions of m Q given by We also find that equation (6.11) ensures g 1 as well, which is preferred for validity of the perturbation series. 5

Allowed parameter regions
Equations (6.1)-(6.11), together with m Q > √ 2, give the set of constraints we employ to define the regions in B -m Q plane that are interesting for future CMB experiments. Figures (12)- (14) show the allowed regions for 3 different choices of H or equivalently r vac = 10 −4 , 10 −3 and 10 −2 . As r vac increases, the allowed parameter space shrinks. This is because the upper bound on r implies that the power in the sourced tensor modes cannot be very large if r vac is large.
The bottom right corners in these figures (i.e. regions with a large m Q and small B ) are shaded as the parameter spaces with non-negligible backreaction, although this does not mean that these regions are excluded. Rather, it indicates that one needs to perform numerical calculations to take into account backreaction, to study this parameter space [20] (see ref. [43] where the backreaction is numerically incorporated).
We find that there's a general trend in the constraining power of tensor power spectrum and bispectrum. While the power spectrum is better at constraining small m Q regions, the bispectrum is better at constraining large m Q regions. This happens because B h is exponentially more sensitive to m Q (it has an extra factor of e 2π(2m Q +m −1 Q ) compared to r), and so, a small change in m Q can easily change B h by a large factor (∼ e 4π∆m Q ). This also has interesting consequences for detectability of the tensor bispectrum, as a large range of bispectra can be generated even for the small range of values of interest, √ 2 < m Q 4, making a detection of the tensor bispectrum (in this model) possible in the near future, even if r is small (see figures (12)-(14)). While even the current constraints on f tens NL are useful for ruling out the top right corners in the figures (i.e. regions with a large m Q and large B ), there remains parameter space in which the tensor bispectrum can be observed in the future. It is then natural to ask what range of parameters can be probed in upcoming CMB missions.
To that end we also plot the line for σ(f tens NL ) = 1 in figures (12)- (14), which is expected to be the target sensitivity of LiteBIRD (M. Shiraishi, private communication). We see that this improved sensitivity will allow us to probe a significant portion of the parameter space with large m Q and small B , which is inaccessible to measurements of r, even if we can measure r = 10 −4 (figure 12, bottom right). Although our present calculation does not ensure that this conclusion stays unchanged when we account for the backreaction, it might still be true when backreaction is included. We also show the line corresponding to r source = r vac . Regions to the left of this line denote regions where the amplitude of the sourced tensor modes is smaller than the amplitude of vacuum tensor of the metric. From figure 14 we see that if r vac = 10 −2 , there is a region of intermediate B and m Q values for which f tens NL > 1, also if r source < r vac . This regime is particularly interesting because one can learn about both vacuum fluctuations of the metric and spectator fields during inflation, by combining the power spectrum and bispectrum. On the other hand, if r vac is smaller, a small r source is accompanied by a small tensor bispectrum as well.

Conclusion
In this paper we have calculated the bispectrum of tensor perturbations sourced by spectator SU(2) gauge fields during inflation [20]. The primary contribution to the bispectrum comes from the self-interaction of the SU(2) gauge fields; thus, it is unique to non-Abelian gauge theory. We find that the amplitude of the bispectrum parametrised by its ratio to the (squared) power spectrum, B h /P 2 h , is very large, ∼ 1/ B [30]. Since B 1, this is much larger than ∼ 1 which is predicted for quantum fluctuations of the metric [31,32].
We also explored parameter space of the model relevant to future CMB missions. Even with an r vac as low as 10 −4 , large parameter space remains consistent theoretically as well as with the current CMB observations. However, the exponential sensitivity of the power spectrum and bispectrum on model parameters makes it difficult to completely eliminate all the parameter space of the model on the basis of just these observations.
Upcoming CMB missions such as LiteBIRD [44] will measure the CMB polarisation to unprecedented accuracy. This will allow us to not only detect B-modes but also to characterise them, hence testing one of our most ambitious claims about our origins. If the primordial Bmodes arise from quantum fluctuations of the metric, we will find them to be parity invariant, near scale-invariant, and weakly non-Gaussian. If not, we will use the deviations to constrain the fraction of energy density in spectator gauge fields in the inflationary Universe [25,30,45,46]. , respectively. The light red shaded region is not necessarily ruled out but a significant backreaction requires a dedicated numerical treatment to obtain the predictions. In the orange shaded region, the system confronts a strong coupling problem, if one considers SU (2) charged particle. We also show f tens NL = 1 as the dashed blue line, because an error of order unity σ(f tens NL ) ∼ 1 would be achieved by upcoming CMB B-mode missions. The solid lines denote r = 10 −2 (blue), 10 −3 (green) and 10 −4 (yellow). Although the general expression for e L/R ij (k) is rather complicated, we can fix θ in the current case. This is because we calculate three polarisation tensors with three different wavenumbers, e R ij (k 1 )e R kl (k 2 )e R nm (k 3 ) whose indices are somehow contracted, and these wave vectors are on the same plane due to momentum conservation, δ(k 1 + k 2 + k 3 ). In that case, we can set θ = π/2 and let these vectors, k 1 , k 2 , k 3 , move only on the x-y plane. For θ = π/2, the polarisation tensors become Now we have three angles, ϕ 1 , ϕ 2 , ϕ 3 , associated with wavenumbers, k 1 , k 2 , k 3 , respectively. Without loss of generality, we can set ϕ 1 = 0. Furthermore, these trigonometric functions of ϕ 2 and ϕ 3 can be rewritten as functions of r 2 ≡ k 2 /k 1 and r 3 ≡ k 3 /k 1 . Using k 1 +k 2 +k 3 = 0, we find k 2 3 = |k 1 + k 2 | 2 = k 2 1 + k 2 2 + 2k 1 k 2 cos ϕ 2 =⇒ cos ϕ 2 = r 2 3 − r 2 2 − 1 2r 2 . (A.8) In the same way, we also find cos ϕ 3 = (r 2 2 − r 2 3 − 1)/(2r 3 ). With this notation, we find

B Equilateral Shape
To measure similarity of the shapes of bispectra, the cosine between two shapes is introduced as [40], where the dot product is defined as X · Y ≡ Here F ref is the reference template to which the similarity is measured. In this paper we use the equilateral template [47] F eq (k 1 , k 2 , k 3 ) = −  Figure 15 shows the shape of this template as a function of r 2 and r 3 .