Chiral instabilities&the onset of chiral turbulence in QED plasmas

We present a first principles study of chiral plasma instabilities and the onset of chiral turbulence in QED plasmas far from equilibrium. By performing classical-statistical lattice simulations of the microscopic theory, we show that the generation of strong helical magnetic fields from a helicity imbalance in the fermion sector proceeds via three distinct phases. During the initial linear instability regime the helicity imbalance of the fermion sector causes an exponential growth(damping) of magnetic field modes with right(left) handed polarization, for which we extract the characteristic growth (damping) rates. Secondary growth of unstable modes accelerates the helicity transfer from fermions to gauge fields and ultimately leads to the emergence of a self-similar scaling regime characteristic of decaying turbulence, where magnetic helicity is efficiently transferred to macroscopic length scales. Within this turbulent regime the evolution of magnetic helicity spectrum can be described by an infrared power-spectrum with spectral exponent $\kappa$ and dynamical scaling exponents $\alpha,\beta$, which we determine from our simulations.

One important aspect of anomalous transport concerns the question how chirality is transferred between gauge fields and fermionic degrees of freedom. While in QCD plasmas chirality transfer can be efficiently accomplished by sphaleron transitions between different topological sectors of the non-Abelian gauge theory [20], the situation is markedly different in abelian plasmas. Quantum Electrodynamics (QED) is topologically trivial and a different mechanism has to be invoked to convert fermionic chirality into magnetic helicity (and vice versa). In this context, a novel type of 'chiral' plasma instability has been suggested as a viable mechanism, whereby a chirality imbalance in the fermion sector can generate helical magnetic fields that exist on macroscopic length scales [21,22].
In this article we present a comprehensive study of chiral instabilities using microscopic real-time lattice simulations of strongly-coupled QED plasmas. Starting from a helicity imbalance in the fermion sector, we employ a classical-statistical description [48][49][50][51][52] to simulate the subsequent non-equilibrium evolution of the system from first principles. We demonstrate that chiral instabilities in QED-like theories follow a characteristic pattern of quantum many body systems subject to instabilities [53][54][55][56][57][58], where the exponential growth of unstable modes leads to the emergence of turbulent behavior. Specifically, the evolution follows a typical sequence, where in the initial regime of linear instability, a well-defined set of helical magnetic field modes exhibit exponential growth. Non-linear effects induce the secondary growth of a wide range of left and right handed magnetic fields modes, until the unstable growth saturates when the magnetic field spectrum exhibits non-perturbatively large occcupation numbers, with a significant left (L) versus right (R) helicity imbalance. Subsequently, the plasma enters a regime of decaying turbulence where the spectrum of magnetic field modes shows a self-similar scaling behavior associated with an inverse cascade of magnetic helicity. Based on our microscopic simulations, we are able to characterize the entire sequence of events, starting from the extraction of growth-(and decay-) rates of helical gauge field modes in the primary and secondary instability regimes all the way to the turbulent scaling regime for which we extract the relevant far-from-equilibrium scaling exponents for the first time.
Simulation technique. We perform real-time simulations of N f degenerate flavors of quantum Dirac fermions, coupled to classical-statistical gauge fields [59][60][61][62][63][64]. By numerically solving the coupled set of of Dirac and Maxwell's equations for electro-magnetic fields E i and is the covariant derivative. We include the effects of electro-magnetic fields on the fermions sector as well as the non-linear back-coupling of fermion currents j ν = 1 2 [Ψ x γ i , Ψ x ] on the dynamical gauge fields.
By discretizing the theory on a N 3 s spatial lattice, the fermion field operatorΨ x (t) becomes finite dimensional, and the solution to the operator Eq. (1) can be constructed from the solution of Dirac equation for a set of 4N 3 s independent wave-functions [59,64]. We use a compact Hamiltonian lattice formulation of QED [65] with O(a 3 ) tree-level improved Wilson-fermions, which is crucial for non-pertubative studies related to chiral anomaly, as detailed in [64]. Our classical-statistical simulations are accurate to leading order in e 2 , but all orders in the coupling between gauge and matter fields e 2 N f . We specifically work in this regime of strong interactions between gauge and matter fields e 2 N f = 64, which is significantly larger than in single flavor QED (e 2 N f ≈ 0.09) to resolve all relevant scales by including quantum fluctuations in our lattice simulations. We prepare the system as a charge neutral but chirally imbalanced Fermi gas, by specifying initial fermion occupation numbers according to the Fermi-Dirac distribution n u/v h (p) = 1/(e (ωp−µ h h)/T + 1) with helicity chemical potential µ h , where ω p = ± p 2 + m 2 labels the energy levels for particle and anti-particles respectively. We work in the low-temperature regime (µ h T ) and use T /µ h = 1/8 in our simulations. Similarly, we consider vacuum initial conditions for the gauge field sector, which are respresented by a classical-statistical ensemble of fluctuating fields (see [66] for details). Our simulations are performed on different lattices with N 3 s = 32 3 and 48 3 sites and µ h a s = 2/3, 1, 1.25, 1.5 and close to the chiral limit ma s = 5 · 10 −4 , where a s is the lattice spacing. We express our results in terms of dimensionless quantities in units of the initial helicity chemical potential µ h . Chiral instabilities Starting from an initial helicity imbalance (µ h > 0) in the fermion sector, the chiral plasma  instability triggers an exponential growth of gauge fields with right-handed (circular) polarization, as can be seen from Fig. 1, where we present the evolution of the left (L) and right (R) handed components of the magnetic fields n L/R B (t, p) ≡ ±|B L/R (t, p)| 2 /|p| for selected momentum modes |p|/µ h ≈ 0.1 − 1. Separating the field into left and right handed components according to the helicity projection B L/R (t, x) = |p|±ip× 2|p| B(t, p) (see supplementary material for details of the lattice definition), one finds that at early times only the right-handed components experience exponential growth within a narrow momentum range; all other modes show a damped oscillatory behavior as can be qualitatively expected from the interplay of electric-magnetic fields and currents in a conducting medium. We further quantify this behavior in Fig. 2, where we show the growth rates (blue symbols) for the primary unstable and damping rates (orange symbols) of initially stable modes, extracted from an exponential fit to the evolution of the occupation number n L/R B . Shown results for the growth and damping rates are quantitatively consistent across different lattices sizes N 3 s = 32 3 , 48 3 and lattice spacings µ h a s = 2/3, 1, 1.25, 1.5, with a maximal growth rate γ 0 ≈ 0.07µ h for right handed modes with |p| ≈ 0.5µ h .
While the exponential growth (damping) of right-(left-)handed modes sets in almost directly, after a short delay of µ h t ≈ 20 − 50 due to the initial quench [67], the evolution continues in this fashion until non-linear interactions between unstable modes induce secondary instabilities [57,68]. During this second phase, which in Fig. 1 occurs around µ h t ≈ 250, a large range of left and right handed momentum modes starts to exhibit exponential growth with strongly enhanced growth rates γ secondary ∼ 2 − 3γ 0 , until around µ h t * ≈ 350 the instability saturates and the exponential growth terminates.
Energy & Helicity transfer Conserved quantities are shared and transferred between fermions and gauge field throughout the evolution of the system. Clearly, the total energy is conserved and receives contributions from the and fermion sector One also finds an approximate conservation law for the net chirality of the system. By means of the axial anomaly [69,70] and magnetic helic- is conserved in the chiral limit (m → 0), and we have checked explicitly that for the small values of m considered, dissipative effects due to finite fermion mass in the second line of Eq. (3) are negligible over the time scale of our simulations. However, it is important to point out that lattice discretization effects can spoil the net chirality conservation; this is precisely the reason why we use tree-level improved operators such that the axialanomaly relation in Eq. (3) is in the continuum-scaling regime (see e.g. [64]).
Simulation results for the individual contributions to the energy density (left panel) and net chirality (right panel) are compactly summarized in in Fig. 3. Different points in each panel show the results for different lattice sizes and spacings, and we have shifted the horizontal axis of the individual data sets to account for the residual discretization dependence in the time µ h t * where exponential growth saturates [71]. While initially the dominant contribution to energy density and net chirality resides in the fermion sector, the chiral plasma instability leads to an exponential growth of electric and magnetic components of the energy density. Growth rates of volume averaged quantities E 2 (t), B 2 (t) and n h (t) are dominated by the growth rate γ 0 of the maximally unstable mode as indicated by the dashed line ∝ e γ0t . Despite the exponential increase, only a small fraction of the total energy density e 2 ε tot ≈ 0.033 e 2 N f µ 4 h is transferred from fermions to electro-magnetic fields. It is also interesting to observe that throughout the evolution, the magnetic field strength exceeds the electric one by at least one order of magnitude, B 2 E 2 , indicating the presence of strong interactions between gauge and matter fields.
When considering the balance of the net chirality in the plasma, a manifestly different picture emerges. Despite the fact that the (continuum) anomaly relation in Eq. (3) is violated at finite lattice spacing, our use of operator improvements significantly reduces discretization effects allowing us to perform controlled continuum extrapolations, which are consistent with conservation of the net chirality, as is indicated by the shaded bands in the right panel of Fig. 3. While the continuum extrapolation is subject to relatively large uncertainties due to the available size lattices, we can safely infer that a substantial amount of the axial charge density of fermions n 5 (t = 0) ≈ 0.039µ 3 h is transferred to magnetic helicity density over the course of the evolution. Specifically, we find that for the two largest lattices available, the mag- Chiral Turbulence. Subsequent to the saturation of unstable growth, the plasma enters a turbulent regime characterized by a much slower time evolution, which we analyze in terms of the magnetic field and fermion spectra depicted in Fig. 4. Spectral distributions of the net-helicity in the gauge field sector, i.e. the difference between occupation numbers of left-and right-handed modes ∆n B (t, k) = n R B (t, k) − n L B (t, k), at different times µ h t = 100, 250, 350, 400, 500, 600 are presented in the top panel. The lower panel shows the spectra of left-and right-handed fermions extracted from equal time correlation functions in Coulomb gauge by performing the appropriate projections onto positive/negative helicity states (see supplementary material for details). Starting from the linear instability regime for µ h t 250, the nethelicity in the gauge sector shows an exponential growth within a limited range of wave numbers |p| 0.8µ h , while left and right handed fermion spectra remain essentially unchanged with distinct sharp Fermi surfaces separated by the helicity chemical potential. Secondary growth of instabilities between µ h t ≈ 250 − 350 leads to a strong population of magnetic field modes at low and high wave-numbers. Over the same period of time the rapid changes in the gauge field sector are accompanied by a significant heating and depletion of the helicity imbalance in the fermion sector, as can be inferred from the softening of the Fermi surface along with narrowing of the gap between left and right handed modes. Eventually, for µ h t > µ h t * ≈ 350, the growth of the chiral instability saturates, and the evolution slows done considerably compared to the rapid changes at earlier times.
In the turbulent regime the spectrum of magnetic helicity ∆n B (t, k) exhibits a self-similar scaling behavior, which is illustrated in the top right panel of Fig. 4. Upon re-scaling, the spectra at different evolution times µ h t = 450 − 600 are all found to collapse onto a single scaling curve f s (|p|). While a detailed characterization of the scaling function f s (|p|) is beyond the scope of this work, we note that for intermediate momenta, the scaling function f s (|p|) features a ∼ p −κ power law behavior with a large scaling exponent κ = 10.2 ± 0.5 illustrated by the gray dashed line. Due to self-similarity, the late time evolution of the spectrum of magnetic helicity can be characterized in a compact form with scaling exponents α, β and scaling function f s , where τ ≡ µ h (t−t * ) is a dimensionless time variable with respect to the reference time µ h t * ≈ 350 for the transition to the turbulent regime. Notably, a self-similiar behavior as in Eq. (4) is characteristic for the late stage evolution of unstable systems, and has been reported previously in a variety of different contexts [54,72,73]. In all of these examples, the initial instability leads to a rapid memory loss of the initial conditions, such that the subsequent turbulent evolution is universal and entirely characterized by α, β and f s , which describe the transport of a conserved quantity across a large separation of scales. Based on a statistical scaling analysis, following the procedures outlined in [73] (see supplementary material for details), we obtain the following estimates for the scaling exponents α = 1.14 ± 0.50 , β = 0.37 ± 0.13 .
Since the magnetic helicity N h can be equivalently expressed as an integral over the net helicity spectrum, N h (t) = d 3 k (2π) 3 ∆n B (t, k), one finds that the approximate validity of the scaling relation α ≈ 3β implies the conservation of the magnetic helicity of the plasma at late times, which is consistent with the behavior seen in Fig. 3. One therefore concludes that the self-similar behavior in Eq. (4) should be associated with an inverse cascade of magnetic helicity, which is transported Visualization of magnetic field lines at times µ h t = 350 (left) and µ h t = 600 (right). Coarsening of the magnetic field due to the inverse cascade of magnetic helicity is observed.
from microscopic ( ∼ µ −1 h ) to macroscopic length scales ( ∼ µ −1 h τ β ). Strikingly, the inverse cascade of magnetic helicity also manifests itself directly in the spatial structure of the field configurations, as illustrated in Fig. 5 where we show stream tracing plots of the magnetic fields. Starting from a significant number of small scale swirls at µ h t = 350 where the unstable growth saturates, one observes a clear coarsening of the magnetic fields towards later times µ h t = 600, where a few swirls fill the entire simulation volume. It is also evident from Fig. 5 that the QED plasma develops sizeable inhomogeneities over the course of the evolution of the chiral instability, which also manifest themselves in other observables such as e.g. vector/axial charge densities which are not shown here but will be discussed in a forthcoming publication [74].
Conclusions & Outlook. We presented an ab-initio study of chiral instabilities and chiral turbulence, based on microscopic real-time lattice simulations of strongly coupled QED. Chirality transfer through chiral instabilities and the subsequent generation of macroscopic helical magnetic fields proceeds in a three stage sequence. Initial primary growth, is followed by secondary growth until the instability saturates, when the gauge field occupation numbers become non-perturbatively large. During the unstable phase, the fermion chirality is significantly depleted and transferred into magnetic helicity, while most of the energy is still carried by fermions. Subsequently, the system enters a turbulent regime, where magnetic helicity is transported to large distances by an inverse cascade. In this regime, the helicity spectrum shows a steep power-law behavior, with exponent κ = 10.2 ± 0.5 and the dynamics follows a self-similar scaling behavior with scaling exponents α = 1.14 ± 0.50 and β = 0.37 ± 0.13. While our study represents the first observation of a self-similar (inverse) cascade of magnetic helicity in non-equilibrium plasmas from the underlying microscopic theory, previous works have reported similar findings based on macroscopic near-equilibrium calculations in anomalous magneto-hydrodynamics [32,75,76]. However, the observed scaling exponents in the non-equilibrium plasma appear to be manifestly different from the corresponding results obtained in nearequilibrium calculations based on anomalous magnetohydrodynamics, which will require further investigation in the future.
While our current study established the dynamics of chiral turbulence close to the chiral limit for strongly coupled QED plasmas, one important next step would be to explicitly verify the universality of our results by varying the coupling strength e 2 N f and further explore the impact of dissipative effects due to finite fermion mass on the chiral turbulent regime. With regards to the dynamics of the Chiral Magnetic Effect in Quantum Chromodynamics (QCD), it would also be interesting to investigate and compare the analogous dynamics in non-Abelian gauge theories, where one ultimately expects the chirality imbalance in the fermion sector to be absorbed into a non-trivial topology of the non-Abelian gauge fields [20]. Nevertheless, on transient time scales there could be an interesting competition between effectively abelian phenomena such as the chiral plasma instability and genuinely non-abelian effects related to the topology of the gauge fields which should be explored further.
Acknowledgements. We gratefully acknowledge the opportunity to participate in the 2018 Intel Knights Landing Hackathon hosted by Brookhaven National Laboratory (BNL) and would like to extend our special thanks to P. Steinbrecher for his help in optimizing the simulation software. The BNL Hackathon was partially supported by the HEP Center for Computational Excellence (KA24001022) and the SOLLVE Exascale Computing Project (17-SC-20-SC). BNL is supported by the Office of Science of the U.S. Department of Energy under Contract No. de-sc0012704. We would like to thank the Institute for Nuclear Theory at the University of We use a compact Hamiltonian formulation of U (1) gauge theory in temporal axial gauge, where we define lattice gauge links U i and lattice electric fields E i as with E i = E a i t a where t a = 1/ √ 2 of the U (1) gauge group. Based on the lattice Hamiltonian, the update rules for the lattice gauge links and lattice electric fields take the form where J a i (t, x) denotes the back-reaction current in the presence of dynamical fermions (c.f. Eq. 13). Gauge links and electric fields are initialized by a gaussian random superposition of plane waves to describe their vacuum fluctuations as detailed in [73].
We use a mode decomposition of the fermion field operator, whereb s andd † s are the time independent fermionic creation and annihilation operators for particles and antiparticles, and Φ s (t, x) are c-numbered wave functions. We initialize the s = 1, · · · , 4N 3 fermion wave functions Φ s (t, x) as plane wave helicity eigenstates, such that s collectively labels the particle/antiparticle nature, helicity and momentum of each mode. Φ s (t = 0, x) = e +ipsx u hs (p s ) particles e −ipsx v hs (p s ) anti-particles (10) where u h (p) and v h (p) denote helicity spinors for Wilson fermions (see Appendix A of [64] for explicit expression). Evolution equations for the fermion wave functions Φ s (t, x) take the form where C 1 = 4/3 and C 2 = −1/6 for leading order treelevel improvement, r W = 1 is the Wilson parameter and we denote Back-reaction currents J a i (t, x) are given by where J a i(n) (t, x) denotes the individual contributions from n-link terms By inserting the operator decomposition in Eq. (9), the expectation values of bi-linear operators are evaluated according to where n s (t = 0) is the initial mode occupancy, given by a Fermi-Dirac distribution n s (t = 0) = 1 exp β(|ω p | − sµ h ) + 1 (16) with helicity chemical potential µ h . Here, the lattice dispersion is ω p = ± p 2 +m 2 for particles and antiparticles respectively, where the lattice momenta are ) and effective mass is given bym = m + 2r W i n Cn as sin (πnq i /N ). Eqns. (6,5,11,13) form a closed set of coupled evolution equations, which we solve numerically using a leap-frog scheme. We note that the coupled set of evolution equations for gauge fields and fermions automatically lead to a conservation of the Gauss law constraint in the sense that G a (t+a t , x) = G a (t, x) such that Gauss law constraint G a (t, x) = 0 is satisfied automatically as long as it is satisfied by the initial conditions. Due to the presence of dynamical fermions, the initial conditions described above do not satisfy the Gauss law constraint automatically. Hence, in order to restore Gauss law at initial time, we follow previous works and subsequently project the electric fields on the constraint surface using the algorithm of [79], and we have checked explicitly that violations of Gauss law stay at the level of machine precision over the course of the entire simulation.
Based on the solution of the gauge links, electric fields and fermion wave-functions, bulk observables, such as individual components of the energy density or vector/axial currents, can directly be calculated from a gauge invariant operator definition (see e.g. [64]). Since in abelian gauge theories the field strength tensor is gauge invariant, we can directly calculate spectra by decomposing the magnetic fields into left/right handed components. Defining the magnetic field B(t, k) as the Fourier transform of the magnetic field in position space B a i (t, x) = ijk ReTr[it a U jk ](t, x), the projection onto left/right handed modes can be performed directly in momentum space according to where the left/right handed helicity projector is given bŷ P ij L/R (k) = |k|δ ij ±ik l lij 2|k| . We employ the following lattice discretization where the eigenvalues of the lattice curl operator are i=x,y,z we denote the Fourier representation of the downward lattice averaging operator and D i = [1 − exp (−ik i )]/a i is the backward derivative, such that in coordinate space the above definitions amount to evaluating the operators B(t, x) and ∇ × B(t, x) at the same lattice point by averaging over neighboring plaquettes. Based on this decomposition, the occupancy of left/right-handed magnetic field modes is simply given by Vice versa, when extracting fermion spectra, we exploit the residual gauge freedom to fix Coulomb gauge at the time of the measurement (see. e.g. [73]) and calculate the helicity projected spectra according to where φ h,p (x) are the free wave functions for particles and antiparticles with helicity h = ±1 and momentum p in Eq.(10).

SCALING EXPONENT EXTRACTION
We perform a statistical analysis of self-similarity following the procedure outlined Appendix D of [73] to determine the most likely values of t * , α, and β. Based on the self-similar behavior of the spectrum of magnetic helicity in Eq.(4), we first define a reference function where τ ref = µ h (t ref − t * ). In order to determine the quality of scaling for a each combination of t * , α, and β, we compare the reference function to the scaled spectrum at a number of different 'test' times, defined analogously as f test (t,p = τ β test p) = log τ −α test ∆n B (τ β test p) , where µ h t test = 450, 475, 500, 525, 550, 575, 600 typically. We quantify the deviations from ideal scaling in terms of withp integration restricted such that the physical momenta p are in the infrared regime 0.1 < p/µ h < 0.5. We then quantify the likelihood that a given set of parameters gives the best overall rescaling by the function where χ 2 min = 1.7 · 10 −3 denotes the minimal value of χ 2 (α, β, t * ) obtained for α ≈ 0.9, β ≈ 0.3, µ h t * ≈ 375 employed in Fig. 4. We visualize the results of the statistical scaling analysis in Fig. 6, where the different panels show the reduced likelihood distributions W (α, t * ), W (β, t * ) and W (α, β) which are obtained by integrating the other variables, e.g. W (α, t * ) = dβ W (α, β, t * ) and similarly for the other components. Based on the selfsimilar scaling behavior, the evolution of the net helicity of the system can be evaluated according to N h (t) p ∆n B (p, t) = τ α p f s (τ β p) = τ α−dβ p f s (p) const × τ α−dβ , (27) such that for α − dβ = 0 the net magnetic helicity is conserved. We find that this scaling relation is approximately satisfied, as indicated by the solid line in the bottom panel of Fig. 6.