Structure of longitudinal chromomagnetic fields in high energy collisions

We compute expectation values of spatial Wilson loops in the forward light cone of high-energy collisions. We consider ensembles of gauge field configurations generated from a classical Gaussian effective action as well as solutions of high-energy renormalization group evolution with fixed and running coupling. The initial fields correspond to a color field condensate exhibiting domain-like structure over distance scales of order the saturation scale. At later times universal scaling emerges at large distances for all ensembles, with a nontrivial critical exponent. Finally, we compare the results for the Wilson loop to the two-point correlator of magnetic fields.


I. INTRODUCTION
Heavy ion collisions at high energies involve non-linear dynamics of strong QCD color fields [1]. These soft fields correspond to gluons with light-cone momentum fractions x 1, which can be described in the "Color Glass Condensate" (CGC) framework. Because of the high gluon occupation number the gluon field can be determined from the classical Yang-Mills equations with a static current on the light cone [2]. It consists of gluons with a transverse momentum on the order of the density of valence charges per unit transverse area, Q 2 s [3]. Parametrically, the saturation momentum scale Q s separates the regime of non-linear color field interactions from the perturbative (linear) regime. It is commonly defined using a two-point function of electric Wilson lines, the "dipole scattering amplitude" evaluated in the field of a single hadron or nucleus [4] as described below.
Before the collision the individual fields of projectile and target are two dimensional pure gauges; in light cone gauge, where m = 1, 2 labels the projectile and target, respectively. Here V m are light-like SU(N c ) Wilson lines, which correspond to the eikonal phase of a high energy projectile passing through the classical field shockwave [5,6]. The field in the forward light cone after the collision up to the formation of a thermalized plasma is commonly called the "glasma" [7]. Immediately after the collision longitudinal chromo-electric and magnetic fields E z , B z ∼ 1/g dominate [7,8]. They fluctuate according to the random local color charge densities of the valence sources. The magnitude of the color charge fluctuations is related to the saturation scale Q 2 s . The transverse gauge potential at proper time τ ≡ √ t 2 − z 2 → 0, is given by [9] A i = α i 1 + α i 2 .
Note that while the fields of the individual projectiles α i m are pure gauges, for a non-Abelian gauge theory A i is not. Hence, spatial Wilson loops evaluated in the field A i are not equal to 1. The field at later times is then obtained from the classical Yang-Mills equations of motion, which can be solved either analytically in an expansion in the field strength [9,10] or numerically on a lattice [11][12][13]. The Wilson loop, and the magnetic field correlator, provide an explicitly gauge-invariant method to study the nonperturbative dynamics of these fields, complementary to studies of the gluon spectrum [14].
Spatial Wilson loops at very early times τ have recently been studied numerically in Ref. [15], using the MV model [2] for the colliding color charge sheets. It was observed that the loops effectively satisfy area law scaling for radii > ∼ 1/Q s , up to a few times this scale. Furthermore, Ref. [16] found that two-point correlators of B z over distances 1/Q s correspond to two dimensional screened propagators with a magnetic screening mass a few times Q s . This indicates that the initial fields exhibit structure such that magnetic flux does not spread uniformly over the transverse plane (like in a Coulomb phase) but instead is concentrated in small domains.
The present paper extends this previous work as follows. We perform lattice measurements of spatial Wilson loops over a much broader range of radii to analyze their behavior at short (R 1/Q s ) and long (R 1/Q s ) distances. We also implement the so-called JIMWLK [3,17,18] high-energy functional renormalization group evolution which resums observables to all orders in α s log(1/x). High-energy evolution modifies the classical ensemble of gauge field configurations (4), (5) to account for nearly boost invariant quantum fluctuations at rapidities far from the sources. Finally, we also solve the Yang-Mills equations in the forward light cone to study the time evolution of magnetic flux loops.
The calculation of the initial conditions and the numer-arXiv:1401.4124v1 [hep-ph] 16 Jan 2014 ical solution of the classical boost-invariant 1 Yang-Mills fields in the initial stages of a heavy ion collision have been documented in the references given below, so here we will only describe them very briefly in Sec. II before moving on to show our results in Secs. III and IV.

II. LATTICE IMPLEMENTATION
We work on a two dimensional square lattice of N 2 ⊥ points with periodic boundary conditions and consider color sources that fill the whole transverse plane. The lattice spacing is denoted as a, thus the area of the lattice in physical units is The calculations are performed for N c = 3 colors. In this work we only consider symmetric collisions, where the color charges of both colliding nuclei are taken from the same probability distribution.
In this work we compare three different initial conditions for the classical Yang-Mills equations: the classical MV model (5) as well as fixed and running coupling JIMWLK evolution. We define the saturation scale Q s (Y ) at rapidity Y through the expectation value of the dipole operator as Throughout this paper we shall use Q s defined in this way from the light-like Wilson lines V (x T ) in the fundamental representation. The saturation scale is the only scale in the problem and we attempt to construct the various initial conditions in such a way that the value of Q s a is similar, to ensure a similar dependence on discretization effects. In the MV model the Wilson lines are obtained from a classical color charge density ρ as where P denotes path-ordering in x − . The color charge density is a random variable with a local Gaussian probability distribution The total color charge dx − µ 2 (x − ) ∼ Q 2 s is proportional to the thickness of a given nucleus.
In the numerical calculation the MV model initial conditions have been constructed as described in Ref. [13], discretizing the longitudinal coordinate Y in N y = 100 1 The YM equations are solved in terms of the coordinates τ = √ t 2 − z 2 , η = 1 2 ln t+z t−z and x T ; hence steps. For the calculations using the MV model directly for the initial conditions (1), (2) we have performed simulations on lattices of two different sizes: N ⊥ = 1024, with the MV model color charge parameter g 2 µL = 156 which translates into Q s a = 0.119; and with N ⊥ = 2048, using g 2 µL = 550, which results in Q s a = 0.172. The MV model also provides the configurations used as the initial condition for quantum evolution in rapidity via the JIMWLK renormalization group equation, starting at Y = log x 0 /x = 0. Performing a step ∆Y in rapidity opens phase space for radiation of additional gluons which modify the classical action (4), (5). This process can be expressed as a "random walk" in the space of light-like Wilson lines V (x T ) [18][19][20]: where the Gaussian white noise η i = η i a t a satisfies η a i (x T ) = 0 and, for fixed coupling, Here the equation is written in the left-right symmetric form introduced in [20,21]. The fixed coupling JIMWLK equation is solved using the numerical method developed in [19,20,22]. For the smaller lattice size N ⊥ = 1024 we start with the MV model with g 2 µL = 31 and without a mass regulator, which corresponds to an initial Q s a = 0.0218. After ∆y = 1.68/α s units of evolution in rapidity 2 this results in Q s a = 0.145. For a N ⊥ = 2048-lattice we again start with g 2 µL = 31, corresponding to Q s a = 0.0107, and after ∆y = 1.8/α s units of evolution end up with Q s a = 0.141.
For running coupling the evolution is significantly slower. We use the running coupling prescription introduced in [20], where the scale of the coupling is taken as the momentum conjugate to the distance in the noise correlator in Eq. (7). For the smaller N ⊥ = 1024 lattice we again start with g 2 µL = 31, i.e. Q s a = 0.0218 and evolve for ∆Y = 10 units in rapidity, arriving at Q s a = 0.118. For the larger N ⊥ = 2048 lattice we test a configuration that is farther from the IR cutoff, starting the JIMWLK evolution with g 2 µL = 102.4, i.e. Q s a = 0.0423 and evolve for ∆Y = 10 units in rapidity, arriving at Q s a = 0.172. In the rc-JIMWLK simulations the QCD scale is taken as Λ QCD a = 0.00293 and the coupling is frozen to a value α 0 = 0.76 in the infrared below 2.5Λ QCD .
As already mentioned above, RG evolution in rapidity resums quantum corrections to the fields α µ m of the individual charge sheets to all orders in α s log 1/x, with leading logarithmic accuracy. In other words, the effective action at Y is modified from that at Y = 0, written in Eq. (5).
Once an ensemble of Wilson lines V (x T ) at a rapidity Y is constructed, separately for both projectile and target, these configurations define α i 1 and α i 2 in light-cone gauge as written in Eq. (1); the initial field A i of produced soft gluons at proper time τ = +0 corresponds to their sum, Eq. (2). The evolution to τ > 0 follows from the real-time Hamiltonian evolution described in Ref. [11]. This has been used in many classical field calculations, e.g. in Refs. [12], or more recently for the first study of the effects of JIMWLK-evolution on the gluon spectrum [14], and in the IP-glasma model for the initial conditions for hydrodynamics [23]. On the N ⊥ = 2048 lattices we evolve the fields up to Q s τ = 5 and on the smaller N ⊥ = 1024 ones to Q s τ = 10.

III. WILSON LOOP
In the continuum the spatial (magnetic) Wilson loop is defined as the trace of a path ordered exponential of the gauge field around a closed path of area A in the transverse plane: On the lattice this is easily discreticed as the product of link matrices around a square of area A. For N c ≥ 3 colors any particular Wilson loop is complex but the ensemble average is real.
We have measured the expectation value of the Wilson loop in the glasma field, with different initial conditions and at different times Q s τ . The results of the calculation are shown in Fig. 1. As expected, the magnetic flux through a loop generically increases with its area. Focusing first on the curves corresponding to the initial time τ = 0 we observe that the resummation of quantum fluctuations (JIMWLK evolution) increases the flux through small loops of area AQ 2 s < 1. This can be understood intuitively as due to emission of additional virtual soft gluons in the pure gauge fields of the colliding charge sheets. On the other hand, the flux through large loops, AQ 2 s > ∼ 2, decreases. This indicates uncorrelated fluctuations of magnetic flux over such areas and is consistent with the suggestion that the flux is "bundled" in domains with a typical area ∼ 1/Q 2 s [15]. Accordingly, loops of area ∼ 1.5Q 2 s are invariant under high-energy evolution. Moving on to finite times we see that the flux through loops of fixed area decreases with τ . This is, of course, a consequence of the decreasing field strength in an expanding metric. For small loops the ordering corresponding to the different initial conditions (MV, rc-JIMWLK, fc-JIMWLK) is preserved even at later times. However, for large loops one observes a striking "universality" emerging at Q s τ ∼ 5 as the curves for all initial conditions fall on top of each other.
The data from Fig. 1 shows an approximately linear dependence of ln(− ln W )) on ln(AQ 2 s ), with different slopes in the regime of small AQ 2 s 1 vs. large AQ 2 s 1. Based on this observation we fit the data to with separate parametrizations for the IR and UV regimes: IR : e 0.5 < AQ 2 s < e 5 (10) UV : e −3.5 < AQ 2 s < e −0.5 .
In addition to limiting the fits to the quoted ranges we also restrict them to the region where W > 0.01 and where W > 0 by a 5σ statistical error margin; beyond these limits the data exhibits too large fluctuations for a meaningful fit. Figures 2 and 3 show the time dependence of the exponents γ in the IR and UV regions. The "string tension" σ naturally decreases as ∼ 1/τ because of the longitudinal expansion of the glasma, which leads to B z ∼ 1/ √ τ . We therefore show, in Figs. 4 and 5, the time dependence of the combination τ σ/Q s , where this leading effect is scaled out. The values of σ/Q 2 s for Q s τ = 0 are given in the captions 3 .
3 For the MV model, at τ = 0 we find σ/Q 2 s = 0.44 in the IR region which is about four times larger than the value reported in ref. [15]. Our present results refer to Nc = 3 colors while ref. [15] considered Nc = 2; also, our current definition of Qs via Eq. (3) leads to smaller values for this quantity than the definition used in [15]. Finally, σ is extracted from fits over a somewhat different range of areas.  The results in the "UV"-regime probed by Wilson loops of small area are shown in Figs. 2 and 4. They are easily understood from the differences in the initial condition. The gluon spectrum in the MV model falls steeply as a function of k T , leading to a steep dependence of short-distance correlators on the distance. Our result for the UV exponent in Fig. 2 is close to the A 2 -scaling obtained analytically in a weak field expansion [24]. The difference is probably due to a combination of logarithmic corrections and lattice UV cutoff effects. For the JIMWLK ensembles the gluon spectrum is much harder [14], especially for fixed coupling. This manifests itself in smaller values of both γ and σ. In addition, the UV exponents are remarkably time independent at Q s τ > 1: this is consistent with the expectation that at such time the UV modes can be viewed as noninteracting gluons whose spectrum is close to the expectation from a perturbative k T -factorized calculation [25].
The behavior in the IR regime (Figs. 3, 5) probed by large Wilson loops points to a very different picture. At τ = 0 the exponents γ and, to a lesser extent, the values of σ depend very much on the initial conditions. As already alluded to above, the scaling exponents γ IR < 1 obtained for the JIMWLK fields indicate that quantum emissions increase magnetic flux fluctuations at the scale ∼ 1/Q s , much smaller than the area of the loop. It is interesting to note that for the rather strong fixed-coupling evolution the initial scaling exponent is not too far above γ IR = 1/2 corresponding to perimeter scaling.
At times Q s τ 3, however, one observes a remarkable universality in the IR as the curves corresponding to different initial conditions collapse onto a single curve in Fig. 1. The string tensions in Fig. 5 are within 10% of each other at late Q s τ , and the exponents γ in Fig. 3 are very close to each other, with values around γ IR ≈ 1.2 . . . 1.3. The exponent gradually decreases with τ , potentially approaching the area law γ = 1 at late times. The initial evolution points at a rapid rearrangement of "magnetic hot spots" to some universal field configurations at later time, Q s τ 3.
We stress that the universal behavior of large magnetic loops, characterized by a nontrivial power-law dependence on the loop area, sets in at rather early time scales of a few times 1/Q s , independent of initial conditions. Actual area law scaling γ = 1 is approached only later. This behavior mirrors a similar universality between MV and JIMWLK results seen in the IR part of the gluon spectrum (determined from correlators of gauge fixed fields) in Ref. [14]. Since the structure of the fields does not seem to depend on the initial conditions, we infer that this universality in due to stong interactions in the glasma phase. This universal behavior of the Wilson loop for different initial conditions at Q s τ 3 and AQ 2 s 1 is the main result of this paper.

IV. MAGNETIC FIELD CORRELATOR
In this section we analyze gauge-invariant two-point magnetic field correlators of the form 4 The points x T and y T are separated in the x or y direction by a distance r = |x T − y T |, and the Wilson line U x T →y T is the ordered product of links along the straight line separating these points.
The magnetic field B z = t a B a z on the lattice is defined as the traceless antihermitian part of the plaquette as where the transverse plaquette is Here U i (x T ) denotes the link matrix in the i-direction based at x T and i T , j T are unit vectors. The resulting magnetic field correlator r C B (r) at time τ = 0 and τ = 10/Q s is plotted in Fig. 6. We have multiplied by r to better expose the behavior around r ∼ 1/Q s . At the initial time there is a significant anticorrelation at intermediate distances. It shows the domain structure of the field such that B z is likely to flip sign 5 (or direction) over distances of order 1/Q s . At later time the correlator exhibits a more regular behavior: the anti-correlation has disappeared and the correlator shifts to larger Q s r.
The measurements at larger Q s r are dominated by fluctuations. Given the clear scaling behavior of the Wilson loop one might expect to see a similar phenomenon for the magnetic field correlator. A very naive scaling argument would assume that if C(r) ∼ r −α then the Wilson loop should scale as The area integrals in the first line extend over |x T |, |y T | < R. Q s τ = 10 would give C B (r) ∼ r −1.55 . On a logarithmic scale C B (r) does indeed qualitatively resemble such behavior as shown in Fig. 7. However, this kind of scaling is less conclusive than for the Wilson loop (see also appendix); this could be an indication for the presence of higher cumulants in the expansion of the spatial Wilson loop [26].

V. SUMMARY
In this paper we have provided some insight into the fields produced initially in a high-energy collision of dense color charge sheets. We have focused, in particular, on the structure of the longitudinal magnetic field B z .
We consider both purely classical as well as JIMWLK RG evolved gauge field ensembles on which we measure expectation values of spatial Wilson loops and two-point correlation functions of B z . These show that the initial fields exhibit domain-like structure over distance scales of the order of the inverse saturation scale 1/Q s . Classical YM evolution to later times leads to universal scaling, for all ensembles, of the magnetic loop with area, with a nontrivial critical exponent. Also, the anti-correlation of B z (x T ) over distances ∼ 1/Q s disappears, which we interpret as rearrangement of the magnetic field domains.
The emergence of a color field condensate in highenergy collisions of dense hadrons or nuclei is a very interesting phenomenon, and its dynamics remains to be understood in more detail. In closing we only draw at-tention to recent arguments that the presence of such a condensate might have important implications for the process of (pre-) thermalization in high multiplicity collisions [27].

Appendix A: Relation between the Wilson loop and magnetic field correlator
In an abelian theory there is a simple relation between the Wilson loop and the magnetic field due to Stokes' theorem: If we assume that in the nonabelian case the magnetic field in each color channel a is independent, and that it consists of uncorrelated domains which are much smaller than the area A and distributed as Gaussian random variables, we obtain the following estimate for the Wilson loop: In Fig. 8 we compare the result of a numerical integration of the r.h.s. of Eq. (A2) using the measured magnetic field correlator, to the direct measurement of the Wilson loop. It can be seen that the two are in a relatively good agreement. This consistency check supports the interpretation of B z as independent field domains of area ∼ 1/Q 2 s .