{'fraction_non_alphanumeric': 0.03573322023115258, 'fraction_numerical': 0.012750503658148659, 'mean_word_length': 4.528282532239156, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gravitational Instability of Yang-Mills Cosmologies

The gravitational instability of Yang-Mills cosmologies is numerically studied with the hamiltonian formulation of the spherically symmetric Einstein-Yang-Mills equations with SU(2) gauge group. On the short term, the expansion dilutes the energy densities of the Yang-Mills fluctuations due to their conformal invariance. In this early regime, the gauge potentials appear oscillating quietly in an interaction potential quite similar to the one of the homogeneous case. However, on the long term, the expansion finally becomes significantly inhomogeneous and no more mimics a conformal transformation of the metric. Thereafter, the Yang-Mills fluctuations enter a complex non-linear regime, accompanied by diffusion, while their associated energy contrasts grow.


Introduction
The spherically symmetric Einstein-Yang-Mills (EYM) equations with SU(2) gauge group has been extensively studied in the past decades, from cosmological solutions (cf. [1,2,3,4,5]) to static configurations (cf. [6,7]). The homogeneous and anisotropic cosmological models with Yang-Mills (YM) fields were considered in [8,9,10] and exhibit chaotic features. For a detailed review of the litterature prior to 1999, see [7] and references therein. More recently, some authors have investigated non-abelian Born-Infeld cosmology [11] in which the YM-Born-Infeld lagrangian is no more scale invariant and therefore yields to different behaviour from classical YM fields, especially near the singularity.
However, most of the work done so far provides (semi-)analytical and/or numerical solutions in the case of only one independent variable: time (cosmological solutions) or radius (gravitating non abelian solitons and black holes). Nevertheless, some authors [5] have proposed time-dependent inhomogeneous analytical solutions of EYM theory but at the cost of strong restrictions on the form of the metric and the gauge potentials. As well, perturbation methods have been used several times to study the instability of self-gravitating non-abelian solutions [7]. Here, we use the hamiltonian formulation of spherically symmetric EYM equations with SU(2) gauge group to build a numerical method that allows us to find numerical time-dependent inhomogeneous solutions under a few general assumptions. In this paper, we apply it to the question of gravitational instability of YM cosmologies, by studying collapsing shells, which was done for Minkowski space in [12].
In the case of homogeneous and isotropic spacetimes, the conformal invariance of YM fields allows to solve separately Einstein and YM equations as the resulting solution is a radiation dominated universe. The general ansatz for open, flat and closed cosmologies filled with su(2)−valued YM fields was proposed in [1], for which the solution in [3] is a special case. In this paper, we will investigate departures from this ansatz, consequent to some initial perturbations, in order to characterize their gravitational instability. Although, our approach differs from usual perturbation theory in the sense that we do not linearize the equations, we consider here time-dependent inhomogeneous numerical solutions around the general ansatz. One will argue that the conformal invariance of YM fields implies that any primeval excitation of these fields will be diluted as the universe expands. We will see that this is true only in the first times of the evolution of the fluctuations, when the expansion is not too inhomogeneous and looks like a conformal transforma-tion of the metric. Thereafter, interesting phenomena such as diffusion and self-interaction arise which results in growing instabilities.
In section 2, we remind the reader about the hamiltonian formulation of the spherically symmetric EYM system with SU(2) gauge group, firstly introduced in [13], then we recall the Gal'tsov-Volkov ansatz that we will use as a background solution. Thereafter, the numerical method used to integrate the hamiltonian formulation is introduced in section 3, a test on the homogeneous ansatz being provided. The section 4 is devoted to a qualitative discussion of the gravitational instability from the numerical results presented in the figures at the end of the paper. Finally, we conclude by some comments on the application of our results to cosmology as well as some perspectives.

Hamiltonian Formulation
The hamiltonian approach of the gravitational field equations is the wellknown Arnowitt-Deser-Misner formalism [14] which can be completed by a hamiltonian treatment of YM equations describing a gauge interaction [13,15]. The spherically symmetric inhomogeneous vacuum space-times were studied in [16,17] while the spherically symmetric EYM equations with gauge group SU(2) were first discussed in [13]. We briefly recall here the hamiltonian equations of the EYM system and fix the gauges in which we will work for the rest of this paper. The coupled EYM action can be written as where κ = 8πG (c = 1), g is the determinant of the 4-metric, R is the scalar curvature, F a µν are the components of the YM field strength tensor 1 . For a particular gauge group G locally specified by its Lie algebra [T a , T b ] = if abc T c , the field strength tensor is given through its components by s are the components of the one-form connection A = A µ dx µ = A a µ T a dx µ · Here, we will focus on the gauge group SU(2) with generators T a = 1 2 τ a where the τ a are the usual Pauli matrices and the structure constants are given by the completely antisymmetric tensor f abc = ǫ abc . The normalization is T r(T a T b ) = 1 2 δ ab · Assuming spherical symmetry of the spatial slices of the universe, we write the metric in a 3+1 decomposition as follows : where dΩ 2 = dθ 2 + sin 2 θ dϕ 2 is the solid angle element, N and N χ are the lapse and shift function respectively 2 , all these components N, N χ , µ, λ being functions of the coordinates χ and t.
If the gauge fields A µ exhibit spherical symmetry up to a gauge (see [13,18,4]), they can be written thanks to the so-called W itten ansatz [19] (with the same conventions as in [6]): where a, b, c, d are all functions of both coordinates χ and t. In Eq. (3), only three components are of physical relevance since the potential d can be gauged away by a transformation that preserves the SO(3) symmetry (see [1,11] and references therein for a detailed discussion).
If we now let the spatial components of the metric g ij and the connection A a i (i = 1, 2, 3) be canonical variables and if we consider the lapse and shift functions N, N χ as well as the electric components A a 0 of the YM connection as being Lagrange multipliers, it is possible to define some momenta conjugated to these variables and then to rewrite the coupled EYM action Eq.(1) under constrained hamiltonian form (cf. [13,14,15] and references therein): where π ij = (3) g g ij K l l − K ij (K ij being the extrinsic curvature tensor), In Eq.(4), the generators of normal and tangential deformations (H ⊥ and H i , respectively) and those of internal gauge transformations G a are given by 3 with R the scalar curvature of the 3-metric and π = π i i , | denoting the covariant derivative w.r.t. the 3-metric and the vector B i a reading Let us now focus on the spherically symmetric hamiltonian formalism. For the gravitational part of the EYM system, the canonical variables are chosen as 4 g ij = diag e 2µ , e 2λ , e 2λ while their conjugate momenta are given by π ij = diag e −2µ π µ 2 , e −2λ π λ 4 , e −2λ π λ 4 . Variation of Eq.(4) w.r.t. the lapse and shift functions N and N χ will provide the hamiltonian and supermomentum constraints while the variation w.r.t. the canonical variables µ, λ and their conjugate momenta π µ and π λ will lead to the Hamilton equations. As we are interested in perturbations of homogeneous solutions and as the YM lagrangian in Eq.(1) is conformally invariant, we will work in the conformal gauge for the gravitational field: N(χ, t) = R(t) where R(t) is the scale factor given by the Friedmann equation and N χ (χ, t) = 0. After the variation and the choice of gauge, we now have for the dynamics of the gravitational field the following constraints: and the following Hamilton equations, where T ν µ are the components of the YM stress-energy tensor (see further), and where a prime and a dot denote derivations w.r.t. position χ and time t respectively. On the YM side, we have, according to Eq.(3), only one electrical component A 3 0 = a and three couples of variables and momenta (A 3 1 , π 1,3 Varying the YM action w.r.t. these variables and fixing the gauge d(χ, t) = 0 so that π d = e µ N ac leads to the following constraint and to the following Hamilton equationṡ With these variables, we can now write the components of the YM stressenergy tensor T µν = −2F a µα F αa in order to achieve the coupling between gauge and gravitation sectors. Eqs. (9)(10)(11)(12)(13)(14)(15)(16)(17) are the hamiltonian formulation of the well-known spherically symmetric EYM equations with gauge group SU(2). In the next paragraph, we recall the solution for those fields in a Friedmann-Lemaître universe before describing the method used to study gravitational instability of such solutions.

Yang-Mills Cosmologies
Non−abelian YM fields can fill a non trivial, radiation dominated, Friedmann-Lemaître universe due to their increased number of degrees of freedom, and the couplings between them, compared to the abelian Maxwell theory. In particular, the case of gauge group SU(2) has been extensively studied, so we refer the reader to [1,3,2,4,11,7] and others references therein for a review.
With the definitions of previous paragraph, we can write the homogeneous solutions we will use for the background in which the fluctuations will evolve. First, the geometry is given by (in the conformal gauge N = R(t), upper script B stands for background): where Σ = sin χ, χ, sinh χ for closed, flat and open cosmologies respectively. The conformal invariance and the consequent vanishing trace of the YM stress-energy tensor result in a radiation-dominated universe for which the Friedmann equation Eq.(10) reads (using Eqs. (18)): with E some constant to be precised in one moment. This allows us to write the following solutions for the scale factor R(t) (see also [24]).
for flat, elliptic and hyperbolic universes, respectively. The constant C = ρR 4 is a first integral of the conservation equations ∇ µ T µ 0 = 0.
The gauge fields and their conjugate momenta are given by the following ansatz [1,11]: where k = +1, 0, −1 for closed, flat and open cosmologies respectively, and σ is a function of time. By substituting the ansatz for the gauge fields Eqs. (21) into Eqs. (12)(13), the YM equations now reduce tö or, expressed as a first integral where E 2 = 2 3 C = 2 3 ρR 4 . The conformal invariance of YM fields is now manifest as the scale factor R(t) does not appear in Eqs. (22,23). Therefore, in the conformal gauge, the time part of the gauge fields σ(t) oscillates in the potential well (σ 2 − k) 2 . In the synchronous gauge N = 1, it can be shown that the frequency of the oscillation decreases with expansion. In the rest of the paper, we will focus on the case of flat space as the observed features of gravitational instability can be deduced from the general properties of the homogeneous solution presented in this paragraph and can be therefore extrapolated to curved cosmologies.

Overview of the Method
In order to integrate the hamiltonian EYM equations, we will use a twostep numerical procedure often referred as free evolution method. First, we solve the three constraints Eqs.(9-10) and Eq.(12) -this consists of the initial value problem (IVP) -then we propagate the initial data of the fields imposed or found by the IVP with the Hamilton equations Eqs. (11)(12)(13) and check the accuracy and stability of evolution by evaluating the error made on the constraints. The growing violation of the constraints in such procedures and the development of formalisms reducing it are very hot topics in numerical relativity (see for example [20], and [21] for an introduction to numerical relativity). Here, we choose the canonical approach to the EYM system (original ADM formalism [14] for gravitation and classical treatment for YM [13]) which we will rewrite in section 3.2 in order to improve stability and accuracy.
But before going further, let us give some comments on the particular choice of gauges we made for both the gravitational and YM fields. For the first one, we chose a slightly modificated version of the famous geodesic slicing: N(χ, t) = 1, N i (χ, t) = 0, which we called the conformal gauge: N(χ, t) = R(t) and N i (χ, t) = 0. Even though the geodesic condition is well-known to produce coordinate singularities (cf. [22]), it is the most efficient for weak gravitational fields such as instabilities in a homogeneous background. The reason why we chose a conformal gauge with N = R(t) rather than a synchronous one N = 1 is because the YM fields oscillate at a fixed frequency in the first while their period is stretched with time in the second. Therefore, using a synchronous gauge would have increased the time over which the fluctuations around the homogeneous solution significantly change. And, as we used an explicit scheme, it was crucial to reduce this time in order to avoid numerical discrepancies. Furthermore, our choice of a gauge dependent formalism (d was set to 0) was prefered for two reasons. The first is that it will sweep away one degree of freedom in the hamiltonian formalism, which simplifies the equations and therefore lowers truncation errors as well as avoiding those related to the couple (d, π d ) of unknowns. The second is because it simplifies the algorithm. Indeed, with a gauge-invariant formalism, the electric potential A 3 0 = a which represents the freedom of making arbitrary gauge transformations in the YM fields during the evolution of the system has to be determined by an additionnal gauge condition. By fixing the gauge, we are left with a Hamilton equation for this component. However, by reducing in this way the number of degrees of freedom in the analysis, we cancel an additional scalar field in the theory, in the context of the two-dimensional abelian Higgs mechanism that consitutes the spherically symmetric EYM system with gauge group SU(2) (cf. [18,7]). But now, we have learned from a more deep analysis of our results that these scalar fields and the interactions related to it bring a good part of the energy of the fluctuations. Part of our future work will be to examine how energy spread out into the different modes of the fully gauge invariant YM equations and how it affects the results presented here.

Initial Value Problem (IVP)
In this problem, we have to determine the initial distribution of nine fields (µ 0 , λ 0 , π µ,0 , π λ,0 , a 0 , b 0 , c 0 , π b,0 , π c,0 ) 5 that verifies the three constraints (two for the gravitational field and one from the gauge sector) of the EYM hamiltonian system. Here, we will only write the formulas for a perturbation of a flat Friedmann-Lemaître universe as they can easily be generalized to closed and hyperbolic cosmologies. In order to do this, we fix six fields as follows: • We start working some time after the Big Bang to avoid primeval singularities : R(0) = 1. Please note that the formula given in section 2.2 for the scale factors should be translated accordingly.
• We set λ 0 = ln χ so that the initial circonferential distances are taken as references for the coordinate χ.
• We perturb the following fields around their homogeneous values : where ρ B 0 is the initial background density, taken as a parameter. The functions ǫ i (χ) are also provided and will be taken so as to have convenient boundary conditions for the fields. Typically, we will study the gravitational instability of shells in a Friedmann-Lemaître background, i.e. we will impose that the ǫ i 's vanish at the boundaries of the χ coordinate interval. For example, we can take for the ǫ i 's some gaussian functions.
From there, we solve the super-momentum constraint Eq.(9) w.r.t. to π λ,0 and the YM constraint Eq.(12) w.r.t. µ 0 and we substitute the results into the hamiltonian constraint to get the following differential equation for π b,0 : 5 Lower script zero means at the time t = 0, for example µ 0 = µ(χ, 0)· where we put φ = π b,0 (χ) = π b (χ, 0) and where the coefficients are The fields with lower script 0 mentionned in Eq. (26) are those given by hypothesis in Eq. (24). From the solution φ = π b,0 of Eq.(25), one can complete the set of initial data with the following relations We used a finite difference scheme to solve Eq.(25) as follows: defining χ i = χ M IN + ∆χ(i − 1) (i = 1, · · · , N χ + 1), φ i = φ(χ i ) and analogous definition for other fields, Eq.(25) becomes a non-linear system of N χ − 2 algebraic equations for determining the φ i 's: with the boundary conditions , σ 0 becoming another parameter setting the initial position in phase space for the time part of the gauge fields filling the flat background. We choose a Newton-Raphson method (cf. [23]) to solve Eq.(28) with the initial guess : φ i = π B b,0 (χ i ) =σ 0 χ 2 i and stop the relaxation when the error on the equation is saturated.

Propagation of Hamilton Equations
In this paragraph, we first rewrite the EYM constraints and Hamilton equations in order to render the numerical integration more accurate and stable. Indeed, a direct integration of Eqs. (11) and Eqs. (13) will have to recover first the homogeneous solution Eqs. (18)(19)(20)(21) before focusing to its perturbation. After having first directly integrated Eqs. (11) and Eqs.(13), we found it more interesting to rewrite the equations so that they would be more adapted to the study of the gravitational instability and less dependent on the homogeneous solution. The change of variables explained below results in a faster and more accurate method.
Let us again assume a flat Friedmann-Lemaître background and proceed to the following change of variables: so that the homogeneous solutions (m B , l B , π B m , π B l , α B , β B , γ B , π B β , π B γ ) are now a set of constants w.r.t. χ : (1, 1,Ṙ,Ṙ,σ, σ 3 , σ 2 ,σ, σσ)· In terms of these new variables, the hamiltonian formulation of EYM system Eqs. (9)(10)(11)(12)(13)(14)(15)(16)(17) now become for the constraints and for the Hamilton equationṡ with the following definitions for the components of the stress-energy tensor is the same expression as T 0 0 except for the sign of the last two terms). Finally, the YM equations can be writteṅ In order to integrate Eqs.(32-34), we used a second order difference scheme in both χ and t (leapfrog-like) and, for the first step only, a forward time centered space scheme (F T CS, cf. [23]). This explicit scheme is easy to implement but yields a very strong Courant condition in this case. Anyway, for the study of gravitational instability during a few oscillations of σ in Eq. (22), it will be sufficient. Another crucial point is the question of boundary conditions, already mentionned in the previous paragraph. As we focus on shells and all our initial perturbations have been chosen to vanish rapidly on the edges of the considered interval in χ, we impose that the boundary conditions follow the homogeneous evolution. For example, a field f (χ, t) ≡ f i,j will have as boundary conditions where the upper script B stands for the background solution Eqs.(18-21).

Test on the homogeneous solution
In this paragraph, we briefly discuss the precision of the algorithm presented before. Figures 1 to 3 show the evolution of the mean quadratic error on the various fields used in the computation. The general divergent behaviour is due to the simple explicit scheme we have used. Therefore, the scheme is unstable, partially due to the violation of the boundary conditions that introduce truncation errors on the edges of the grid (see next section). We found a very strong Courant condition, for example ∆χ ≈ 10 −1 ; ∆t ≈ 10 −6 to have an acceptable precision on the time scale of an oscillation period of the function σ B (measuring the period of the gauge fields). For the gauge variables α, β and γ, we see in Figure 1 that E α > E β > E γ , while for the metric variables E m > E l (Figure 2). For the constraints (Figure 3), the boundary effects dominate the quadratic error on the hamiltonian and super-momentum constraints H and H 1 though the errors are much lower away from the edges of the χ−interval (see also Figure 15). These boundary effects are present in all the fields in the computation (see also other 3D plots). In Figure 3 (iii), we see that the error on the gauge constraint G 3 is much lower than those of Einstein's part.

Discussion
In this section, we use the numerical results obtained from the method introduced before in order to present some interesting features of the gravitational instability of YM cosmologies. For the sake of simplicity, we used for the initial pertubations a gaussian function of the form: where ǫ i , χ F and w are parameters indicating respectively the amplitude, position and width of the fluctuation. The last is related to the coherence length λ P of the fluctuation, the length over which the perturbation changes significantly. The figures were computed with only one field initially perturbed amongst the set (a 0 , b 0 , c 0 ,ċ 0 ,λ 0 ). Perturbing several fields will change the amplitude of density and pressures fluctuations on the initial space-like slice specified by the IVP and will affect the developments of the resulting instabilities in the long term. Indeed, as we will see further, the first stages of the evolution are characterized by a conformally-invariant regime maintained by the too less inhomogeneous expansion. The instabilities arise only after some time, depending on the initial conditions, when the fields reach a sufficient amount of inhomogeneity. Figures 5, 6, 7, 8, 11 and 15 show a typical evolution for the metric, the gauge potentials and the YM stress-energy tensor components as well as the corresponding values of the three constraints under an initial perturbation of the gauge field c. It is useful to notice that only 50 iterations have been plotted in each 3D plot.
For the gravitational sector, the perturbations of the spatial coefficients of the metric grow up as a power of conformal time (see Figure 5). This is typical of a radiation dominated universe and has been studied for a long time by several authors such as Lemaître, Tolman, Lifshitz and Bonnor (for a review, see [25]). From there, one should expect that the conformal invariance of YM lagrangian should result in a rather independent dynamics of the gauge fields regards to the expansion. The last should thus dilute the energy densities of the YM initial fluctuations. In fact, this is true when space-time is not significantly inhomogeneous i.e. when the perturbations of metric coefficients are negligible. In this case, the slightly inhomogeneous expansion mimics a conformal transformation of the metric. However, when the metric perturbations become important, we observe that the gauge fields start diffusing and the pressures and density contrasts start growing.
Before going further, let us examine more precisely the evolution of the gauge potentials. One can see in Figures 6 to 8 that the gauge fields seem to roll on the potential well V = σ 4 of the homogeneous solution but this potential is now distorded by the perturbations. For example, the electric component a, related to the speed of the oscillation, reaches its maxima at the middle of the phase space orbit σ B = 0 (see Figure 4), while magnetic components b and c oscillate at almost the same frequency of the background solution σ B · Let us give an idea of the distortion of the potential. Discarding the kinetic terms into Eq. (14), this potential can be written, in terms of the variables introduced in section 3.2: In Figure 9, we have represented the potential V associated to the results in Figures 5,6,7,8,11 and 15. The oscillation has been unfolded along the σ−axis and the resulting plot has been split in two for sake of readiness. The distortion has also been scaled up to be visible (the function plotted is σ 4 + 6 × 10 3 (V − σ 4 ), the other parameters are same as in Figure 4). The potential of inhomogeneous YM cosmologies appears to be an implicit function of time.
Furthermore, the gauge fields diffuse along the χ−axis, the perturbations spreading out with time. This diffusion is due to the increasing importance of spatial gradients of the gauge fields and the inhomogeneity of the metric in the YM equations. Another way to see this is to draw the evolution of the gauge fields during half a period of σ B for the same amplitude of the initial perturbation but for different values of the background initial expansion rate (parametrized here by ρ B 0 ). This is done in Figure 10 for an initial perturbation in the electric component a. For small values of the expansion rate ( Figure 10 (i)), the diffusion is important but it rapidly disappears for stronger expansion (Figure 10 (ii) and (iii)). A much clearer way to express this is to compare the size of the fluctuations, given roughly by their coherence length λ P , to the Hubble length L H = 1/H = t + 3 κC (c = 1). At the beginning of the evolution (t = 0), the last is given by Therefore, by varying the parameter ρ B 0 and keeping the same coherence length through the same parameter w, we just modify the ratio between both quantities. The corresponding initial Hubble lengths have been specified in the caption of Figure 10 with a coherence length λ P ≈ 1. We can thus say that the diffusion is initially important for the short wavelengths λ P ≈ L H and, for wavelengths greater than the Hubble length, becomes important only when the spatial inhomogeneity is strong enough, as in Figures  5, 6, 7, 8 and 11.
The observable quantities are represented in Figure 11 under the form of the density contrast: and similar relations for the radial (T 1 1 ) and tangential (T 2 2 ) pressures contrasts. It should also be noted that, according to the conformal invariance of YM fields, the scalar curvature R vanishes exactly and therefore the components of the stress-energy tensor are directly proportional to those of the Ricci tensor through Einstein equations and we have therefore a direct information on the spatial inhomogeneities. These contrasts oscillates in general, reflecting the underlying oscillations of gauge potentials, even though the amplitude of the oscillation is damped in the first stages of the evolution and then grow after some time, with a different regime for density and pressures contrasts. For example, in the case illustrated in Figure 11, the density constrast first decays, then enter a growing regime while the amplitude of the (radial and tangential) pressures contrasts slowly decreases.
According to the radiation nature of YM, one should have expected that the evolution of the density contrast would have been a superposition of a growing mode in squared power of the conformal time 6 and a decaying one in inverse squared power of conformal time. This is well-known from a simple first order perturbation computation (see [25]). In fact, when the perturbations of the geometry are small and their gradients negligible, the expansion can be considered as roughly homogeneous and it dilutes the energy of the conformally invariant YM fields. This case is a pretty good approximation of the early universe and is represented in Figure 12: the amplitudes of the density and pressures contrasts decay with conformal time while the gauge potentials still oscillate in their distorded interaction potential, insensitive to an expansion that is too close to a conformal transformation. The growing modes of the linear theory simply do not appear immediately due the scale invariance of YM fields. This means also that the YM cosmology can be much more inhomogeneous than it looks at the level of the observable quantities as the non observable ones (the metric components and the gauge potentials) do not decay in general, as we have seen before.
By the way, when space inhomogeneity is too strong, growing modes for the contrasts and the perturbations of the gauge potentials (see Figure 13 for illustration) appear, along with the diffusion already mentionned above. In fact, the details of this regime seem to be quite complex as it strongly depends on the interactions between the gauge components and on the field(s) initially perturbed. More work is under way to explore that regime and to characterize it in terms of the two dimensional abelian Higgs model of the theory.
To summarize the evolution of Figures 5,6,7,8,11, if one starts with some fluctuations in a universe dominated by YM fields, the density and pressures contrasts will enter an oscillating regime, firstly damped by the quasi-homogeneous expansion. Then, when the expansion becomes significantly inhomogeneous, the constrasts grow, at a time depending on the nature, the size and the amplitude of the initial perturbation. The transition to growing modes (and diffusion) seems indeed to occur when the perturbation of the metric coefficients δ m = (m − 1)R or δ l = (l − 1)R 2 χ 2 reach the order of magnitude of the scale factor (here a few percents, see Figure 5). Therefore, the transition to a non conformal expansion regime is to be linked directly with the emergence of the instabilities.
If the existence of growing modes is due here to the fact that the inhomogeneous expansion looks no longer like a conformal transformation of the metric and therefore yields to non trivial interactions with the inhomogeneous gravitational field, it would be interesting to see what happens to the growing modes with other gauges than the conformal one (N = R(t)) used here, for example, in the case of maximal slicing (cf. [21] and references therein). Indeed, this will lead to the study of a true gravitational collapse of YM fields instead of a shiny gravitational instability. Now, let us be more precise about how the parameters of the simulation affect the evolution of the observable quantities. A run is specified by the following parameters: ρ B 0 , the initial background density ; w, related to the coherence length λ P of the fluctuations ; ǫ, the relative amplitude of the perturbation ; σ 0 , the initial position of the Gal'tsov-Volkov particle in the interaction potential, and the nature of the perturbation (what fields were initially perturbed). The first two define the ratio between the coherence and Hubble lengths λ P L H . Here, we mostly focus on perturbations larger than the Hubble length: in Figure 11, λ P ≫ L H (λ P ≈ 3, L 0 H ≈ 0.09); in Figure  12, λ P >>> L H (λ P ≈ 2, L 0 H ≈ 0.009) and finally in Figure 13, λ P > L H (λ P ≈ 2, L 0 H ≈ 0.3). In the first case, we see that the density contrast has grown after half a period of σ ; in the second case and for the same period, we do not see yet any instability arising and finally, in Figure 13, the contrasts have grown quite significantly during half a period of σ. A more complete analysis shows that, all other parameters being equal, the large fluctuations λ P ≫ L 0 H evolves more slowly than the small λ P ≈ L 0 H . This can be intuitively deduced from the previous comments on Figures 11-13 and is clearer by examining Figure 14 where the evolutions of contrasts for two different coherence lengths are represented, all other parameters being equal. The evolution of the contrasts is qualitatively the same, although it is slower for large fluctuations (Figure 14 (iv), (v), (vi)). Please note the different duration and amplitudes of the two simulations. Therefore, the details of the evolution, e.g. the growing time, the period and/or the amplitudes of the contrasts, depend on the coherence length, the absolute amplitude of the perturbations (defined by for a certain field f ) which is related to ǫ and σ 0 , and the nature of the perturbation. However, the scheme damped-then-growing oscillations seems to be general. We even performed some computations at short wavelengths λ P ≪ L H but we did not recover the usual Jeans criterion for radiation (cf. [25]). We think that this criterion could be modified somehow due the locally anisotropic nature of YM pressures 7 (cf. [26] and references therein for some properties of locally anisotropic fluids). More work is under way to fix that point.
In Figure 11 (iv), we see that the non-diagonal component of the stressenergy tensor T 1 0 becomes more and more homogeneous with expansion. This was observed in all simulations and can be deduced from Eqs. (17 -33) where one can see that this quantity evolves as the inverse fourth power of the scale factor R(t).
For hyperbolic and closed cosmologies, we have performed some computations that lead essentially to analogous features as those presented here in the flat case, for fluctuations of the order of the Hubble length. G 3 (see Figure 15). The central bump in each constraint is related to the truncation errors when evaluating the spatial gradients, this can be improved by increasing the number of discretization points in the variable χ. The interesting information is the time-evolution of the constraints: we see that the violation increases almost linearly during the computation and accelerates at the end. This is typical of the very simple explicit scheme we used here. However, it is important to note that the constraints remain several orders of magnitude below the perturbations of the related fields thanks to the very strong Courant condition we used. Another crucial point is the boundary effect as the careful reader will have already noticed from previous graphs. This can partly be improved by refining the accuracy of the IVP and partly by using more localized fluctuations. Some compactification of the χ domain could help in verifying more exactly the boundary conditions. In order to produce a more accurate method, one should also focus on improving the stability by including implicit features in the scheme as well as a better control of boundary conditions.

Conclusion
By using a numerical method based on the hamiltonian formulation of the EYM system with SU(2) gauge group, we have checked that the primeval excitations of YM fields have been diluted by strong expansion in the very early universe, due to the scale invariance of YM equations. The gauge fields oscillate in a distorded potential well inherited from the homogeneous solution while the metric suffers an ever increasing departure from homogeneity. Meanwhile, the corresponding density and pressures contrasts undergo firstly damped oscillations.
However, the YM fields are gravitationally unstable in the long term. Indeed, the metric perturbations grow monotonically with time and, when the inhomogeneity of the metric and the gauge potentials cannot be neglected anymore, the last start diffusing non-linearly. This yields growing oscillation modes for the density and pressures contrasts, with a transition from damping depending on the interactions between the perturbed gauge fields, the nature, size and amplitude of the primeval fluctuations. For example, the larger the fluctuation w.r.t. the Hubble length, the slower it will grow. A fruitful comparison can be made with the gravitational instability of the inflaton which lead us to the preliminary conclusion that the growth rate of the YM fluctuations should depend also on the strength of the gauge coupling constant and should be more important for strong interaction.
Nevertheless, in the real universe, we can conjecture that the phase transitions that made YM fields short-ranged occur probably long before their inhomogeneities reach a sufficient size. Therefore such fields could be of minor importance for cosmology. There are, however, two main objections to that affirmation. First, a more physical and quantitative analysis, including precise considerations about the amplitude and size of the primeval excitations of the gauge fields, other components of matter, precise cosmological parameters as well as quantum considerations on the status of the gauge fields after the phase transitions should be performed before concluding. In second, considering a Born-Infeld type modification of the YM lagrangian, deriving from string theory, which breaks the conformal invariance can be of interest. In particular, an interesting work would be to study the gravitational instability of non-abelian Born-Infeld cosmology which was examined in details in [11].
Furthermore, another interesting point for cosmology is that the conformally invariant YM fields can hide for a while the deep inhomogeneity of space-time. Indeed, in a universe dominated by such fields, the observable quantities tend first to mimic a homogeneous universe while the metric is getting more and more inhomogeneous and the perturbations of the gauge fields oscillate. The question is to know which role such inhomogeneities could have played after the phase transitions that left the YM fields short-ranged as other types of matter fields feel differently the riddles of space-time hidden during the gauge dominating era than the conformally invariant gauge fields did.
For all these reasons, we think that our results are of interest at least for the investigation of mathematical properties of the pure EYM theory, if not for a tiny step in the understanding of the very early universe .           (i) ρ B − 1 (iv) T 1 0 (same parameters as in Figure 4).     Figure 4)