Non-Abelian chiral instabilities at high temperature on the lattice

We report on an exploratory lattice study on the phenomenon of chiral instabilities in non-Abelian gauge theories at high temperature. It is based on a recently constructed anomalous Langevin-type effective theory of classical soft gauge fields in the presence of a chiral number density $n_5=n_{\rm R}-n_{\rm L}$. Evaluated in thermal equilibrium using classical lattice techniques it reveals that the fluctuating soft fields indeed exhibit a rapid energy increase at early times and we observe a clear dependence of the diffusion rate of topological charge (sphaleron rate) on the the initial $n_5$, relevant in both early universe baryogenesis and relativistic heavy-ion collisions. The topological charge furthermore shows a drift among distinct vacuum sectors, roughly proportional to the initial $n_5$ and in turn the chiral imbalance is monotonously reduced as required by helicity conservation.


Introduction
Elucidating the properties of chiral plasmas is central to deepen our understanding of the physics of the electroweak era of the early Universe [1,2], the quark-gluon plasma created in heavy-ion collision experiments [3,4], as well as electron and neutrino media in supernovae [5,6], to name a few. What distinguishes a chiral plasma from conventional plasmas is the presence of anomalous transport, such as the chiral magnetic effect (CME) [4,[7][8][9]. In this particular phenomenon a current parallel to magnetic fields is induced if a chiral imbalance, characterized by the chiral density n 5 ≡ n R − n L , is present. When the participating gauge fields fluctuate, the CME leads to a new type of plasma instability christened the chiral plasma instability (CPI), or simply, the chiral instability [10,11]. The time evolution of phenomenologically relevant chiral plasmas exhibiting such an instability (both electroweak or quark-gluon) has not yet been fully understood, partly due to the nonlinear nature of non-abelian gauge fields. (See Refs. [12][13][14] for recent attempts on Abelian gauge theories with a chiral imbalance.) Take as example the the yet-to-be-settled question of whether the CME can be observed in relativistic heavy-ion collision at RHIC and LHC via event-by-event charge separation measurements [15][16][17]. Its solution will require a thorough quantitative understanding of anomalous transport including its physical time scales and thus the effect of the chiral instability. In this paper, we report on a first study of the nonlinear evolution and the fate of the chiral instability for high-temperature non-Abelian chiral plasmas using non-perturbative real-time lattice methods.
In a heavy-ion collision both color and electromagnetic fields are at play. The two incident nuclei traveling along the light cone can be described in the theory of the Color Glass Condensate [18,19] as highly occupied gluon fields at small Bjorken x and a collection of classical color sources at large x [20,21]. At the collision instance the interplay between the two sectors leads to the emergence of strong boost-invariant color electric and color magnetic fields, parallel to each other, spanning between the receding collision fragments [22]. The fact that both fields lie parallel to each other corresponds to the presence of non-vanishing chiral charge. These initially strong and classical (or nonfluctuating) fields subsequently decay either via classical Yang-Mills evolution or via particle production facilitated by the Schwinger mechanism. There have been indications from model studies [23,24] that quark production might happen quite rapidly, such that chiral fermions are present already at very early times: after a short but finite time the collision region can fragment in many subdomains of finite topological charge possessing a corresponding chiral imbalance as would be implied by the Atiyah-Singer index theorem in a Euclidean setting. 1 The fact that we cannot cite definitive timescales in fm for these processes is a reminder of the need for more quantitative numerical studies elucidating these dynamics from first principles at early times, currently underway in several groups [26][27][28].
If the incident nuclei are viewed as classical point-sources of electromagnetic charge [29,30], Maxwell's equations in the form of the Líenard-Wiechert potentials tell us that they too generate electric and magnetic fields that lie parallel to each other. With field strengths of the order of the pion mass squared they may be capable of influencing the subnuclear composition of the collision bulk region. While the recession of the nuclei quickly diminishes such magnetic fields, the presence of quark degrees of freedom early on, i.e., their contribution to a non-equilibrium electric conductivity (for a state-of-theart equilibrium estimate see Refs. [31,32]) can potentially elongate the lifetime of these U(1) magnetic fields. The interplay between the U(1) magnetic fields and a possible chiral imbalance is a possible candidate to produce charge separation in relativistic heavy-ion collisions [33,34].
We are interested in the system at a point in time shortly after the initial collision, at which the external classical fields generated by the collision fragments have become negligible. The collision region is assumed to host several distinct regions of CP-odd matter corresponding to a finite chiral density n 5 and an accompanying finite Chern-Simons number, or topological charge. Overall the bulk will possess vanishing chiral charge. An important aspect of this idealized setting is that helicity is at this point a conserved quantity, and so changes in the chiral imbalance have to be necessarily accompanied by opposite changes in the topological charge of the dynamical gauge fields.
We now take one of the CP-odd domains and zoom into its interior, being interested in the physics of fluctuations of the soft gauge fields, which operate on length scales much smaller than the extend of the CP-odd domain. Hence we will take n 5 in the following to be spatially homogeneous from the point of view of the gauge sector. The quantum back-reaction of the chiral imbalance onto the gauge fields will then lead to some of the gauge field modes to become unstable and to exhibit rapid growth over time, the chiral instability. How these modes evolve and eventually saturate and how the presence of the chiral imbalance will change the transitions between topological sectors within the the CP-odd domain is what we wish to shed light on in the following sections.
We are able to investigate these phenomena related, to the quantum back-reaction of the chiral fermions onto the gauge field sector by deploying a recently constructed Langevin-type effective theory for soft gauge fields that captures the chiral instability [11]. Here the label soft refers to the scale g 2 T , where g is the coupling constant of the gauge theory and T is the temperature. In the following, we will use the label hard for the scale T and semi-hard for the scale gT . Using classical statistical simulation techniques we numerically determine the diffusion constant for Chern-Simons number [see Eq. (2.12) for its definition], or the so called sphaleron rate. As will be shown in Fig. 12 below, we find that the presence of n 5 together with the chiral instability quantitatively modifies the sphaleron rate with potential phenomenological consequences for relativistic heavy-ion collisions but also for the question of baryogenesis in the early Universe. 2 The organization of the paper is as follows. In Sec. 2 we review the continuum formulation of the chiral Langevin theory constructed in Ref. [11]. Section 3 discusses our implementation of the chiral Langevin theory in the context of classical statistical lattice simulations in thermal equilibrium and the technical challenge associated with the determination of topological charge. Our main numerical results on the chiral instability is subsequently presented in Sec. 4., while Sec. 5 is devoted to a discussion of these results and concludes the paper.

Chiral Langevin theory in continuum
We start here with a brief summary of the Langevin-type effective theory for soft gauge fields (with momentum k ∼ g 2 T ) at high temperature in the presence of a chiral imbalance [11]. It is assumed that the system temperature is sufficiently large compared to the dynamical scale of the non-Abelian gauge theory so that g 1. The first relevant scale below the hard scale T is the semi-hard electric scale k ∼ m D ∼ gT at which electric fields still propagate, but the corresponding gluons receive a thermal mass m D , the Debye mass. At even lower momenta, k ∼ g 2 T the electric fields have become fully screened and do not propagate any longer. At this soft or magnetic scale the color magnetic fields remain as dynamical degree of freedom possessing a magnetic screening mass m M ∼ g 2 T , which is generated via fully nonperturbative dynamics at this scale.
Note first that the soft degrees of freedom of can be regarded as classical, because these modes with the momentum k T are highly occupied in equilibrium as Note also that the soft gauge fields with A ∼ gT and k ∼ g 2 T interact strongly among themselves via the classical Yang-Mills equation, because the kinetic term (k) and the interaction term (gA) in the covariant derivatives are of the same order in g. In addition to their self interaction, the key dynamics of the soft gauge fields arises from the color currents j hard of hard particles with the momentum p ∼ T . The non-Abelian version of Ampere's law reads where D µ is the covariant derivative. Here and below, E and B denote the color electromagnetic fields and we will often suppress color indices for simplicity.

Langevin-type effective theory without chiral imbalance
To make the paper self-contained, we first consider the Langevin-type effective theory for soft gauge fields without a chiral imbalance [44,45,[54][55][56]. The color current manifests itself as an effective degree of freedom, which in the leading log approximation can be expressed locally in terms of the soft gauge fields, where γ is the damping rate for hard gauge bosons and the coefficient σ c is called the color conductivity in analogy with the Ohmic law for the electric current (even though color is not a conserved number). The normalization of the right-hand side of Eq. (2.5) is chosen such that the noise ζ satisfies the fluctuation dissipation theorem.
In the presence of soft gauge fields, the color current relaxes towards an Ohmic current in Eq. (2.3) due to semi-hard scatterings with the momentum transfer q ∼ gT among the hard particles. The reason why semi-hard scatterings dictate the relaxation of the color current is that the effect of the semi-hard scattering on the momentum of hard particles is small as p ∼ T q, while that on their color is significant. Therefore the non-Abelian nature of the gauge fields is essential to obtain the color conductivity in Eq. (2.4), which allows us to approximate the color current in the local form as Eq. (2.3) in the leading log approximation. 3 To proceed let us take the gauge A 0 = 0. From Eqs. (2.2) and (2.3) we can determine the time scale of the soft gauge fields, τ soft . We assume that τ soft 1/σ c , which will be justified in a moment. Then we can ignore the D t E term and obtain which validates the assumption above. By writing down a three-dimensional effective Hamiltonian of the color magnetic fields as  This Langevin equation (2.9) drives the soft gauge fields towards the equilibrium distribution: (2.10) The soft gauge fields give the dominant contribution to the Chern-Simons number diffusion rate Γ: where f abc is the structure constant. The N CS -changing processes require the condition, g 2 ABR 3 ∼ g 2 B 2 R 4 1, where R ∼ 1/k is the wave length of the mode. The energy for such a mode is E A ∼ B 2 R 3 1/(g 2 R), which should satisfy E A T so as not to be suppressed in the thermal environment. This gives the condition 1/k ∼ R 1/g 2 T , but gauge fields with k g 2 T are suppressed due to the magnetic screening mass m M ∼ g 2 T . Therefore, the soft gauge fields k ∼ g 2 T contribute dominantly to the Chern-Simons number diffusion rate. Since we know that the soft gauge fields evolve at time scale τ soft ∼ 1/[g 4 T ln(1/g)] [see Eq. (2.7)], the Chern-Simons number diffusion rate is parametrically estimated to be The scales of the soft gauge fields and their power counting are summarized in Tab. 2.1.

Langevin-type effective theory with chiral imbalance
So far, we have not considered the nature of the hard modes in the plasma, which couple to the soft gauge fields through the current j hard = σ c E + ζ. In fact the hard modes with the momentum p ∼ T consist of gauge fields and fermions that interact with the soft gauge fields. Without the imbalance of chiral fermions, they would affect the Debye screening mass m D only via the number of colors and flavors. However, in the presence of chiral fermions, the N CS -changing processes in the soft sector induce changes of their chirality through the anomaly relation, (2.14) Therefore, the Langevin equation needs to be amended to take into account nonzero chiral charge [11]. The evolution of the chiral imbalance N 5 = N R − N L itself is driven by the change in the topological chargė where we denote the chiral chemical potential by µ 5 = µ 5 (t). Equation (2.15) means that total helicity N 5 + N CS is conserved as a function of time. The back reaction of the chiral fermions on the gauge sector is implemented via the non-Abelian version of a chiral magnetic current proportional to the color magnetic field B, In the form of Eq. (2.9), the latter modification is compactly expressed by This derivation of the modified Langevin theory utilizes a power counting of g assuming µ 5 ∼ T in the parametric sense. Note also that all the terms of j hard are of the same order g 5 T 3 . This modified effective theory [11] possesses the following two new properties which the original one [44,45,[54][55][56] does not have: 1. Chern-Simons number diffusion receives back reaction from the fermion sector.
2. Soft gauge fields possess an unstable mode in the presence of finite µ 5 as argued in Ref. [10].
Let us consider the time scales of the associated phenomena. The order of the righthand side of Eq. (2.15) is g 10 T 4 ln(1/g), and so the time scale for µ 5 ∼ T is thus estimated to be The time scale of the instability is obtained from Eq. (2.18) as 4 The unstable mode carries the Chern-Simons number so that it reduces the chiral imbalance µ 5 . Due to the separation of the time scales, τ inst τ µ 5 , the instability develops and saturates due to nonlinear interactions among the soft gauge fields during the time in which µ 5 only changes very little.
At the longer time scales τ µ 5 , the fermion chiral charge decreases towards its equilibrium value. The questions here are • What is the fate of the instability? How do the nonlinear interactions saturate the instability?
• What is the fate of the fermion chiral imbalance? In particular, what is the fraction of N 5 and N CS eventually?
In order to understand these questions which are highly nonperturbative and nonlinear, we now turn to numerical lattice simulations in the following.

Lattice discretized effective theory and numerical setup
Let us discuss in detail the numerical implementation of the real-time evolution of the soft gauge fields in the background of chiral fermions. The qualitatively new aspect of this work is the inclusion of anomalous affects related to the quantum back-reaction of a dynamically evolving chiral imbalance onto the gauge sector, described effectively by the term proportional to the B field in Eqs. (2.2), (2.15), and (2.17). Note that we are solving the full Hamiltonian equations of motions and do not neglect the D t E term. While implementing the Hamilton dynamics of the gauge fields is straightforward, the treatment of topological charge on the lattice constitutes a challenge. Since the lattice definition of the Chern Simons number suffers from the fact that it cannot be related to a total derivative, it is known that spurious fluctuations on the level of the lattice spacing contaminate its numerical evaluation (see, e.g., Ref. [57]). In order to robustly identify changes in the topology of our system we combine several techniques, developed in both classical statistical Yang-Mills simulations, as well as Euclidean lattice QCD over the last two decades.

Classical Yang-Mills on the lattice
Let us start with classical Yang-Mills theory on the lattice in the absence of fermionic degrees of freedom. A common approach to the derivation of its real-time Hamilton dynamics is to start from Wilson's formulation [58,59] in terms of an action on a finite hypercubic space-time grid with lattice spacings a s and a t and volume N s × N t . The gauge degrees of freedom are link variables U µ (x, t) that are related to the gauge fields A a µ (x, t) in the adjoint representation via where λ a represents either the Pauli matrices in the case of SU(2) (a = 0, 1, 2) or the Gell-Mann matrices in the case of SU(3) (a = 0, . . . , 7). The appropriately rescaled quantities T a denote the generators of the gauge group. We use a shorthand notation for the lattice spacing a µ=0 = a t and a µ=i = a s , where summation over indices is not implied. From the gauge links one can construct the smallest gauge invariant object, the plaquette U µν , which is related to the field strength tensor F µν as As is common practice, we consider only the simplest formulation of the Wilson action to derive the equations of motion. In order to arrive at an equation of motion of Hamiltonian form, we take the established route of fixing to temporal gauge A 0 (x, t) = 0, i.e., U 0 (x, t) = 1 after evaluating the stationarity condition From the functional derivative with respect to the spatial fields and using a historic choice of normalization E a k = F a 0k ga s /(2 √ N c ) we arrive at the Hamilton-like classical equation of motion, where the second line follows from a Taylor expansion of the temporal plaquette. The staples S j and S j are defined such, that when multiplied with the link U j from the right or left respectively, lead to the closed backward and forward plaquette. The particular form of these equations invites the use of a second order accurate O(dt 2 ) leap-frog scheme for their implementation. The electric fields are first evolved one half time step via a simple Euler update and then E and U are subsequently updated based on their shifted counterpart. While the time derivative for E is naively discretized the update of the links is implemented exactly via matrix exponentiation. For SU(2) this is achieved via the handy relation, while for SU (3) we would resort to the fast and robust matrix diagonalization of Ref. [60] based on Cardano's analytic formula. Classical Yang-Mills fields represent a constrained system, a fact reflected in the functional derivative in the A 0 direction leading to a relation between fields that does not involve a partial derivative in time, the so called Gauss constraint, If we start from an initial field configuration that respects the Gauss constraint, the application of Eq. (3.7) will preserve this property. Only the accumulation of rounding

Leap-frog evolve U(x,t) E(x,t) to equilibrate discard E and repeat until T 00 stabilizes
Euler t=1/2 step Themal Initial cond. errors gradually weakens the constraint, which can be conveniently monitored via the global quantity b,x Tr T b G(x, t) .
Note that while we have started from spatial and temporal links on a hypercubic spacetime grid, we have arrived at a formulation of electric fields and spatial links on individual three-dimensional time-slices, which evolve via (continuous) time derivatives in Eq. (3.7). The Hamiltonian corresponding to these classical equations of motion reads (3.10) In this study we wish to investigate the physics in a thermal equilibrium setting for which we note that the Boltzmann weight can be written in natural units as The properties of a thermal configuration are thus solely determined by the lattice parameter β L = 2Nc g 2 asT . In practice simulations are carried out at a fixed β L using a s = g = 1.
To generate thermal initial configurations we follow the well established strategy of Refs. [37,41,61] with the implementation of Refs. [62,63], which is described in detail in the following and schematically depicted in Fig. 1. It relies on the fact that the electric fields enter quadratically in the Hamiltonian.
1. Preconditioning of spatial links: Using the Wilson plaquette part of the Hamiltonian a standard three-dimensional multihit Metropolis Monte-Carlo implemen-tation (N hits = 10 × 400) brings the spatial links from either hot-or coldstart towards an ensemble member of the Markov chain representing P . This is not yet the correct thermal distribution but the links are already closer to the correct ones than, e.g., a randomly chosen configuration.
2. The electric field distribution in thermal equilibrium does not explicitly depend on the gauge fields and x,j,a E a j (x)E a j (x) and hence we may start out by drawing each field entry separately from a Gaussian distribution (3.12) 3. The electric fields drawn such contain both physical and unphysical modes, as only the former obey the Gauss law. Hence we need to project out those fields that are perpendicular to the constrained directions. A particularly clever way to do so is via the iterative prescription of Ref. [64], where one computes for each E k (x) the Gauss constraint at x and x +k and carries out the following update The choice ofκ = 0.12 has been found to lead to fast and stable convergence, allowing us to enforce the Gauss constraint down to x G(x) = 10 −15 .
4. After bringing the projected electric fields forward by a t /2 via a naive Euler step, they are evolved together with the gauge links according to the Hamilton dynamics of Eq. (3.7) to mutually equilibrate for (N eq. steps = 300). After the evolution phase, the current electric field is discarded and one repeats from step 2 by drawing another round of Gaussian E a k . This cycle is itself repeated N cycl = 40 times and at the same time the approach of the energy T 00 = F 0µ F 0µ + F µν F µν towards a time independent value is monitored.
For the steps involving pseudo-random numbers we deploy Lüscher's ranlux algorithm [66], as it is proven to yield statistically independent chains if started from differing seeds.
Naive classical thermal equilibrium is ill-defined in the continuum due to the Rayleigh-Jeans divergence, 5 nevertheless on the lattice it is regulated by the cutoff π/a s . Since in such a scenario all modes up to the cutoff are occupied we have to expect that lattice artifacts do significantly influence the simulated physics. The physics of the semi-hard scale at m D ∼ gT indeed suffers from lattice artifacts, as can be seen from the explicit a s dependence of the lattice regularized Debye mass parameter [38,[67][68][69] computed in lattice perturbation theory It is the physics at scales ∼ g 2 T and below that is considered to be reproduced quantitatively within the classical lattice theory. The physics of topological transitions is among the phenomena that can be studied in this approach.

Topological charge on the lattice
We have seen in Sec. 2 that the gauge field topology is intimately related to the dynamics of the chiral fermions via the anomaly and vice versa. Hence we need robust means to identify how the topology of the gauge field changes over time. In this study we use the naive continuum integral formula, where the F µν is defined from a clover-type approximation [70] summing over the field strength components extracted from four neighboring plaquettes Instead of the matrix logarithm one may also use the simple relation Since we need also access to the temporal plaquettes at this point we will keep a copy of the fields at the previous time step in memory even though they are not needed to evolve Eq. (3.7). This however means that only the two backward plaquettes are used to define F 0k = −F k0 . While in this study we deploy the standard clover approximation to the field strength tensor, in principle one can improve on the definition [71], which leads to an expression for E · B to which lattice spacing artifacts contribute only at O(a 4 s ). As has been discussed in the literature in great detail (see, e.g., Ref. [74]), the naive evaluation of Eq. (3.15) on a lattice suffers from the fact that it is not related to a total derivative. It hence becomes susceptible to fluctuations at small length scales that are unrelated to the topology. Indeed it was shown that modifying a single link [57] can change the lattice discretized value of dN CS dt . In other words, difference in topological charge between two configurations, depends in general on the path in field space (here along real-time) one takes to connect the two. In a finite temperature setting we should note that N CS does not necessarily has to take on integer values as it does in vacuum, as thermal fluctuations allow it to explore the topological sector it currently resides in. In the presence of these ambiguities, the only truly well-defined piece of information is what topological vacuum sector the current field configuration belongs to. It is this integer value we set out to distill in the following.
Note that focusing on the topological sector alone introduces other conceptual difficulties. It has been argued [40] that gauge fields can spend significant time in the vicinity of a sphaleron configuration at half-integer topological charge, moving slightly from one sector to the neighboring one and back over timescales of 1/(g 2 T ). If we naively counted the number of transitions between sectors this would then overestimate how topology actually diffuses over time. In our case we use the formula for the sphaleron rate Γ, where the rapid changes around a sphaleron configuration will average out in the late time limit. In order to numerically extract the change between topological sectors another possible path might lie in devising a Minkowski space implementation of the overlap operator and to measure the number of exact chiral zero modes, which at least in the Euclidean domain are related to the change in topology by the Atiyah-Singer index theorem. Another strategy laid out by Moore and Turok is the slave-field method [42], which explicitly constructs a representation of the topology of the gauge group.
A more practical way, which has been deployed widely in classical-statistical lattice simulations and which under the moniker of Wilson flow [65] has recently captured also the interest of the Euclidean-time lattice QCD community is dissipative cooling. It acts on the spatial links of the theory and loosely speaking damps away their high frequency modes. Namely, a copy of the current-time links U (x) is evolved in a fictitious cooling time τ according to In the case of four-dimensional Euclidean lattice QCD τ represents the flow-time. Since this procedure is formulated in a gauge invariant fashion, we know that it cannot simply modify field modes outside of a certain momentum shell but also has to influence the low lying modes in a certain fashion. Indeed, in the context of Euclidean Yang-Mills theory it has been shown that cooling with the naive Wilson plaquette Hamiltonian reduces the size of any instanton currently present on the configuration. Eventually it makes the topological object small enough that it cannot be resolved on the lattice and its contribution to N CS will vanish. Namely, with the naive plaquette Hamiltonian one can easily overcool a configuration and hence miss transitions. In response a number of studies were devoted to constructing classical actions that lead also to growth of instantons, so as to balance the overall behavior. In this study we choose to use the improved Wilson action x,µ,ν with = 0 introduced in Ref. [70] and discussed in, e.g., Refs. [72,73]. In practice it gives a much smoother behavior during cooling, i.e., the step size in cooling time can be chosen much coarser than in the case of the naive single plaquette. While it has been suggested [73] to monitor the deviation of each individual link during cooling to avoid overcooling we do not implement this strategy here.
To be able to obtain consistently cooled configurations, we furthermore deploy the fourth order Runge-Kutta update proposed in Ref. [65] to implement the cooling of the form ∂ τ U (τ ) = Z(U (τ ))U (τ ), given by As was proposed first in Ref. [74] one might consider an additional coarsening step at intermediate cooling time, in which one replaces two neighboring links with the product of these two links, effectively reducing the number of points on the lattice by factor 2 D=3 . The idea behind this is that after having carried out sufficient cooling in the first place the field modes that resolve individual lattice spacings are already highly suppressed and thus the blocked lattice will contain essentially the same information as before. While in the literature it has been suggested to use more elaborate schemes, where also staples around the coarsened links are considered in this study we refrain from such an additional smearing step. The practical benefit of coarsening is that any subsequent cooling proceeds much faster. We have checked in several cases that we end up in the same vacuum sector with and without cooling.
In summary, at each step in real-time we prepare a copy of the gauge links, cool them for a certain flow-time, carry out one blocking before finally cooling towards the point at which the behavior of N CS (t) − N CS (0) stabilizes to unit steps. Since even with the improved action overcooling can still occur we have to empirically determine the appropriate cooling depth for the lattices at hand. Note that this dissipative approach should not be confused with the recently introduced gauge cooling in the context of complex Langevin simulations [76,77], as it constitutes more than a simple gauge transformation and actually changes the value of the classical action.
Note also that since we here directly evaluate N CS from F µν , which is defined from plaquettes of spatial links on neighboring time-slices, the explicit electric fields of the hamilton dynamics do not enter. In turn, besides cooling the links there is no need to change the electric fields themselves. For comparison purposes we have implement the electric field cooling from Ref. [41] which did however not reproduce our results, a comparison with the original code was unfortunately not possible anymore.

Effective gauge dynamics in the presence of chiral fermions
Now that we have spelled out how the gauge fields evolve classically and how to determine their topology we can implement Eqs. (2.15) and (2.17) on the lattice. While the update of the gauge links remains unchanged the equation of motion for the electric field changes due to the presence of the Ohmic current, the corresponding thermal noise, and the dissipationless anomalous current, . We explicitly implement the stochastic evolution and do not resort to heat-bath updates [74] nor the electric field refresh method [75] previously deployed in the literature.
At this point we need an explicit expression for the magnetic field, for which we take B i = 1 2 ijk F jk (i.e., B 0 = −F 12 , B 1 = −F 20 , B 2 = F 01 ) with F jk obtained from the same clover approximation used in the determination of the topological charge. . Schematic view of the evolution of the soft gauge fields E a k (x, t) and U k (x, t) in the presence of a dynamical chiral imbalance n 5 (t) coupled effectively via an anomalous current. At each step in the leap-frog-like real-time evolution we prepare a copy of the gauge links and cool them, including a single coarsening step, in order to compute the topological contribution to dN CS dt ∝ F µνF µν . The change in N CS enters the update of n 5 (t) which in turn via µ 5 (t) influences the update of the electric fields at this step. The cooled copy of the configuration is kept in memory for one more time-step, since the electric part of F µν requires the evaluation of temporal plaquettes.
The evolution of the chiral imbalance n 5 (t) is driven by changes in the topology and we naively discretize the continuum equations (2.15). Since no spatial dependence of n 5 (t) is considered here its derivative can be related to the change in the topological charge

Numerical simulations and their discussion
With the lattice discretized formalism in place we are ready to embark on the investigation of the effective gauge dynamics in non-Abelian SU(2) Yang-Mills theory 6 in the presence of a chiral imbalance with N f = 2 flavors. Let us start first with a cross-check, whether our numerical implementation is able to reproduce the well-established results for the Yang-Mills sector. The first test concern the evolution of the topological charge in a thermal ensemble of electric fields and gauge links propagating via pure Hamilton dynamics, as was originally investigated in Ref. [74]. The second testing ground is the driven dynamics of color fields in the presence of an Ohmic current and stochastic noise, for which in Ref. [75] the sphaleron rate has been determined.
We choose our lattice geometry according to the insight gained in Ref. [75]. It showed that the sphaleron rate in both SU(2) and SU(3) can be reliably determined on three-dimensional lattices with an extend of fifteen points, hence we choose N = 20. In order to not having to deal with too rough field configurations we furthermore choose the lattice coupling β L = 8 with a s = 1, g = 1, which corresponds to the choice of aN C g 2 T = 1 in Tab. 1 of Ref. [75]. In their study this choice of parameters already produced a result for the sphaleron rate within 1σ of the continuum extrapolated result. We set the time discretization of the Hamilton dynamics fine enough so that the same a t can be used in all runs, even in the presence of fast dynamics, expected in the presence of a chiral imbalance. A posteriori we found that a t = 0.0375a s provides an appropriate amount of temporal resolution and at the same time allows conservation of magnetic and electric energy on the single percent level.

Thermal Yang-Mills dynamics
Starting from thermal initial field configurations obtained by the prescription laid out in Sec. 3.1, we evolve E a k (x) and U k (x) via Eq. (3.7). Our goal is a robust determination of N CS (t) − N CS (0) for which it is necessary to set up an appropriate cooling procedure. The use of the fourth-order Runge-Kutta scheme allows us to use a relatively large cooling step size of dτ = 0.15a s . It furthermore turns out that a single coarsening after N coarse = 100 cooling steps does not adversely affect the determination of the current topological sector. In order to avoid overcooling we perform the computation of the topological charge for different cooling depths starting with the same initial condition, as shown on the left of Fig. 3. The appropriate cooling depth L cool = N cool dτ is chosen such that N CS (t) − N CS (0) has approached the unit step function. As one can clearly see further cooling eventually reduces the unit step behavior and will ultimately remove all structure. With the optimal setting of N cool = 272 for the pure Hamilton dynamics, we proceed to measure the topological charge on N Ham conf = 40 ensemble realizations, a selection of which is shown on the right of Fig. 3. 7 The sphaleron rate Γ is obtained from the appropriately rescaled variance of the topological charge γ(t) defined in Eq. itself a fluctuating quantity with an uncertainty we need to estimate. To this end we carry out resampling of γ(t) by repeatedly drawing a random selection of twenty of the forty ensemble realizations we computed numerically. This procedure leads to the blue error bands in the left panel of Fig. 4.
The sphaleron rate can be estimated in two similar ways. Both assume that the topological charge diffuses, i.e., its variance will grow linearly with time. In the left panel of Fig. 4 we fit the values of γ(t) using a straight line without intercept. Fixing the function at the origin is a quite restrictive choice, which leads to relatively small error estimates when changing the upper end of the fitting interval. To estimate the inherent systematic error we assign error bars on the slope, such that the meandering orange average lies within 1σ at late times t > 2000,  Figure 4. Numerical results of pure Yang-Mills theory: (Left) Estimating the sphaleron rate Γ from the slope of γ(t). We resample twenty times by drawing a subset of twenty ensemble realizations from which the average (orange) and variance (blue) of γ(t) follows. Assuming a diffusive process we fit with a linear function without intercept (green). The gray curve denotes the value for Γ = 64.5 observed in Ref. [74]. (Right) The behavior of γ(t)/t which will asymptote to the constant Γ at late times. Different curves correspond to a different number of ensemble realizations used. Note that more than ten runs are required obtain convergence.
Γ Ham = 44.5 ± 12, (4.1) to be compared to Γ Ham [74] = 64.5 ± 9.7 in Ref. [74]. The error bands of our estimate overlap with those of the published result from Ref. [74]. We note however that our average value lies systematically lower than Γ Ham [74] . Three possible reasons for this difference exist, the first lies in the cooling prescription, the second in the operator used to define E · B and the third in a different number of configurations used in the ensemble averaging. Instead of the calibrated cooling of Ref. [74] with standard Wilson action, we here use the gradient flow from the improved action to directly reach the vacuum in each time step. When measuring the field strength tensor on the other hand we deploy the standard clover approximation, while Ref. [74] is based on an improved operator. Last but not least the number of ensemble realizations used in the estimate influences the value of Γ. While Ref. [74] does not explicitly mention the number of ensembles used, in an older study [41] to which it refers, 20 to 60 realizations have been generated.
In the right panel of Fig. 4 we show the values of γ(t)/t computed from a different number of ensemble members. It seems necessary to include at least twenty of them before convergence to the correct late time result is possible. Fitting a constant to this quantity at t > 2000 constitutes the second way of estimating Γ. We obtain a similar value as in (4.1) consistent within error bars. If Γ is obtained from a smaller number of ensemble members it shows slightly larger values and thus lies closer to the results reported in Ref. [74].

Langevin-type effective theory without chiral imbalance
The next step is to investigate a scenario where hard modes are coupled to the gauge fields via an Ohmic current and stochastic noise as discussed recently, e.g., in Ref. [75]. Starting from the same initial thermal configurations as in Sec. 4.1 we evolve the electric fields via Eq. (3.22), neglecting the anomalous term at this stage. The value of the color conductivity is set to σ c = 1, as its influence can be simply absorbed into a rescaling of time [75]. The interplay between dissipation and the driving noise by construction leaves the energy content of the magnetic sector unchanged, which is fulfilled in our implementation within statistical errors, as can be seen form the two top curves in Fig. 5. There is a slight systematic shift of the energy to lower values in the driven case, which we take as contribution to our systematic uncertainties.
The electric fields on the other hand show a pronounced difference in energy reflected in the two lower curves in Fig. 5. The energy in the electric sector exhibit a sharp rise at early times before settling at a stationary value at t 10. Note that in accordance with the assumptions in the derivation of the effective dynamics (see Tab. 2.1) the electric energy never grows beyond its magnetic counterpart E el < E mag . The initial rise in E el is not particular to our numerical implementation, as the stochastic update deployed in Ref. [75] also leads to a similar increase. It tells us that if we are interested in investigating possible fast dynamics from anomalous effects in the following, we should first let the driven system equilibrate to its stationary state before switching on a chiral imbalance.
The very different fluctuation content compared to pure Yang-Mills theory also requires us to choose a different cooling scheme for the robust determination of N CS (t)− N CS (0). With the same dτ = 0.15 as before and a single coarsening after N coarse = 100 steps, we find from Fig. 6 that a cooling depth of L cool = 242dτ leads to a stable unitstep behavior in the topological charge. The behavior of the cooled N CS (t) − N CS (0) for a selection of ten different initial conditions is given in the right panel of Fig. 6.
We estimate the value of the sphaleron rate in this scenario again from a linear fit to the resampled values of γ(t) as shown in the left panel of Fig. 7. The systematic errors are assigned such that the fluctuating sample average (orange) lies within the 1σ cone of the fit. The value we obtain here is  Figure 5. Numerical results of nonchiral Langevin theory: Magnetic and electric energy density T 00 of the pure Yang-Mills system compared to Yang-Mills with dissipatively coupled fermions. By construction the energy density in the magnetic sector does not deviate significantly between the two cases, while the presence of the stochastic noise leads to a marked increase in T 00 el before settling to a new steady state at around t 10. Note that T 00 el < T 00 mag is always fulfilled.  which agrees with the result Γ Ohm [75] = 22.10 ± 0.62 in Ref. [75] within our relatively large error bars. Note that the amount of cooling required in the driven case is actually less than for pure Hamilton dynamics and also the dependence of the sphaleron rate on the number of used ensemble members is weaker, as can be seen from the right panel of Fig. 7, where γ(t)/t is shown.

Langevin-type effective theory with chiral imbalance
The genuinely new ingredient in this study is the inclusion of the chiral imbalance n 5 (t) as an effective but nevertheless explicit degree of freedom. Its back-reaction onto the gauge fields arises from the anomalous current proportional to the magnetic field entering Eq. (3.22). We start the real-time evolution of the system from the same initial configurations as before but at first t < t 0 = 15 with no chiral imbalance present. As we have seen in Sec. 4.2 this is necessary for the system to settle into the stationary state, where the electric energy of the dissipatively coupled system does not change with time anymore. At t 0 = 15 we then switch on n 5 (t 0 = 15) = n 0 and observe the subsequent behavior of the system. Note that since the anomalous current is dissipationless the character of the stochastic noise has not changed and we thus find that the same cooling as in Sec. 4.2 with L cool = 272dτ suffices for a robust determination of N CS (t) − N CS (0).
The first question to ask is whether the interplay of the gauge sector with n 5 (t) via the anomalous current influences the sphaleron rate of the system even without a as shown in Fig. 8. Since the error bars are still relatively large, we cannot distinguish the two values statistically. It would however be interesting to reduce the uncertainty in a future study to understand whether this apparent agreement is accidental or a reflection of physics. As discussed in the introduction section, analytic considerations have lead to the conclusion that in the presence of a chiral imbalance fluctuating gauge fields can develop instabilities which drive rapid dynamics in the magnetic sector at early times [10]. The back-reaction of the system on the unstable modes will eventually lead to a saturation of their energy and the system is expected to enter a new steady state.
In Fig. 9 we show the electric and magnetic energy of the system at early times for different values of n 0 = 0, 12.5, 25, 50, 100, 500, 1000. They indeed exhibit a rapid increase instantly after switching on the imbalance at t 0 = 15. Looking at the magnetic sector at intermediate values of n 5 , it appears as if three different regimes can be identified. A very rapid early rise, followed by a slightly longer linear phase and eventually saturation to the new steady state. Not only does the value at which the energy saturates increase with larger values of n 0 is also reached a slightly earlier times indicating a more rapid instability at larger imbalances. The electric fields also show a rise over the same timescales, which however is smaller in magnitude. Interestingly the  Figure 9. Numerical results of chiral Langevin theory: (Left) Magnetic energy density of the gauge fields at early times with the chiral imbalance switched on at t 0 = 15. One observes a rapid rise, clearly ordered with increasing values of n 0 . A simple inspection by eye hints at three different regimes, a very brief but rapid one shortly after t 0 followed by a linear increase until saturation is reached; (Right) The corresponding electric energy density, which has stabilized at t 10 from its initial rise before following the magnetic energy density with a rapid increase after t 0 .
electric energy never overtakes its magnetic counterpart, an a posteriori justification of the assumption that E el < E mag made in the derivation of the effective dynamics (see Tab. 2.1).
Let us have a look at the behavior of the topological charge in the presence of a chiral imbalance. As an example we showcase N CS (t) − N CS (0) at n 0 = 25 in the left panel of Fig. 11. Before t 0 the same changes in topology occur as in the case without anomalous current but once n 5 is present a clear drift among N CS (t)−N CS (0) is induced. The drift is topology driven, as it occurs via well separated transitions of unit steps, which when taken together lead to a linear increase of the topological charge over time.
For larger values of n 5 even more transitions per time occur, so that eventually our temporal resolution will not be able to accommodate all of them and we will end up with a saturation of the drift. This is exactly what we find in the left panel of Fig. 11, where the averaged N CS (t) − N CS (0) is shown. Our choice of a t = 0.0375 appears to give reliable results up to n 0 = 100, beyond which saturation effects become significant.
One expects that the induced changes in the gauge field will diminish the chiral imbalance, in line with the requirement of helicity conservation. We have explicitly checked it to be fulfilled in our numerical implementation. And indeed the increase in N CS (t) − N CS (0) thus leads to a decrease of n 5 (t) over time as can be seen on the right Estimating the sphaleron rate Γ from the slope of γ(t). We resample twenty times by drawing a subset of nineteen ensemble realizations from which the average (orange) and variance (blue) of γ(t) follows. Assuming a diffusive process we fit with a linear function without intercept (green).
of Fig. 11. It is our particular choice of σ c = 1, β L and n 0 that yields a rather slow decrease in n 5 (t) here. Interestingly, when rescaled to unity at t 0 the relative decrease is weaker, the larger the initial n 0 is.
While currently too costly, an investigation of the very late time behavior of the system would be highly desirable. In particular it would be interesting to observe whether, the chiral imbalance is fully depleted in the end and to what state the system relaxes if n 5 = 0 is eventually reached. As we saw that the system is stable for n 0 = 0 and no replenishment of the imbalance occurs spontaneously, once n 5 (t) is depleted, we expect to eventually return to a state of vanishing drift, diffusing around a large but constant value of the topological charge. Now we can proceed to inspect the sphaleron rate in the presence of a chiral imbalance, the relevant plots are given in Fig. 12. For completeness we plot the values for γ(t)/t in the left panel to qualitatively illustrate that there is an ordering in the sphaleron rate with n 0 . The actual values of Γ are obtained as before from linear fits to γ(t), summarized in Tab. 2 and plotted in the right hand pane of Fig. 12.
While it is not easy to distinguish from the left panel of Fig. 12, the values for Γ     below this trend. A reason for this deviation might be that at n 0 = 100 the drift in N CS (t) − N CS (0) has become too fast to be adequately resolved by our temporal step size. I.e. it is related to the saturation of N CS (t) − N CS (0) itself, observed in Fig. 11.
As all previous results have been obtained at a single lattice spacing we need to make sure that the phenomena we observe are not simply artifacts of the finite lattice spacing discretization. To probe the small spatial lattice spacing limit within a classical statistical simulation in thermal equilibrium we can increase β L . In order to make sure that we retain the relevant physics of the infra-red regime, the lattice volume has to increase at the same time so that β L /N s = const.
In the normalized magnetic energy shown on the left of Fig. 13 we find that as the chiral imbalance is switched on, the onset of initial fast dynamics occurs independent of β L . The relative value at which the instabilities saturate on the other hand grows slightly with β L but seems to saturate already at β L = 20. In absolute terms this limiting value decreases, as the energy in the system is diminished by lowering the temperature.
Looking at the average of the Chern-Simons number on the right panel tells us that the onset of the linear rise seems to shift slightly to later values. The slope becomes steeper as we go from β L = 8 to β L = 20 but it appears to have saturated at the larger value. Note that we have not adapted the cooling prescription for different β L , so that a quantitative interpretation of the drift in Fig. 13 should not be attempted. What we can however infer is that the phenomenon of chiral instabilities survives in the a s → 0 limit.
The other limit concerns the temporal step size. Unfortunately in the context of this study we do not have the computational resources to perform a thorough analysis of the a t dependence of the late time sphaleron rate. On the other hand, according to the findings of Ref. [75] the difference in Γ between a t = 0.1 and a t = 0 is below 2%. This gives us confidence that with a t = 0.0375 the time discretization error is insignificant compared to the other systematic uncertainties.

Discussion and conclusion
We started out with two sets of questions. The first one concerned the existence and early-time evolution of the chiral instability. The second one aimed at understanding the fate of fermion chiral imbalance: how n 5 and N CS behave at late times. For both cases this exploratory numerical study has provided qualitative insight.
The basis for our simulation of the real-time evolution is a chiral Langevin theory [11] in which the soft gauge sector is coupled to hard chiral fermions via both an Ohmic current, as well as a chiral magnetic current. The interaction between the chiral imbalance and the topology of the gauge sector is dynamically included. We have tested the numerics of cooling, required for the lattice determination of N CS , as well as those for solving Hamilton equations of motion in the pure Yang-Mills sector for N s = 20 and β L = 8. While the values we obtain for the diffusion of topological charge lie quite close to those reported in Ref. [74], we note that our results are systematically smaller. The reason for this deviation might be related to our choice of fully cooling to the vacuum at each time-step, the use of the standard clover approximation for the field strength tensor and a different number of ensemble realizations used. On the other hand our sphaleron rate in the driven system without anomalous effects, agrees within 1σ with the corresponding reference value obtained in Ref. [75].
In the chiral Langevin theory we find that for an initially depleted n 5 the system is stable and on average no chiral imbalance is generated. The system diffuses around the topological sector it started from and shows within errors the same diffusion constant as if the anomalous term is switched off by hand. Increasing the statistics for these two cases in the future is highly desirable, as it will allow us to reduce the still relatively large error bars on the sphaleron rate. In turn we will be able to settle, whether the 1σ difference between the Ohmic and n 0 = 0 rate diminishes or indeed different diffusion occurs. We actually expect that the picture of conventional diffusion of Chern-Simons number needs to be modified even for n 0 = 0. Once an individual realization of the system evolves for long enough times that changes in topological charge and thus fermion chiral imbalance can accumulate, the effects of the anomaly can become relevant and reduce the growth of (N CS (t) − N CS (0)) 2 . The systematically smaller central value of the sphaleron rate in the presence of the anomalous term compared to the purely Ohmic theory might be an indication of this mechanism.
In the presence of a chiral imbalance, we find that modes in the soft gauge fields indeed become unstable. A chiral instability ensues with a rapid increase of magnetic and electric energy, which is fed by the hard sector. For intermediate values of initial n 5 (t 0 ), an inspection by eye hints at the existence of three different regimes. A very fast initial rise of T 00 mag is followed by a slower linear increase that saturates to an apparent constant. Note that n 5 at that time is still appreciably different from zero and thus will continue to feed topology changing processes in the gauge sector. Conservation of energy in the overall system implies that further energy should be accumulated in the gauge fields. Since for our choice of parameters however the decrease of n 5 proceeds very slowly (see Fig. 11) the accompanying minute changes in energy might still be masked by statistical noise. The value of energy at which the instability saturates increases with the size of the initial imbalance n 5 (t 0 ). Not only is more energy deposited in the gauge fields but the time it takes to reach the saturation also diminishes. When decreasing the spatial lattice spacing towards zero by increasing β L we find that the instability prevails and hence is not just an artifact of discretization.
Interestingly the presence of a finite n 5 results in a clear linear drift over time in the average Chern-Simons number. Inspecting the individual ensemble realization affirms its topological nature, as it is composed of well separated unit step transitions (see. Fig. 10). A deviation of N CS (t) − N CS (0) from zero is not surprising if n 5 = 0. The linear time dependence exhibited here on the other hand is puzzling, as the non-Abelian character of the gauge fields does not immediately suggest such a simple relation. Actually if we consider as model of the evolution of N CS and n 5 a simple coupled random walk, the linear dependence on the initial value is just a manifestation of the early time rise of an exponential behavior that will eventually asymptote to its equilibrium value. Comparing the rise in energy (Fig. 9) and the changes in topology (Fig. 10) we see that the onset of the drift is located just after the saturation of the rise in the magnetic and electric field energy.
Due to helicity conservation the changes in N CS act opposite to the imbalance by which they are fed. This reduces n 5 (t) monotonously. Even though we did not have the computational resources to follow the time evolution up to the full depletion of n 5 , the stability of the system in the absence of n 5 tells us that eventually the drift in N CS will abate and we do not expect a subsequent spontaneous replenishment of n 5 to occur: The linear rise observed at early time cannot be maintained forever and the system will begin to diffuse around a different topological sector than the one it started from. Whether the sphaleron rate at these times reverts to the one expected for the n 5 (t 0 ) = 0 case needs to be ascertained in a future study.
At intermediate times (t = 100 − 4500), i.e., after the saturation of the initial instability but before the value of n 5 has changed appreciably we find that the sphaleron rate shows a clear dependence on n 5 (t 0 ). Once the initial imbalance becomes too large and the drift between topological sectors becomes too rapid to be resolved by our temporal spacing, a saturation of the rate is expected and also observed (n 5 (t 0 ) ≥ 100). For smaller values of the initial imbalance we find a clear linear dependence of Γ on n 5 (t 0 ). The interpretation of this finding and its possible phenomenological consequences in the context of heavy-ion collisions and early universe baryogenesis are work in progress.
Rapid developments in the simulations of nonequilibrium systems in the classical statistical regime bodes well for further insight into the questions investigated in our exploratory study. On the one hand state-of-the-art simulations (see, e.g., Ref. [78]) now incorporate by default the expanding geometry of the heavy-ion collision region and make close contact with the saturation picture for the incident nuclei. The fact that the gauge field modes do not fill up the spectrum up to the cutoff in these simulations contrary to the thermal case, makes them systematically better controlled. On the other hand developments are underway to further develop the treatment of explicit fermionic degrees of freedom from first principles [27,28,[79][80][81] and to eventually include quantum corrections to classical statistical simulations. Progress in those directions most surely will provide us with new tools for understanding anomalous effects at early times.