Backreaction of axion-SU(2) dynamics during inflation

We consider the effects of backreaction on axion-SU(2) dynamics during inflation. We use the linear evolution equations for the gauge field modes and compute their backreaction on the background quantities numerically using the Hartree approximation. We show that the spectator chromo-natural inflation attractor is unstable when back-reaction becomes important. Working within the constraints of the linear mode equations, we find a new dynamical attractor solution for the axion field and the vacuum expectation value of the gauge field, where the latter has an opposite sign with respect to the chromo-natural inflation solution. Our findings are of particular interest to the phenomenology of axion-SU(2) inflation, as they demonstrate the instability of the usual trajectory due to large backreaction effects. The viable parameter space of the model becomes significantly altered, provided future non-Abelian lattice simulations confirm the existence of the new dynamical attractor. In addition, the backreaction effects lead to characteristic oscillatory features in the primordial gravitational wave background that are potentially detectable with upcoming gravitational wave detectors.


Introduction
Natural inflation was originally proposed in order to address the UV sensitivity problem of inflation [1,2].By identifying the inflaton with an axion, the pseudo-scalar field possesses a shift symmetry that protects the inflationary potential against radiative corrections.The inflationary potential is generated through instanton effects and in the early versions of natural inflation the potential is a sinusoidal function of the field.
The shift symmetry also dictates the possible couplings of the inflaton to other fields, since an axion can only couple derivatively to other fields.In particular, the lowest order coupling to fermions ψ and gauge bosons A µ are f −1 (∂ µ ϕ) ψγ µ γ 5 ψ and f −1 ϕF µν F µν , respectively, where F µν is the field strength tensor and Fµν = ϵ µναβ F αβ is the dual field strength tensor (see e.g.[3][4][5][6]).
The coupling of an axion inflaton to gauge fields through a Chern-Simons term has attracted significant attention in the literature.In the case of an Abelian field, the parityviolating nature of the coupling leads to the two helicities developing a different effective frequency.One of them can even become tachyonic, when the velocity of the inflaton is high enough.After the end of inflation, tachyonic production of gauge fields can lead to instantaneous preheating [6].Identifying the gauge field with the hypercharge sector of the Standard Model can lead to the generation of observationally relevant cosmological magnetic fields [7,8].During inflation, the production of gauge fields can lead to observable non-Gaussianity.
Depending on the axion-gauge coupling strength, the tachyonic amplification of the gauge fields can arise during inflation.In this case, the generation of gauge fields leads to a significant backreaction on the inflaton, leading to a sudden drop of its velocity.Once the gauge fields are diluted by the expansion of space-time, the backreaction term subsides and the inflaton starts rolling again.This can lead to periodic bursts of gauge field production during inflation [9][10][11][12][13][14][15][16][17], as has been shown analytically and numerically.However, recent lattice simulations showed that the inclusion of inhomogeneous backreaction and a larger dynamical range can significantly change the resulting dynamics [18].
Despite the interesting backreaction dynamics that occurs for large axion-gauge coupling, the rolling of the axion in the linear regime is determined by the potential and Hubble friction.By replacing the Abelian field with a non-Abelian one, this ceases to be true.The fact that SU(2) fields (see ref. [19] for the generalization to SU(N) fields) possess a non-trivial vacuum structure leads to a new inflationary attractor, in which the dominant source of friction for the rolling axion is not the Hubble term, but the non-Abelian field VEV [20][21][22][23][24].
This family of models, collectively named chromo-natural inflation, allows for slow-roll inflation even in steep potentials.Due to the parity-violating Chern-Simons coupling, one of the tensor modes of the SU(2) sector experiences a similar instability to one of the helicities in the Abelian case.However, the fact that the SU(2) tensor mode is linearly coupled to the gravitational sector leads to a similar enhancement of chiral gravitational waves (GWs).
While chromo-natural inflation with a cosine potential has been shown to be incompatible with CMB observations [24], spontaneous breaking of the SU(2) symmetry or going beyond the cosine potential can bring the model in agreement with current data [25,26].A further generalization of CNI was proposed in ref. [27], where the axion-SU(2) action was treated as a spectator sector.This separates the inflationary sector, which is responsible for scalar fluctuations, from the chromo-natural sector, which can produce detectable B-modes, while remaining subdominant in both scalar fluctuations and energy density during inflation.This family of models can be described as "spectator chromo-natural inflation" (SCNI) and their GW spectra are directly related to the shape of the axion potential [28].Currently, the dynamics of spectator non-Abelian gauge fields during inflation is an active research direction [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].
Previously, the effects of backreaction in spectator axion-SU(2) inflation were estimated by comparing the backreaction contributions to other terms in equations of motion [27,31,48,49], or more generally in ref. [50] for the case when the backreaction integral is regularized 1 .The authors of the ref.[50] studied the stability of slow-roll dynamics during axion-SU(2) inflation and found that the slow-roll solutions become unstable when the backreaction terms dominate in the equations of motion.In this work, we go beyond previous studies of axion-SU(2) dynamics during inflation by considering the effects of backreaction without the regularization of backreaction integrals.We use the linear equations of motion for the tensor SU(2) fluctuations and self-consistently solve the background equations for the axion and SU(2) VEV, including the homogeneous (averaged) backreaction from tensor fluctuations.This is in spirit similar to the analysis performed in ref. [16] for the Abelian case.We must note that ref. [15] largely validated these calculations, while more recent simulations [18] point out the importance of inhomogeneous effects during the strong backreaction regime.Our analysis can therefore be considered an important and necessary next step into the uncharted backreaction regime of axion-SU(2) dynamics during inflation.Our use of the linearized equations of motion for the gauge fields does not allow us to capture possible non-Abelian effects.While we can safely conclude that backreaction leads to the destabilization of the "standard" chromo-natural attractor, full lattice simulations are required to validate the existence of the new dynamical late-time attractor described below.This is left for future work.
Our manuscript is organized as follows.In section 2 we review spectator chromo-natural inflation and provide the necessary equations and analytical solutions.The numerical procedure is described in section 3, followed by the results and semi-analytical analysis of the solution.We conclude in section 4. In Appendix A we present the set of parameters used for the numerical computations.In Appendix B we discuss possible numerical artifacts.In Appendix C we give details of the backreaction dynamics and in Appendix D discuss the behavior of tensor perturbations on the dynamical attractor solution found in this work and provide a comparison with axion-U(1) inflation.

Review of spectator axion-SU(2) inflation
In this section, we review the spectator axion-SU(2) inflation or the spectator chromo-natural inflation model, outlining the background and perturbation analysis based on previous works [22,24,27].

Model and background evolution
The action for spectator axion-SU(2) inflation is given by [27] ) where R is the space-time Ricci scalar, ϕ(t) and V (ϕ) are the inflaton field and its potential, respectively, χ(t) and U (χ) are the axion field and its potential, ν is the field strength of the SU(2) gauge field A a µ , F a µν = ϵ µνρσ F a ρσ / 2 −det g µν is its dual (ϵ µναβ is the antisymmetric tensor and ϵ 0123 = 1), g is the gauge field coupling, λ is the coupling constant between the gauge and axion sectors, f is the axion decay constant, and M pl is the reduced Planck mass.
In this work we use the axion potential of the form where µ is a constant that sets the energy scale of the axion field.In this convention, the axion field takes values χ ∈ [0, πf ].The potential for the inflation field, V (ϕ), is left unspecified.We work with the FLRW metric where i, j indicate the spatial directions.An isotropic solution for the background is given by the gauge field configuration [20,21] which is also an attractor [51].For this ansatz, the closed system of equations for the vacuum expectation value (VEV) of the gauge field Q(t) and the Hubble parameter H(t) is given by where V ϕ (ϕ) = ∂V (ϕ)/∂ϕ, V χ (χ) = ∂U (χ)/∂χ, and an overdot denotes a derivative with respect to cosmic time t.The Hubble slow-roll parameters are defined as which are much smaller than unity during inflation.The first slow-roll parameter ϵ H contains contributions from the inflaton field ϕ and the spectator sector that consists of the axion and gauge fields The various contributions are defined as (2.12) For the axion-gauge sector to remain a spectator, their energy densities must be subdominant to that of the inflaton, i.e., where the energy densities are given by The original chromo-natural inflation model in the slow-roll approximation has an attractor solution [22,24] The VEV of the gauge field that minimizes the axion effective potential is which is a solution of equation (2.15) when Q is constant.It is convenient to introduce the parameters where the dimensionless mass parameter m Q characterizes the mass of the gauge field fluctuations and controls their amplification.On the chromo-natural inflation attractor, m Q and ξ are related via ξ ≃ m Q + 1/m Q .

Perturbations
Let us now review the perturbations in the spectator axion-SU(2) model.We adopt the gauge choice and decomposition for field fluctuations following ref.[52] of the form together with (2.19) The perturbations consist of seven scalar modes (δϕ, δχ, Y, δQ, M, φ, B), six vector modes ) and four tensor modes (t + , t × , h + , h × ).At the linear level, all perturbations are decoupled from each other.Vector perturbations decay on superhorizon scales and at the linear level metric fluctuations can be neglected [52].The scalar perturbations have been studied in detail at the linear [27] and nonlinear levels [49,53].It was shown that for m Q < √ 2, the scalar perturbations are tachyonically unstable [52].The combination from linear and nonlinear analyses leads to the constraint [49] √ on the parameter m Q , where P ζ,CMB = 2.1 • 10 −9 .
Let us now turn to the discussion of tensor perturbations.It is convenient to express the plus and cross polarizations of tensor perturbations via the left-handed and right-handed polarizations as We canonically normalize them by introducing We will work with conformal time defined as which, with the near de Sitter expansion, leads to for the scale factor.In the following, derivatives with respect to η are denoted by primes.
In conformal time and to leading order in slow-roll, the equations of motion for the tensor perturbations are Here k is the wave number.The spectator axion-SU(2) model is known to have a transient growth of one of the polarizations of the gauge field tensor modes that leads to the production of chiral GWs.The produced GW background is enhanced with respect to the standard single-field slow-roll models of inflation with predictions potentially observable by near-future B-mode experiments.The total GW power spectrum is defined as where P tot h (k) can be expressed in terms of left and right polarization modes as where P (s) h is the late-time sourced GW power spectrum, defined as (2.29) It is convenient to introduce a parameter that characterizes the enhancement of the GW background with respect to the vacuum prediction, where the vacuum prediction for the power spectrum is given by In order to have sizable GW production, the parameter range of the model has to be chosen such that R GW ≳ 1.This requirement leads to the constraint [49] g ≲ 1.8 The growth of tensor modes results in the backreaction [27,31,[48][49][50] on the background equations of motion (2.5)-(2.9).Taking into account the contribution from backreaction, the background equations of motion in conformal time take the form with H = a ′ /a.The backreaction terms are evaluated in the Hartree approximation, leading to the integrals over mode functions (2.36) It is worth noting that in this work we consider homogeneous backreaction, where the spatial gradients of inflation and axion fields are neglected, keeping these fields homogeneous during inflation.
In the small-backreaction regime with an approximately constant m Q parameter, the spectator axion-SU(2) model can be solved analytically [23].The regime of small backreaction is achieved with the constraint [49] (2.37) In figure 1, we plot the three constraints (2.20), (2.32), and (2.37)2 on the parameter range of the theory, indicating the fiducial parameters used in the current work.

Numerical treatment of the backreaction
To study the backreaction in axion-SU(2) inflation, we now solve the perturbation equations along with those for the background numerically.We begin by describing the numerical method and then discuss the results.

Numerical implementation
We solve equations (2.26) and (2.25) for both the left-and right-handed components of T R,L and ψ R,L .For each perturbation variable, we solve the equations for the real and imaginary parts and represent them on a k mesh to compute the integrals in (2.35) and (2.36).For most of our studies, we use the logarithmic wave number along with conformal time as the independent coordinates.In that case, we use n k points in ln k that are separated by uniform intervals in ln k in the range To solve the background (2.33) and (2.34), we compute the integrals (2.35) and (2.36) up to second-order accuracy.We advance the solution in conformal time using a third-order time-stepping scheme.The initial conformal time is η i and the final one is η f .In practice, we choose η i = −H −1 i with H i = a i H and a i = 1 along with η f = −10 −5 , which corresponds to a total duration of N ≈ 25 e-folds.The length of the conformal time step is then usually chosen to be ∆η = 10 −6 /H.
It is convenient to use the compute and data management infrastructure provided by the Pencil Code [54], which allows for efficient parallelization using the Message Passing Interface library.In some exploratory cases, we also solved the equations on a mesh where n min and n max grow in time such that the main contributions to the integral are captured during the entire evolution.
We initialize the perturbation variables with the Bunch-Davies initial condition.Specifically, we set the initial conditions for the real and imaginary parts of the perturbation variables as follows: The same initial conditions are used for ψ R,L .To discard the contributions from quantum vacuum fluctuations of T R,L in the calculation of the integrals in (2.35) and (2.36), we use the criterion that for wave numbers where |T R,L | 2 < 1/2k, the value of |T R,L | 2 is replaced by zero.

New late-time attractor solution
We have performed a range of simulations with different values of µ, g, and λ; see Appendix A for a summary.In figure 2, we show the time evolution of Q and χ for runs µ1, µ3, µ4, and µ5, which corresponds to m Q = 2.44, 3.29, 3.58, and 4.19, respectively.When m Q is large enough, the backreaction of the perturbations becomes important and the system undergoes a transition to the new dynamical attractor with negative values of Q and a reduced velocity for χ after about N = 2-10 e-folds.
As a test of our numerical implementation, we choose the initial value of µ and g for one of the runs where the backreaction of the perturbations on the background evolution is negligible.Run µ1 with m Q = 2.44 is an example of such a case.It is evident from figure 2   green curve shows the analytic solution obtained using the homogeneous solution of T R in equation (2.26).This solution is given by [23,24], Here, W β,α (2ikη) is the Whittaker function with The Whittaker function provides a good analytical approximation for T R for a particular wave number approximately until the Hubble horizon crossing.However, the solution starts to differ in the deep superhorizon regime due to the contribution from metric tensor perturbations.It is evident from the upper panel of figure 3 that our numerical solution matches well with the analytical solution in the regime where it is valid.
In the lower panel of figure 3 we show the evolution of √ 2k(x|T R,L |) and √ 2k(x|ψ R,L |) for the case when backreaction is important.We see that the gauge field tensor perturbations as well as the metric tensor perturbations are getting amplified for larger values of m Q .However, as we have seen in figure 2, such amplification drastically changes the behavior of the background dynamics 3 .
It is worth discussing the evolution of the backreaction integrals T Q BR and T χ BR defined in equations (2.35) and (2.36).We show it in figure 4. From figure 4, we conclude that most of the contribution to the backreaction comes from a fixed narrow range of wave numbers.This range is different for different values of m Q .For run µ4 with m Q = 3.58, the range is around ln k/(a 0 H) ≈ 6; see figures 4(b) and (e), while for run µ5 with m Q = 4.15, it is ln k/(a 0 H) ≈ 4; see figures 4(c) and (f).For run µ1 with m Q = 2.44, there is no backreaction; see figures 4(a) and (d).We also performed variable k range integrations for the same runs by calculating the backreaction integrals for modes around the comoving horizon.However, such a comoving integration scheme leads to unphysical oscillatory features in the background evolution, as we show in Appendix B.

Semi-analytical modelling
In this section, we provide a semi-analytical analysis to approach the new dynamical attractor solution.In order to provide some intuition on the dynamics in the backreaction regime, it is instructive to investigate the evolution of each contribution to the equations of motion, (2.33) and (2.34).The time dependence of contributions is shown in figure 5. Following the evolution of the backreaction terms T Q BR , T χ BR , three distinct phases of dynamics can be distinguished.From the top left and bottom left panels of figure 5 one can see that the backreaction contributions (solid black curves) grow exponentially in absolute value up to around 8 e-folds.We refer to the stage of exponential growth of backreaction as Stage I.This leads to the destabilization of the "standard" chromo-natural attractor and pushes the system toward a new regime.When the backreaction contributions become comparable to one of the terms in the equations of motion, Stage II begins.At this stage, the backreaction terms change their behavior, start to decrease, and eventually cross zero.Following this stage, the system converges to the final solution (see the top right and bottom right panels of figure 5) which we refer to as Stage III 4 .The dynamics during different stages is described in more Let us turn right away to the discussion of the new late-time dynamical attractor.In Stage III, all the contributions to the equations of motion become nearly constant.In addition, the term −(gλ/f ) aχ ′ Q 2 becomes nearly equal to the contribution 2H 2 Q in the equation of motion (2.33) 5 , which leads to This solution resembles the original chromo-natural attractor solution given in equation (2.15) with Q = const and just with the second term present that has an opposite sign.It follows the late-time attractor ξ ≃ −1/m Q .In figure 6, we show that equation (3.4) does indeed hold.With (3.4), the equations of motion in Stage III become where we have taken into account that terms with derivatives in the last stage are negligible.
scope of the present work.We use the linear equations as an approximation to the full system. 5We used H ′ = a 2 Ḣ + H 2 and Ḣ = 0 in numerical simulations.Let us now take a closer look into the time dependence of each component of the backreaction integrals (2.35) and (2.36).Backreaction terms may be written in conformal time as where we have denoted the integrals over wave numbers by We have checked that in all our cases, the contributions from |T L | to the backreaction are negligible.In figure 7, we show the time evolution of the different contributions to the T Q

and T χ
BR integrals.There is a clear correspondence between the background dynamics and backreaction integrals, and vice versa.We refer the interested reader to the detailed discussion in Appendix C.
The backreaction integrals include crucial information that governs the evolution of the whole system.It is therefore convenient to introduce a new parameter that quantifies the ratio of the two backreaction integrals (3.10) From equations (3.5), one can find the relation between χ and Q on the final attractor At the final stage it holds T Q BR ≃ T Q 1 and T χ BR ≃ T χ 1 +T χ 3 , so the parameter α can be expressed as where we used ⟨|T | 2 ⟩ ′ ≃ 2aH⟨|T | 2 ⟩ and Q fin denotes the value of Q on the final attractor.The dependence of α on the parameters λ, g, and µ is confirmed in figure for runs of the three families.To sum up, the late-time dynamical attractor is given by equations (3.4), (3.11) and (3.12).

Observational signatures
The GW energy density power spectrum can be approximated as [55]  where Ω rad ≃ h −2 2.47 × 10 −5 is the present radiation density parameter and k eq ≃ 1.3 × 10 −2 Mpc −1 is the wave number entering the horizon at matter-radiation equality.Here, h is defined such that H 0 = 100 h km s −1 Mpc −1 is the Hubble parameter at the present epoch.To express Ω GW as a function of frequency, we use f ≃ 1.5 × 10 −15 (k/Mpc −1 ) Hz. Here, the frequency is not to be confused with the axion decay constant, which is also called f .
In figure 9, we show 2kx 2 |T | 2 vs k/k H and Ω GW vs f /f H for runs µ4, µ7, and µ4 * with m Q = 3.58, 5.15, 3.78 respectively (solid blue, dotted red, and dashed green curves).The run denoted as µ4 * is analogous to the µ4 run but with a Hubble parameter H larger by a factor of 10.In this figure, we have normalized the wave number and GW frequency using the values corresponding to the Hubble horizon size (k H ) at the end of inflation and the corresponding frequency (f H ).These quantities are given by Here, H = 1.04 × 10 −6 M pl and we assume an adiabatic evolution of the Universe to calculate a e /a 0 , given by a e a 0 = g 0s g rs In the above expression, g * s and g 0s denote the effective degrees of freedom in the entropy at the end of inflation and the present epoch, respectively, and T r denotes the reheating temperature assuming instantaneous reheating.We estimate T r by using the relation It is worth noting that we normalize the frequency in figure 9 by assuming that the end of the simulation is also the end of inflation.However, if the end of the simulation is N e e-folds before the end of inflation, the wave number and corresponding frequency must be multiplied by e −Ne .
As is evident from figure 9, the modes that are amplified around horizon crossing give the largest contribution to Ω GW .The oscillations observed in Ω GW are related to the oscillations in T R , as illustrated in the left panel of figure 9.The higher initial value of m Q in run µ7 results in a higher peak value of Ω GW .As demonstrated in figure 9, the higher value of H leads to a higher value of Ω GW .
In figure 10, we show h 2 Ω GW for runs µ4, µ7, and µ4 * (m Q = 3.58, 5.15, and 3.78, respectively), along with sensitivity curves for various GW detectors.The sensitivity curves are obtained by using the "strain noise power spectra" file available in the Zenodo online repository [56] associated with ref. [57].Our simulation covers approximately 25 e-folds.To provide a comparison of the obtained h 2 Ω GW with the sensitivity of upcoming GW detectors DECIGO and BBO, we rescale the GW frequency such that the end of our simulation occurs 14 e-folds before the end of inflation.For run µ4 * , the peak in Ω GW at f = 0.1 Hz is potentially detectable.

Summary and discussion
In this work, we simulated an axion-SU(2) sector, which is a spectator during inflation, meaning that its energy density is subdominant to the inflaton and that both the Hubble rate and the density perturbations are unaffected by its presence.In our simulations, the fluctuations are computed using the linearized equations of motion.Their effect on the axion and gauge field VEV are computed self-consistently, but averaged over the simulation domain.This method is similar in spirit to the one followed in ref. [16] in the case of axion-U (1) inflation, where the gauge fields are computed using linear equations of motion and their collective effect is considered a background quantity and added to the corresponding background equations.In the case of the Abelian model, an oscillatory result was found that can be understood in the following way: an increase in the axion velocity leads to an increase in the gauge field amplification.This results in an increased backreaction on the rolling axion through the ⟨E • B⟩ term; see Appendix D. This backreaction leads to a slow-down of the axion, which, in turn, reduces the subsequent amplification of gauge fields.Since the gauge fields red-shift after their production, their backreaction will also reduce, leading to a speedup of the axion and the whole process will start anew, thus leading to periodic bursts of gauge field production.
In the non-Abelian case, the initial stage is similar to the Abelian case: as the axion picks up speed, the gauge fields (a tensor mode in this case) are exponentially amplified.This leads to a backreaction on the equations of motion that define the VEV of both the axion field as well as the SU(2) sector.This leads to both a slow-down of the axion, as well as a completely new sign-flipped value for the gauge field VEV.Furthermore, in this new regime, the superhorizon tensor modes of the gauge field evolve as T R,L ∝ 1/η.This leads to the absence of the periodic behavior found in the Abelian case, because the backreaction of the gauge field fluctuations onto the background quantities does not diminish with time (to lowest order in slow-roll).
Having revealed this new regime, several questions remain to be answered.An intriguing relation was found between the backreaction terms in the axion and gauge VEV, leading to a universal relation between the parameters of the potential and the late-time value of the gauge VEV.However, we are not able to predict the gauge VEV itself in this new attractor.We believe that the initial value of the gauge VEV (in the original chromo-natural attractor) plays a role in determining its late-time value.
Furthermore, our analysis neglects spatially dependent backreaction effects that can lead to mode-mode coupling of the gauge field, as well as the excitation of scalar fluctuations in the axion sector.It has been shown in ref. [18] that space-dependent backreaction effects can be very important in the Abelian case.More importantly, large values of the gauge fields may lead to non-Abelian interactions becoming important and invalidate the linear approximation.While we expect the initial growth of the fluctuations and destabilization of the "standard" chromo-natural solution to hold, providing conclusive proof on the existence and stability of the late-time attractor requires solving the full system on a lattice, without making any linear or Hartree-type approximations.In addition, non-Abelian fields coupled to a rolling axion have been argued to lead to warm inflation [58,59].It is an intriguing possibility to examine if a system can transition from the chromo-natural inflation attractor to warm inflation, either in the original chromo-natural inflation formulation or in spectator models.Furthermore, our calculation was performed with a constant Hubble scale, in an exact de-Sitter background.While this can be an excellent approximation for several inflationary models, it does not allow us to probe the evolution of this new attractor close to the end of inflation, when Finally, the flipped sign of the gauge field VEV provides the possibility of amplifying the subdominant helicity of gauge tensor modes.Further analysis of this is left for future work, as it can provide interesting scale-dependent observables.For run µ4 * , we have used H = 1.04 × 10 −5 M pl , which leads to a higher value of Ω GW ; see section 3.4.

B Artifacts from not resolving the superhorizon modes
One might have expected that it is important to resolve the modes around the comoving horizon.Looking at figure 4, this is not obvious, however.Once backreaction becomes important, most of the contributions to the backreaction come from a fixed band of wave numbers.It is instructive to examine the results where we allow for the possibility to move the range of integration to a comoving strip with The result of numerical simulation using a comoving strip of wave numbers is shown in figure 11.
As we see from the insets of figure 11, the corresponding evolution of Q is different in cases where the modes in the proximity of the comoving horizon are resolved at the expense of not capturing any more the strongly superhorizon modes.We can conclude that this causes numerical artifacts that look like periodic bursts of gauge field production during inflation.

C Detailed description of backreaction stages
At the initial stage, referred to here as Stage I, the solution follows the chromo-natural attractor solution (2.15), where the three terms −(gλ/af )χ ′ Q 2 , 2g 2 Q 3 , and 2H 2 Q/a 2 balance each other in (2.33).Note that, in the present case of H = const, we have H ′ + H 2 = 2H 2 .When the backreaction T Q BR becomes important, the contribution −(gλ/af ) χ ′ Q 2 becomes more dominant compared to the rest of the terms (see purple dashed curve in the top left and bottom left of figure 5).To compensate for the increase of the sum T Q BR and −(gλ/af ) χ ′ Q 2 terms, the Q ′′ contribution becomes negative.It causes the change in the sign of Q ′ .This changes the behavior and turns on the Q ′ term, T χ 4 , in the integral T χ BR of equation (3.8), which produces a bump in ξ; see figure 12.This happens around N = 8 e-folds.The change in ξ makes the two terms T Q 1 and T Q 2 in (3.6) almost cancel each other, see figure 7.As a result, the T Q BR term becomes first negative and is then close to zero.This causes a decrease of Q.The steps at Stage I can be described by the following chain sequence: The three stages of dynamics in axion-SU(2) inflation with backreaction.Grid lines are the same as in figure 5. Right: Evolution of the ξ parameter, defined in equation (2.17) with respect to the number of e-folds, for the same run as in the left panel.
in ξ and causes the inflection of the T Q 1 contribution.As a result, the solution arrives at the final attractor (3.4) with At Stage III, we observe the following chain sequence: The three stages of evolution are summarized in Table 2.

D Dynamical attractor, gravitational waves and the Anber-Sorbo solution
It is interesting to compare the dynamical attractor solution found in this work to the Anber-Sorbo (AS) solution [61], which applies to axion-U(1) inflation.We start by briefly reviewing the AS solution.While it has recently been shown to be unstable [17], its analytical derivation can still be illuminating when compared to the dynamical attractor found in the present work.
We assume an exact de-Sitter expansion for simplicity, meaning that conformal time η is η = −1/Ha = −e −Ht /H.The equation of motion for the background axion field is where the gauge fields follow the linear equations (in conformal time) A + is approximated as For any finite range of comoving wave numbers, this contribution redshifts as a −4 .However, in the presence of a slow rolling axion field ϕ, new wave numbers are amplified all the time, meaning that k max ≃ 2ξaH, leading to ⟨E • B⟩ ∝ e 2πξ H 4 , which is constant in time.Simply put, the largest wave numbers are dominating the backreaction.
Let us now move to the case of the new dynamical attractor in the axion-SU (2) system.The linearized equation for the gauge field fluctuations at late times can be approximated as where we ignored terms proportional to T R,L that grow slower than η −2 at late times, as well as terms proportional to ψ R,L or ψ ′ R,L , which are suppressed by slow-roll quantities.During the new attractor, m Q ξ ≃ −1 (see equation (3.4)), so the above equation has a growing late-time solution T R,L ∝ 1/η, as seen for example in figure 3 (see also ref. [24] for a similar late-time evolution of T R,L in the original chromo-natural inflation model).
The backreaction terms T Q BR , T χ BR given in equations (2.35) and (2.36) can be calculated similarly to ⟨E • B⟩ to scale as T Q BR , T χ BR ∝ k max for late times.Thus, contrary to the AS solution, the new non-Abelian attractor does not require continuous amplification of new, ever-larger wave numbers.Instead, it is supported by a fixed range of comoving wave numbers.This is shown in figure 4.
Before we conclude, it is worth exploring the late-time superhorizon evolution of gravitational wave modes ψ R in the presence of the growing mode T R .By keeping only the dominant terms of eq.(2.25) for η → 0, we arrive at where we used T R = T 0 R /η at late times.This can be analytically solved to give

Figure 1 .
Figure1.Constraints on the spectator axion-SU(2) model similar as in ref.[49] with indication of parameters used in the current work.The stars correspond to the parameters (from left to right) m Q = 2.44, 3.29, 3.58, 4.19 with g = 0.011 (runs µ1, µ3, µ4, and µ5 from table1).
that Q remains constant in time for run µ1, as expected from the analytical results.Furthermore, we have compared our numerical results for √ 2k(x|T R,L |) and √ 2k(x|ψ R,L |) with the analytical solution and show the comparison in the upper panel of figure 3 for k = 10 −4 .Here, x = −kη is the dimensionless time variable.In this figure, the solid red and dotted blue curves show the numerical result for the right-and left-handed polarizations, respectively, and the dashed Figure 3. √ 2kx|T R,L | and √ 2kx|ψ R,L | vs N for k = 10 −4 for runs µ1 with m Q = 2.44 (upper panel), µ4 with m Q = 3.58 (lower panel).The solid red and dotted blue curves represent the numerical results for the right and left-handed polarization respectively.The dashed green curve shows the Whittaker solution for T R given by equation (3.3).

Figure 5 .
Figure 5. Top left: Contributions to the equation of motion for the VEV of the gauge field, Q, with respect to the number of e-folds at the initial stage and when backreaction occurs for the run µ3 with m Q = 3.29.The three vertical gray grid lines correspond to the moments when Q ′′ = 0 , Q ′ = 0, and then again Q ′′ = 0 (from left to right) respectively.Top right: Contributions to the equation of motion for Q during the transition to the final attractor solution for the same run.Vertical gray grid lines correspond to the moments when Q = 0, Q ′′ = 0, Q ′′ = 0 (from left to right) respectively.Bottom left: Contributions to the equation of motion for the axion field, χ, with respect to the number of e-folds at the initial stage and when backreaction is turned on for the same run as the top panels.Grid lines are the same as in the top left panel.Bottom right: Contributions to the equation of motion for χ during the transition to the final attractor solution for the same run.Grid lines are the same as in the top right panel.

Figure 6 .
Figure 6.Verification of equation (3.4), showing (λχ ′ /af ) fin versus −2H2 /gQ fin , where Q fin denotes the value of Q at the final attractor for multiple series of runs from table 1.

11 Figure 7 .
Figure 7. Left: Backreaction terms T Q BR and their contributions with respect to the number of e-folds for the run µ3 with m Q = 3.29.The dashed curve shows when T Q BR is negative.Right: Backreaction terms T χ BR with contributions with respect to the number of e-folds for the same run.Grid lines are the same as in figure 5. Dashed curves represent negative values and solid curves show when functions are positive.

Figure 8 .
Figure 8. Dependence of the parameter α on Q fin for three series of runs.
II, is characterized by a continuous decrease of Q.The backreaction T Q BR ≈ 0 remains small, and T χ BR ≈ const.The last stage, Stage III, begins when the gauge field VEV reaches zero, i.e., Q = 0.This changes the sign of terms with m Q , i.e., T χ 1 and T χ 3 from (3.8).This governs the change

Figure 12 .
Figure 12.Left: VEV of the gauge field Q versus the number of e-folds for run µ3 with m Q = 3.29.Grid lines are the same as in figure5.Right: Evolution of the ξ parameter, defined in equation (2.17) with respect to the number of e-folds, for the same run as in the left panel.

2 ) 2 Figure 13 .
Figure 13.The normalized mode-functions √ 2kx|T R | (left) and √ 2kx|ψ R | (right) as a function of e-folds N for the run µ4 (m Q =3.58) for different wave numbers.The solid red, blue, and green curves correspond the numerical results for k = 2.5 × 10 −4 , 3.7 × 10 −4 , and 5.0 × 10 −2 respectively, while the dashed black, brown, and magenta curves in the right panel represent the corresponding analytical results given in equation (D.8) in the superhorizon limit.

1 / 4 e πξ e − 2 √ 4 .
2ξk/aH , (D.3) where we take ξ > 0. The backreaction is ⟨E • B⟩ = − consider the amplified polarization and ∂A + /∂η ≃ √ 2ξkaHA + .By integrating from k = 0 to k max , we arrive at ⟨E • B⟩ ∝ e 2πξ k max a Evolution of Q and χ with N is shown for different initial parameters.The blue solid, green dot-dashed, black dotted, and red dashed curves correspond to runs µ1, µ3, µ4, and µ5 respectively from table 1 with g = 0.011 and m Q = 2.44, 3.29, 3.58, 4.19 respectively.For the blue curve, the backreaction is negligible therefore the value of Q remains constant.However, for a larger value of m Q (green, black and red curves), the backreaction effects are important.Using equation (3.4) from section 3.3, the velocity of χ is dχ/dN = −2Hf /(gλQ), which can be evaluated using the data of table 1, leading to dχ/dN = −8.2× 10 −6 , −7.7 × 10 −6 , and −6.9 × 10 −6 for the green dot-dashed, black dotted and red dashed curves respectively.These estimates agree with the χ evolution shown in this figure.
2/gQ fin , where Q fin denotes the value of Q at the final attractor for multiple series of runs from table 1.