Instability of cosmic Yang-Mills fields

There exists a small family of analytic SO(4)-invariant but time-dependent SU(2) Yang-Mills solutions in any conformally flat four-dimensional spacetime. These might play a role in early-universe cosmology for stabilizing the symmetric Higgs vacuum. We analyze the linear stability of these"cosmic gauge fields"against general gauge-field perturbations while keeping the metric frozen, by diagonalizing the (time-dependent) Yang-Mills fluctuation operator around them and applying Floquet theory to its eigenfrequencies and normal modes. Except for the exactly solvable SO(4) singlet perturbation, which is found to be marginally stable linearly but bounded nonlinearly, generic normal modes often grow exponentially due to resonance effects. Even at very high energies, all cosmic Yang-Mills backgrounds are rendered linearly unstable.


Introduction: a tale of three anharmonic oscillators
Classical Yang-Mills fields play a central role in various areas of theoretical physics, from QCD confinement to spin-orbit interactions in condensed matter theory and early-universe cosmology. Concerning the latter, scenarios have been proposed and analyzed for isotropic inflation driven by non-Abelian gauge fields, such as gauge-flation or chromo-natural inflation, by employing a homogeneous and isotropic Yang-Mills background in a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe (for a review, see [1]. More minimalistically, Friedan has recently put forward an evolution of the electroweak epoch based on the Standard Model and general relativity alone, where an oscillating isotropic SU (2) gauge field stabilizes the symmetric Higgs vacuum in a spatially closed FLRW spacetime [2]. It is therefore of natural interest to establish the stability features of such "cosmic Yang-Mills fields" against classical and quantum perturbations. While this has been done using cosmological perturbation theory in the context of gauge-flation, the issue remains unclear in the pre-inflation scenario of Friedan, despite some early partial analysis of Hosotani [3].
Here, we address this matter by performing a complete stability analysis of the only known family of analytic SU(2) Yang-Mills solutions in a closed FLRW universe. The situation differs from that of gaugeflation, since we do not break conformal invariance, and our homogeneous and isotropic gauge background would conformally map to an inhomogeneous Yang-Mills configuration in spatially flat FLRW spacetime. Our analysis, however, is partial in that we investigate only gauge-field fluctuations while keeping the metric fixed, for two reasons. Firstly, the background FLRW dynamics does not influence the gauge field, since the gauge sector is conformally invariant and we conformally map to a static metric. Secondly, the time scale of the gauge-field dynamics in the electroweak epoch is supposed to be hugely shorter than that of the spacetime geometry, thus we do not expect (non-conformal) metric fluctuations to exert a sizeable influence on the stability of the gauge-field configuration. Our results should therefore be relevant to Friedan's scenario. Nevertheless, a full cosmological perturbation theory of the combined Einstein-Yang-Mills system requires turning on also metric perturbations, something we reserve for a future task.
Even without coupling to gravity and independent of potential cosmological applications, the perturbation theory around any analytic classical field configuration is of general interest for assessing its relevance for quantum properties, since the determinant of the second variation of the action yields the leading quantum correction to a saddle point in the semiclassical analysis of the path integral. This aspect has been investigated for many classical Yang-Mills solutions but not yet for the ones studied here, to our knowledge.
It is generally impossible to find analytic solutions to the coupled Einstein-Yang-Mills system of equations, in part because they are coupled both ways. However, in a homogeneous and isotropic universe, where the metric is conformally flat, the Yang-Mills equations decouple due to their conformal invariance in four spacetime dimensions. Thus, if one can find isotropic Yang-Mills solutions on Minkowski, de Sitter, or anti de Sitter space, then their energy-momentum tensor will be compatible with any FLRW metric (of the same topology) and allow for an analytic computation of the scale factor from the Friedmann equation.
Fortunately, for the de Sitter case and a gauge group SU (2), such Yang-Mills configurations with finite energy and action are available [4,5,6]. They are most easily constructed on the cylinder (0, π) × S 3 , which is related to de Sitter space by a purely temporal reparametrization and Weyl rescaling [3,7,8,2], ds 2 dS4 = −dt 2 + 2 cosh 2 t dΩ 2 3 = 2 sin 2 τ −dτ 2 + dΩ 2 3 for t ∈ (−∞, +∞) ⇔ τ ∈ (0, π) , (1.1) where dΩ 2 3 is the round metric on S 3 , and is the de Sitter radius. This provides an explicit relation between co-moving time t and conformal time τ , and it fixes the cosmological constant to Λ = 3/ 2 . One employs the identification SU(2) S 3 to write an S 3 -symmetric ansatz for the SU(2) gauge potential A µ , which produces a solvable ODE for a parameter function ψ(τ ) of conformal time. This ODE has the form of Newton's equation for a mass point in a double-well potential yielding a first anharmonic oscillator. A two-parameter family of (in general time-dependent) solutions can be given in terms of Jacobi elliptic functions and describe SO(4) invariant Yang-Mills fields on the Lorentzian cylinder over S 3 .
These Yang-Mills solutions then exist (at least locally) in any conformally related spacetime, but the conformal transformation will ruin isotropy unless we restrict ourselves to spatially closed FLRW metrics, where we impose a big-bang initial condition a(0)=0, so that dτ = dt a(t) with τ (t=0) = 0 and τ (t=t max ) =: The lifetime t max of the universe can be infinite (big rip, a(t max )=∞) or finite (big crunch, a(t max )=0). Bouncing cosmologies as in (1.1) are also allowed but will not be pursued here. Since the energymomentum tensor of our Yang-Mills configurations is SO(4) symmetric, their gravitational backreaction will keep us inside the FLRW framework and merely modify the cosmic scale factor a(τ ). The latter is fully determined by the Friedmann equation in the presence of the Yang-Mills energy-momentum and a cosmological constant Λ, whose value may be dialed. It is well known that the Friedmann equation takes the form of another Newton equation. Its (cosmological) potential for the case at hand reads which is our second anharmonic oscillator (although inverted). Each pair (ψ, a) of solutions to the two systems (1.2) and (1.5) yields an exact classical Einstein-Yang-Mills configuration. One parameter in ψ is the conserved mechanical energy E in the potential V , which in turn determines the mechanical energyẼ for a in the potential W . This one-way coupling is the only relation between the two anharmonic oscillators.
For a physical embodiment of the cosmological constant, we may add a third player, for instance a complex scalar Higgs field φ in the fundamental SU(2) representation. The standard-model Higgs potential where v/ √ 2 is the Higgs vev and λv is the Higgs mass, gives us a third anharmonic oscillator. The dictate of SO(4) invariance, however, allows only the zero solution, φ ≡ 0, which provides us with a definite positive cosmological constant of where κ is the gravitational coupling. The full Einstein-Yang-Mills-Higgs action (in standard notation), where g is the gauge coupling, and the overdot denotes a derivative with respect to conformal time.
For large enough "gauge energy" E, the universe undergoes an eternal expansion, which is accompanied by rapid fluctuations of the gauge field. The latter's coupling to the Higgs field stabilizes the symmetric vacuum φ ≡ 0 at the local maximum of U as a parametric resonance effect, as long as a is not too large. Eventually, when a exceeds a critical value a EW , the Higgs field will begin to roll down towards a minimum of U , breaking the SO(4) symmetry. The corresponding time t EW signifies the electroweak phase transition in the early universe. This scenario was put forward recently by D. Friedan [2].
The goal of the current paper is a stability analysis of these classical oscillating "cosmic" Yang-Mills fields. To begin with, Section 2 describes the geometry of S 3 and reviews the classical configurations (A µ , g µν ) in terms of Newtonian solutions (ψ, a) for the anharmonic oscillator pair (V, W ). To investigate arbitrary small perturbations of the gauge field departing from the time-dependent background A µ parametrized by the "gauge energy" E, Section 3 linearizes the Yang-Mills equation around it and diagonalizes the fluctuation operator to obtain a spectrum of time-dependent natural frequencies.
To decide about the linear stability of the cosmic Yang-Mills configurations we have to analyze the longtime behavior of the solutions to Hill's equation for all these normal modes. In Section 4 we employ Floquet theory to learn that their growth rate is determined by the stroboscopic map or monodromy, which is easily computed numerically for any given mode. We do so for a number of low-frequency normal modes and find, when varying E, an alternating sequence of stable (bounded) and unstable (exponentially growing) fluctuations. The unstable bands roughly correspond to the parametric resonance frequencies.
With growing "gauge energy" the runaway perturbation modes become more prominent, and some of them persist in the infinite-energy limit, where we detect universal natural frequencies and monodromies. A special role is played by the SO(4)-invariant fluctuation, which merely shifts the parameter E of the background. We treat it exactly and beyond the linear regime in Section 5. This "singlet" mode turns out to be marginally stable, i.e. it has a vanishing Lyapunov exponent. Its linear growth, however, gets limited by nonlinear effects of the full fluctuation equation, whose analytic solutions exhibit wave beat behavior. Finally, some explicit data for the first few natural frequencies are collected in an Appendix.

Cosmic Yang-Mills solutions
In order to describe the classical Yang-Mills solutions we need to develop some elements of the spatial S 3 geometry. Taking advantage of the fact that and we introduce a basis {L a } for a = 1, 2, 3 of left-invariant vector fields on S 3 generating the right multiplication on SU(2) and forming the su(2) L algebra It is dual to a basis {e a } of left-invariant one-forms on S 3 , i.e. e a (L b ) = δ a b , subject to de a + ε a bc e b ∧ e c = 0 and e a e a = dΩ 2 3 .

(2.3)
One may obtain this basis by expanding the left Cartan one-form Here, the group element g provides the identification map which sends the S 3 north pole (0, i) to the group identity 1 2 . We shall coordinatize S 3 by an SU(2) group element g. The su(2) R half of the three-sphere's so(4) isometry is provided by right-invariant vector fields R a belonging to the left multiplication on the group manifold and obeying The differential of a function f on I × S 3 is then conveniently taken as Functions on S 3 can be expanded in a basis of harmonics Y j (g) with 2j ∈ N 0 , which are eigenfunctions of the scalar Laplacian, 1 where L 2 = L a L a and R 2 = R a R a are (minus four times) the Casimirs of su(2) L and su(2) R , respectively, The left-right (or toroidal) harmonics Y j;m,n are eigenfunctions of L 2 = R 2 , L 3 and R 3 , i 2 L 3 Y j;m,n = m Y j;m,n and i 2 R 3 Y j;m,n = n Y j;m,n (2.10) for m, n = −j, −j+1, . . . , +j, and hence the corresponding ladder operators The gauge potential is an su(2)-valued one-form on our spacetime. We use the I ×S 3 parametrization and write It has been shown [5,2] that the requirement of SO(4) equivariance enforces the form with some function ψ : I → R, where T a denotes the su(2) generators subject to so that the adjoint representation produces tr(T a T b ) = −8 δ ab . The corresponding field strength reads whereψ ≡ ∂ τ ψ. The Yang-Mills action on this ansatz simplifies to where I = [0, T ] and g here denotes the gauge coupling. Due to the principle of symmetric criticality [9], solutions to the mechanical problemψ + V (ψ) = 0 (2.18) will, via (2.14), provide Yang-Mills configurations which extremize the action. Conservation of energy implies that 1 2ψ and the generic solution in the double-well potential V is periodic in τ with a period T (E).
Hence, fixing a value for E and employing time translation invariance to setψ(0) = 0 uniquely determines the classical solution ψ(τ ) up to half-period shifts. Its explicit form is where cn and dn denote Jacobi elliptic functions, K is the complete elliptic integral of the first kind, and For E 1 2 , we have k 2 → 1 2 , and the solution is well approximated by 2 cos 2 , the unstable constant solution coexists with the celebrated bounce solution, and below it the solution bifurcates into oscillations in the left or right well of the double-well potential, which halfens the oscillation period. The two constant minima ψ = ±1 correspond to the vacua A = 0 and A = g −1 dg. Actually, the time translation freedom is broken by the finite range of I, so that time-shifted solutions differ in their boundary values ψ(0) and ψ(T ) and also in their value for the action. The corresponding color-electric and -magnetic field strengths read which yields a finite total energy (on the cylinder) of 6π 2 E/g 2 and a finite action [8,10] The energy-momentum tensor of our SO(4)-symmetric Yang-Mills solutions is readily found as which is traceless as expected.
The Einstein equations for a closed FLRW universe with cosmological constant Λ reduce to two independent relations, which can be taken to be its trace and its time-time component. In conformal time one gets, respectively, with a gravitational coupling κ = 8πG, a gravitational energy E and a cosmological potential The two anharmonic oscillators, with potential V for the gauge field and potential W for gravity, are coupled only via the balance of their conserved energies, which is nothing but the Wheeler-DeWitt constraint H = 0. Figure 2: Plots of the cosmological potential W (a) for Λ=1 and the double-well potential V .
The Friedmann equation (2.25), being a mechanical system with an inverted anharmonic potential (2.26), is again easily solved analytically, where we abbreviated 2 , we have k 2 → 1 2 , and the solution is well approximmated by We only listed solutions with initial value a(0) = 0 (big bang). There exist also (for E < 3 8Λ ) bouncing solutions, where the universe attains a minimal radius a min = 3 extension in the far past (t=−∞ ↔ τ =0) and the far future (t=+∞ ↔ τ =T ). For E >0 they are obtained by sending dn → −dn in (2.28) above. The quantity T listed there is the (conformal) lifetime of the universe, from the big bang until either the big rip (for E > 3 8Λ ) or the big crunch of an oscillating universe (for E < 3 8Λ ). The solution relevant to our Einstein-Yang-Mills system is entirely determined by the Newtonian energy E characterizing the cosmic Yang-Mills field: above the critical value of the universe expands forever (until t max =∞), while below this value it recollapses (at t max = T 0 dτ a(τ )). It demonstrates the necessity of a cosmological constant (whose role may be played by the Higgs expectation value) as well as the nonperturbative nature of the cosmic Yang-Mills field, whose contribution to the energy-momentum tensor is of O(g −2 ).

Natural perturbation frequencies
Our main task in this paper is an investigation of the stability of the cosmic Yang-Mills solutions reviewed in the previous section. For this, we should distinguish between global and local stability. The former is difficult to assess in a nonlinear dynamics but clear from the outset in case of a compact phase space. The latter refers to short-time behavior induced by linear perturbations around the reference configuration. We shall look at this firstly, in the present section and the following one. Here, we set out to diagonalize the fluctuation operator for our time-dependent Yang-Mills backgrounds and find the natural frequencies.
Even though our cosmic gauge-field configurations are SO(4)-invariant, we must allow for all kinds of fluctuations on top of it, SO(4)-symmetric perturbations being a very special subclass of them. A generic gauge potential "nearby" a classical solution A on I × S 3 can be expanded as with, using (µ) = (0, a), on which we notice the following actions (supressing the τ and g arguments), where the L a matrix elements are determined from (2.10) and (2.12), and S a are the components of the spin operator. The (metric and gauge) background-covariant derivative reads which is equivalent to The background A obeys the Coulomb gauge condition, but we cannot enforce these equations on the fluctuation Φ. However, we may impose the Lorenz gauge condition, which is seen to couple the temporal and spatial components of Φ in general. We then linearize the Yang-Mills equations around A and obtain with the Ricci tensor R µ0 = 0 and R ab = 2δ ab .
After a careful evaluation, the µ=0 equation yields while the µ=a equations read 11) It is convenient to package the orbital, spin, isospin, and fluctuation triplets into formal vectors, respectively, but they act in different spaces, hence on different indices, such that S 2 = T 2 = −8 on Φ.
In this notation, (3.7), (3.10) and (3.11) take the compact form (suppressing the color index p) A few remarks are in order. First, except for the last term, (3.14) is obtained from (3.15) by setting S = 0, since Φ 0 carries no spin index. Second, both equations can be recast as 17) which reveals a problem of addition of three spins and a corresponding symmetry under Third, for constant backgrounds (ψ=0) the temporal fluctuation Φ 0 decouples and may be gauged away. Still, the fluctuation operator in (3.17) is easily diagonalized only when the coefficient of one of the first three spin-squares vanishes, i.e. for L=0 (j=0), for the two vacua ψ=±1, or for the "meron" (or "sphaleron") ψ=0. The latter case has been analyzed by Hosotani and by Volkov [3,7].
Let us decompose the fluctuation problem (3.13)-(3.15) into finite-dimensional blocks according to a fixed value of the spin j ∈ 1 2 N, and suppress the j subscript. We employ the following coupling scheme, 3 Clearly, U and V act on Φ in su (2) representations j ⊗ 1 and j ⊗ 1 ⊗ 1, respectively. On Φ 0 , we must put S=0 and have just V = U act in a j ⊗ 1 representation. Combining the coupled equations (3.14) and (3.15) to a single linear system for (Φ p µ ) = (Φ p 0 , Φ p a ), we get a 12(2j+1)×12(2j+1) fluctuation matrix Ω 2 (j) , Actually, there is an additional overall (2j+1)-fold degeneracy present due to the trivial action of the su(2) R generators R a , which plays no role here and will be suppressed. Roughly speaking, the 3(2j+1) modes of Φ 0 are related to gauge modes, 4 and we still must impose the gauge condition (3.13), which also has 3(2j+1) components. Therefore, a subspace of dimension 6(2j+1) inside the space of all fluctuations will represent the physical gauge-equivalence classes in the end.
Our goal is to diagonalize the fluctuation operator (3.21) for a given fixed value of j. It has a block structure, whereN and N are given by the left-hand sides of (3.16) and (3.17), respectively. We introduce a basis where U 2 , V 2 and V 3 are diagonal, i.e.
23) and denote the irreducible su(2) v representations with those quantum numbers as v u . On the Φ 0 subspace, u is redundant since u=v as S=0. Working out the tensor products, we encounter the values v u with some representations obviously missing for j<2.
Let us treat theψ T term in (3.22) as a perturbation and momentarily put it to zero, so that Ω 2 (j) is block-diagonal for the time being. Then, it is easy to see from  2 is not diagonal in our basis. Therefore, we have a degeneracy in m. Furthermore, bothN and N decompose into at most three respectively five blocks with fixed values of v ranging from j−2 to j+2 and separated by semicolons in (3.24). Moreover, the N blocks are irreducible and trivially also carry a value of u=v. In contrast, N is not simply reducible; its V representations have multiplicity one, two or three. Only the N blocks with extremal v values in (3.24) are irreducible. The other ones are reducible and contain more than one U representation, hence the u-spin distinguishes between their (two or three) irreducible v subblocks. The only non-diagonal term in N is the ( L+ S) 2 contribution, which couples different copies of the same v-spin to each other, but of course not to any u=v block ofN , and does not lift the V 3 =m degeneracy. As a consequence, the unperturbed fluctuation equations for Φ 0 =Φ (v) and Φ=Φ (v,α) take the form (suppressing the m index) where 1 (v) denotes a unit matrix of size 2v+1, and α counts the multiplicity of the v-spin representation in N (between one and three). According to (3.16) the unperturbed frequency-squares forN are the eigenvaluesω with multiplicity 2v+1, hence we get Considering N in (3.17), we can read off the eigenvalues at v = j±2 because in these two extremal cases ( L+ S) 2 = U 2 is already diagonal in the |uvm basis. For the other v-values we must diagonalize a 2×2 or 3×3 matrix to find , ω 2 (j−1,α) = two roots of Q j−1 (λ) , ω 2 (j,α) = three roots of Q j (λ) , ω 2 (j+1,α) = two roots of Q j+1 (λ) , ω 2 (j+2) = root of Q j+2 (λ) = 2(2j+3) ψ + 2(2j 2 +6j+5) , Let us now turn on the perturbationψ T , which couples N withN , and consider the characteristic polynomial P j (λ) of our fluctuation problem, where we made use of Since T furnishes an su(2) representation (and not an intertwiner) it must be represented by square matrices and thus cannot connect different v representations. Hence the perturbation does not couple different v sectors but only links N andN in a commonv=v sector. Therefore, it does not affect the extremal sectors v = j±2. Moreover, switching to a diagonal basis {|αvm } for N we can simplify to Observing that T |α α| T (v) = −tv ,α T 2 (v) = 8 tv ,α 1 (v) with some coefficient functions tv ,α (ψ), with α tv ,α = 1, we learn that the V 3 degeneracy remains intact and arrive at (v ∈ {j−1, j, j+1}) where Pv = Qv α tv ,α (ω 2 (v,α) −λ) −1 is a polynomial of degree one less than Qv since all poles cancel, and Rv is a polynomial of one degree more. We list the polynomials Q v , Pv and Rv for j≤2 in the Appendix.
We still have to discuss the gauge condition (3.13), which can be cast into the form with a 3(2j+1)×7(2j+1) linear (in ψ) matrix function K. 6 Here the v sum runs over (j−1, j, j+1) only, since the gauge condition We conclude this section with more details for the simplest examples, which are constant backgrounds and j=0 backgrounds. For the vacuum background, say ψ = −1, which is isospin degenerate, one gets for u = j 4(j+1) 2 for u = j+1 (3.39) 6 We have to bring back the m indices because the gauge condition is not diagonal in them.
For a time-varying background, the natural frequencies Ω (j,v,β) inherit a τ dependence from the background ψ(τ ). Direct diagonalization is still possible for j=0, where we should solve It implies the unperturbed frequencies (suppressing the j index) as long asψ is ignored. There are no v-spin multiplicities (larger than one) here. Turning onψ and observing that ( T · Φ) p ∼ δ pb φ b , the characteristic polynomial of the coupled 12×12 system in the |uvm basis reads Specializing the general discussion above to j=0, we find just t 1 =1 so that P 1 =1 and arrive at

(3.47)
We see that the frequencies Ω 2 (0) =ω 2 (0) and Ω 2 (2) =ω 2 (2) are unchanged and given by (3.44), while the gauge modeω 2 (1) gets entangled with the (unphysical) v=1 mode to produce the pair Ω 2 (1,±) = 1 2 (ω 2 (1) +ω 2 (1) ) ± 1 4 (ω 2 (1) +ω 2 (1) ) 2 −ω 2 (1) ω 2 (1) + 8ψ 2 = 3ψ 2 +3ψ+2 ± ψ 2 (ψ−1) 2 + 8ψ 2 (3.48) with a triple degeneracy. There are avoided crossings at ψ=0 and ψ=1. Removing the unphysical and gauge modes in pairs, we remain with the singlet mode Ω 2 (0,0) and the fivefold-degenerate Ω 2 (0,2) . For all higher spins j>0, analytic expressions for the natural frequencies Ω (j,v,β) now require merely solving a few polynomial equations of order four at worst. We have done so up to j=2 and list them in the Appendix but refrain from giving further explicit examples here. Below we display the cases of j=0 and j=2, with similar coloring for like v values, whose curves avoid crossing each other. One can see that some of the normal modes dip into the negative regime, i.e. their frequency-squares become negative, for a certain fraction of the time τ . Because of this and, quite generally, due to the τ variability of the natural frequencies, it is not easy to predict the long-term evolution of the fluctuation modes. Clearly, the stability of the zero solution Φ≡0, equivalent to the linear stability of the background Yang-Mills configuration, is not simply decided by the sign of the τ -average of the corresponding frequency-square.

Stability analysis: stroboscopic map and Floquet theory
The diagonalized linear fluctuation equation (3.34) represents a bunch of Hill's equations, where the frequency-squared is a root of a polynomial of order up to four with coefficients given by a polynomial of twice that order in Jacobi elliptic functions. A unique solution requires fixing two initial conditions, and so for each fluctuation Φ (j,v,β) there is a two-dimensional solution space. It is well known that Hill's equation, e.g. in the limit of Mathieu's equation, displays parametric resonance phenomena, which can stabilize otherwise unstable systems or destabilize otherwise stable ones.
For oscillating dynamical systems with periodically varying frequency, there exist some general tools to analyze linear stability. Switching to a Hamiltonian picture and to phase space, it is convenient to transform the second-order differential equation into a system of two coupled first-order equations (suppressing all quantum numbers), where the frequency Ω(τ ) is T -periodic (sometimes T 2 -periodic) in τ . The solution to this first-order system is formally given by where T denotes time ordering. Because of the time dependence of Ω, the time evolution operator above is not homogeneous thus does not constitute a one-parameter group, except when the propagation interval is an integer multiple of the period T . For τ =T , one speaks of the stroboscopic map [11] M := T exp The linear map M is a functional of the chosen background solution ψ and hence depends on its parameter E or k. This background is Lyapunov stable if the trivial solution Φ≡0 is, which is decided by the two eigenvalues µ 1 and µ 2 of M . Since the system is Hamiltonian, det M =1, and we have three cases:  dτ n H n (τ 1 , τ 2 , . . . , τ n ) Ω 2 (τ 1 ) Ω 2 (τ 2 ) · · · Ω 2 (τ n ) with H n (τ 1 , τ 2 , . . . , τ n ) = (τ 1 −τ 2 )(τ 2 −τ 3 ) · · · (τ n−1 −τ n )(τ n −τ 1 +1) and H 1 (τ 1 ) = 1 .
It is impossible to evaluate the integrals M n without explicit knowledge of ω 2 (x). As a crude guess, we replace the weight function by its (constant) average value and obtain which yields This expression indicates stability as long as ω 2 > 0. However, the result for the j=0 singlet mode ω 2 = Ω 2 (0,0) in (4.11) below already shows that the averaged frequency-squared may turn negative in certain domains thus changing the cos into a cosh there.
To do better, let us look at the individual terms M n in (4.7) for the simplest case of the SO(4) singlet fluctuation, i.e. Ω 2 (0,0) = 6ψ 2 −2 in (3.44). Its average frequency-square is easily computed to be where E(k) and K(k) denote the second and first complete elliptic integrals, respectively. Plotting this expression as a function of the modulus k, we see that it becomes negative only in a very narrow range around k=1, namely for |k−1| 0.00005. We have only been able to analytically evaluate (with k<1 for and which does not suffice to rule out instability. Indeed, numerical studies show that M n as a function of k looses its positivity in a range around k=1 which increases with n, where the series (4.7) ceases to be alternating. Moreover, even in the limit of a very large background amplitude, k 2 → 1 2 , we find that implying that we must push the series in (4.7) at least to O(M 10 T 20 ), even though it turns out that M n < Ω 2 (0,0) n at k 2 = 1 2 for n > 1. For a more complete analysis of linear stability in an oscillating system with time-dependent frequency we can take recourse to Floquet theory. It tells us that a general fundamental matrix solution 14) of our system (4.1) with some initial condition Φ(0) = Φ 0 can be expressed in so-called Floquet normal form as where Q(τ ) and R are real 2×2 matrices, so that the time dependence of the frequency can be transformed away by a change of coordinates, Due to the identity we see that our stroboscopic map M is nothing but the monodromy, and so that its eigenvalues (or characteristic multipliers) define a pair of (complex) Floquet exponents ρ i whose real parts are the Lyapunov exponents. Since µ 1 µ 2 = 1 implies that ρ 1 +ρ 2 = 0, our system is linearly stable if and only if both eigenvalues ρ i of R are purely imaginary (or zero).
Generally it is impossible to find analytically the monodromy pertaining to a normal mode Φ (j,v,β) . 7 However, we can evaluate it numerically for a number of examples. Before doing so, let us estimate at which energies E or, rather, moduli k, possible resonance frequencies might occur. To this end, we determine the period-average of the natural frequency Ω (j,v,β) and compare it to its modulation frequency 2π T . If we model dτ Ω 2 (τ ) and h(τ ) ∝ cos(2πτ /T ) , (4.20) where T = 4 K(k), then the resonance condition is met for Since this model reproduces only the rough features of Ω 2 (τ ), we expect potential instability due to parametric resonance effects in a band around or near the values k .
Our expectation is confirmed by precise numerical evaluation of various monodromies as a function of k. Below we display, together with the would-be resonant values k , the function trM (k) for the sample cases of (j, v) = (2, 0) and (2,2). One sees that, on both sides of the critical value of E= 1 2 (or k=1), corresponding to the double-well local maximum, the k values accumulate at the critical point. But while for k>1 (energy below the critical point) trM (k) oscillates between values close to 2 in magnitude and thus exponential growth is rare and mild, for k<1 (energy above the critical point) the oscillatory behavior of trM (k) comes with an amplitude exceeding 2 and growing with energy. Hence, in this latter regime stable and unstable bands alternate. This is supported by long-term numerical integration, as we demonstrate in Figure 9 by plotting Φ(τ ) for (j, v, β) = (2, 2, 1) with initial values Φ(0)=1 andΦ(0)=0 on both sides of the first transition from instability to stability for tr M (2,2,1) shown in Figure 8 (at the highest value of E or the lowest value of k).
Most relevant for the cosmological application is the regime of very large energies, E → ∞ (or k → 1/ √ 2). In this limit, we observe the following universal behavior. Because the period T collapses with = k 2 −1/2, we rescale 7 An exception is the SO(4) singlet perturbation Φ (0,0) , to be treated in the following section. so that the tilded quantities remain finite in the limit, and find, withω 2 because all j-dependent terms in the polynomials are subleading and drop out in the limit. Factorizing the R polynomials, we find the four universal natural frequency-squares One must pay attention, however, to the fact that the avoided crossings disappear in the → 0 limit. Therefore, the correct limiting frequencies to input into

Singlet perturbation: exact treatment
Even though the Floquet representation helped to reduce the long-time behavior of the perturbations to the analysis of a single period T , it normally does not give us an exact solution to Hill's equation. However, for the SO(4) singlet fluctuation around ψ(τ ), we can employ the fact thatψ trivially solves the fluctuation equation, with a frequency function which is T 2 -periodic. This implies that all fluctuation modes are T -periodic. With the knowledge of an explicit solution to the fluctuation equation we can reduce the latter to a first-order equation and solve that one to find a second solution. The normalizations are arbitrary, so we choose which are linearly independent since For simplicity, we restrict ourselves to the energy range 1 2 <E<∞, i.e. 1>k 2 > 1 2 . Explicitly, we have Φ 1 (τ ) = sn τ , k dn τ , k , where am(z, k) denotes the Jacobi amplitude and E(z, k) is the elliptic integral of the second kind. As can be checked, the initial conditions are which fixes the ambiguity of adding to Φ 2 a piece proportional to Φ 1 . Hence,  Figure 11: Plot of the SO(4) singlet fluctuation modes Φ 1 and Φ 2 over eight periods for k 2 =0.81.
We know that Φ 1 ∼ψ is T -periodic, and so isΦ 1 , but not the second solution, where the integral diverges at the turning points and must be regularized by subtracting the Weierstraß ℘ function with the appropriate half-periods. Since Φ 1 has periodic zeros, Φ 2 does return to −1 at integer multiples of T . It follows that the Φ 2 oscillation linearly grows in amplitude with a rate (per period) of which is always larger than 7.629, attained at k ≈ 0.882. In essence, we have managed to compute the monodromy and thus easily obtain the Floquet representation, Obviously, we have encountered a marginally stable situation, since M is of parabolic type. There is no exponential growth, and Φ 1 is periodic thus bounded, but Φ 2 grows without bound as long as one stays in the linear regime. Note that we never made use of the form of our Newtonian potential. In fact, this behavior is typical for a conservative mechanical system with oscillatory motion.
= δk ∂ k ψ(τ −δ ) + 1 2 (δk) 2 ∂ 2 k ψ(τ −δ ) + 1 6 (δk) 3 ∂ 3 k ψ(τ −δ ) + . . . , (5.17) because ∂ ψ = −ψ. Clearly, a shift in only shifts the time dependence of the frequency and does not alter the energy E, which is not very interesting. Its linear part corresponds to the mode Φ 1 ∼ψ of the previous section. A change in k, in contract, will lead to a solution with an altered frequency and energy. Its linear part is given by Φ 2 , which grows linearly in time. However, due to the boundedness of the full motion, the nonlinear corrections have to limit this growth and ultimately must bring the fluctuation back close to zero. This is the familiar wave beat phenomenon: the difference of two oscillating functions, ψ and ψ, with slightly different frequencies, will display an amplitude oscillation with a beat frequency given by the difference. This is borne out in the following plots. As a result, we can assert a long-term stability of the cosmic Yang-Mills fields against the SO(4) singlet perturbation, even though on shorter time scales an excursion to a nearby solution is not met with a linear backreaction.

Conclusions
We have revisited a cosmological scenario recently put forward by Friedan [2] and based on the Standard Model plus gravity alone. An SO(4) symmetric sector is analytically solvable and reduces to three coupled anharmonic oscillators (for the metric, an SU(2) Yang-Mills field and the Higgs field, the latter being frozen to its vacuum state). We have presented a complete analysis of the linear gauge-field perturbations of the time-dependent Yang-Mills solution, by diagonalizing the fluctuation operator and studying the long-time behavior of the ensuing Hill's equations using the stroboscopic map and Floquet theory. For parametrically large gauge-field energy (as is required in Friedan's setup) the natural frequencies and monodromies become universal, and some unstable perturbation modes survive even in this limit. This provides strong evidence that such oscillating cosmic Yang-Mills fields are unstable against small perturbations, although we have not yet included metric fluctuations here. Their influence will be analyzed in follow-up work.