Numerical study of chiral plasma instability within the classical statistical field theory approach
NNumerical study of chiral plasma instability within the classical statistical field theoryapproach
P. V. Buividovich ∗ and M. V. Ulybyshev
1, 2, † Institute of Theoretical Physics, University of Regensburg,D-93053 Germany, Regensburg, Universit¨atsstrasse 31 Institute for Theoretical Problems of Microphysics,Moscow State University, Moscow, 119899 Russia (Dated: May 23rd, 2016)We report on a numerical study of real-time dynamics of electromagnetically interacting chirallyimbalanced lattice Dirac fermions within the classical statistical field theory approach. Namely, weperform exact simulations of the real-time quantum evolution of fermionic fields coupled to classicalelectromagnetic fields, which are in turn coupled to the vacuum expectation value of the fermionicelectric current. We use Wilson-Dirac Hamiltonian for fermions, and non-compact action for thegauge field. In general, we observe that the backreaction of fermions on the electromagnetic fieldprevents the system from acquiring chirality imbalance. In the case of chirality pumping in parallelelectric and magnetic fields, electric field is screened by the produced on-shell fermions and theaccumulation of chirality is hence stopped. In the case of evolution with initially present chiralityimbalance, axial charge tends to transform to helicity of electromagnetic field. By performingsimulations on large lattices we show that in most cases this decay process is accompanied bythe inverse cascade phenomenon which transfers energy from short-wavelength to long-wavelengthelectromagnetic fields. In some simulations, however, we observe a very clear signature of inversecascade for the helical magnetic fields which is not accompanied by the axial charge decay. Thissuggests that the relation between inverse cascade and axial charge decay is not as straightforwardas predicted by the simplest form of anomalous Maxwell equations.
PACS numbers: 12.38.Aw, 11.15.Tk
I. INTRODUCTION AND BRIEF SUMMARY
Over the past few decades, real-time instability of thesystem of chiral fermions coupled to dynamical gaugefields has been attracting a lot of attention in variousfields of physics, ranging from astrophysics to condensedmatter physics. This instability manifests itself in the de-cay of the initial imbalance between the densities of theleft- and right-handed fermions at the expense of the gen-eration of magnetic fields with nonzero magnetic helicity(or, in other words, winding number of magnetic fluxlines). In the astrophysical context, the phenomenon ofchiral plasma instability is actively discussed as the mech-anism responsible for the generation and enhancement ofprimordial magnetic fields [1–5] as well as for the transferof magnetic field energy from short to cosmological scales[3, 6, 7].In the context of condensed matter physics, chiralplasma instability was initially discussed and experimen-tally detected as the so-called helical instability of liquid He [8]. It was also considered as a mechanism of spon-taneous magnetization of topological magnetic insulators[9]. In experiments in which chirally imbalanced Weylsemimetal states are created from Dirac semimetals byapplying parallel electric and magnetic fields [10–12] chi- ∗ Electronic address: [email protected] † Electronic address: [email protected] ral plasma instability might manifest itself in the spon-taneous emission of circularly polarized terahertz-rangeelectromagnetic radiation [13].In heavy-ion collisions, chiral plasma instability mightlead to enhanced emission of circularly polarized soft pho-tons [13]. It should be also important for the correct esti-mate of the lifetimes of chirality imbalance and magneticfields [14, 15]. However, the estimates of [16] suggest thatin heavy-ion collisions the volume and the lifetime of thequark-gluon plasma might be too small for the instabilityto develop.The origin of this instability of chirally imbalancedDirac fermions is the Chiral Magnetic Effect (CME)[17, 18] - electric current flowing parallel to the magneticfield in the presence of chirality imbalance. Within thelinear response approximation the contribution of CMEto the electric current is (cid:126)j
CME = σ CME (cid:126)B. (1)The commonly quoted value for the chiral magnetic con-ductivity σ CME is σ CME = µ A π , where µ A is the so-calledchiral chemical potential which parameterizes the differ-ence between the Fermi levels of right- and left-handedfermions and hence also the total axial charge of the sys-tem. The value of σ CME , however, strongly depends onfrequency w and wave vector (cid:126)k of electromagnetic field,and, in the limit of constant and homogeneous magneticfield, on the way in which the limits w → k → a r X i v : . [ h e p - t h ] J u l stability, one can insert it into the classical Maxwell equa-tions, along with the conventional Ohmic current (cid:126)j = σ (cid:126)E ,where σ is the electric conductivity. Assuming the unbro-ken translational invariance both in time and space, wecan write these so-called anomalous Maxwell equations[1, 3, 4, 13, 14, 16, 25–27] in frequency-momentum spaceas iw (cid:126)B = − i(cid:126)k × (cid:126)E, iw (cid:126)E = i(cid:126)k × (cid:126)B − σ (cid:126)E − σ CME (cid:126)B. (2)From these equations we find the following four-branchdispersion relation for transversely polarized plane waveswith the wave vector (cid:126)k = k (cid:126)e [16]: w s,r = iσ s (cid:114) k + rσ CME k − σ , (3)where s = ± r = ± (cid:15) r = 2 − / (1 , − ir,
0) for the electric field (cid:126)E cor-respond to circularly polarized waves with opposite he-licities (handedness) for opposite r .While for nonzero electric conductivity σ the imagi-nary part of w in (3) is always positive and hence corre-sponds to decaying plane waves, nonzero chiral magneticconductivity can also lead to exponentially growing solu-tions if the absolute value of the wave vector k is smallerthan σ CME . From (3) it is also easy to see that for agiven wave vector (cid:126)k , only one of two helical mode willexhibit exponential growth. For example, for µ A > Q A > σ CME >
0) and σ CME > k > E = f e κt cos ( kx ) , E = − f e κt sin ( kx ) ,B = − f kκ e κt cos ( kx ) , B = f kκ e κt sin ( kx ) ,E = 0 , B = 0 , (4)where f is an arbitrary amplitude and κ ≡ − iw = − σ + (cid:113) σ − k + σ CME k . It is important to stress thatsince this solution grows monotonously in time, here weuse the terms “circular polarization”, “handedness” and“helicity” to describe the rotation of the vectors (cid:126)E and (cid:126)B along the x axis, rather than in time. The growthof long-wavelength electromagnetic waves and the de-cay of short-wavelength waves predicted by the anoma-lous Maxwell equations (2) is a novel mechanism forthe inverse cascade in relativistic magnetohydrodynam-ics [3, 6], which transfers energy from long- to short-wavelength helical magnetic fields.The fact that the exponentially growing solution (4)has the helical structure of the form (4) also suggests themechanism which can stop the growth of electromagneticfield at later times. Namely, let us recall that for masslesschiral fermions the time evolution of the axial charge isgoverned by the anomaly equation: ∂ t Q A = g π (cid:90) d x (cid:126)E · (cid:126)B, (5) where the axial charge Q A = Q R − Q L is defined asthe difference between the charges Q R and Q L of theright- and left-handed fermions, g is the electromagneticcoupling constant and we have integrated over space toget rid of the spatial divergence of the axial current. Forsimplicity, in this paper we consider only a single flavorof Dirac fermions with electromagnetic coupling g = 1.For the exponentially growing solution (4) the product (cid:126)E · (cid:126)B is negative: (cid:126)E · (cid:126)B = − f kκ e κt [57]. The anomalyequation (5) then dictates that the time derivative ∂ t Q A of the axial charge is negative. Since we have assumed Q A > µ A >
0, we see that the growing helical solu-tion (4) will result in the decrease of Q A and hence of µ A . This depletion of chirality imbalance should even-tually suppress the chiral magnetic conductivity in (1)and hence slow down or stop completely the exponentialgrowth in (4).However, the above analysis of the chiral plasma insta-bility, which follows [3, 4, 13, 14, 16, 25–28], essentiallyrelies on an assumption that the electric current takesthe form (cid:126)j = σ (cid:126)E + σ CME (cid:126)B with constant ohmic andchiral magnetic conductivities. In reality, both σ and σ CME depend on the frequency and wave vector of elec-tromagnetic field in a nontrivial way [18–24]. One canalso expect a strong dependence of σ and σ CME on thespatial and temporal modulation of the axial charge den-sity, which will in general appear at late evolution times[27]. Moreover, as the instability might lead to quitelarge strengths of electric and magnetic fields, nonlineareffects beyond the linear response result (1) might be-come important. Using linear response approximation todescribe the interactions between the fermions and theelectromagnetic fields is in fact similar to the Lyapunovanalysis of the full quantum evolution, which is in generalnonlinear. What concerns inter-fermion interactions, sofar they were taken into account only indirectly, by usingthe relaxation time approximation [15] or the decoher-ence of the fermionic wave functions [29]. A consistentinclusion of all these effects in the anomalous Maxwellequations (2) would be certainly difficult with approachesbased e.g. on the chiral kinetic theory [25–27, 30], chi-ral hydrodynamics [6, 31] or the Langevin-type effectivetheory [32].This situation clearly calls for a more first-principle de-scription of the real-time dynamics of chirally imbalancedplasma which would overcome these limitations and ap-proximations. In this paper, we report on the numericalstudy of the real-time chiral plasma instability within theframework of the so-called classical statistical field theory(CSFT) [33–38], which captures the first nontrivial orderof the expansion of the full quantum evolution operatorin powers of the Planck constant. CSFT is currently thestate-of-the art method for numerical simulations of real-time quantum evolution. The CSFT approximation isjustifiable as long as the characteristic occupation num-bers of the physically relevant gauge field modes are large.That is, the dynamics of the gauge fields should be almostclassical. On the other hand, the real-time dynamics of
FIG. 1: Two ways of introducing initial chiral imbalance forthe many-body Dirac Hamiltonian. On the left: by intro-ducing the chiral chemical potential µ A in the single-particleDirac Hamiltonian. On the right: by filling more right-handedeigenstates and less left-handed eigenstates (or vice versa). fermions is exact in CSFT. Taking into account that inall previous studies gauge fields were also treated classi-cally, the applicability of the CSFT approach is obviouslywider than that of the previously used approaches. Anapproach very similar to CSFT has been recently usedin [39] to study the real-time dynamics of CME. How-ever, in this work the back-reaction of fermions on theelectromagnetic field, which is the origin of the chiralplasma instability, was not taken into account. The real-time dynamics of axial charge was also studied in 1 + 1-dimensional Abelian Higgs model in the pioneering work[33].Our studies are based on the non-compact formulationof lattice quantum electrodynamics, which avoids poten-tial problems with monopole condensation in the strong-coupling phase [40]. For fermions, we use the masslessWilson-Dirac Hamiltonian which has a low-energy chiralsymmetry. At sufficiently high energies, this symmetry isbroken due to the Wilson term. In the condensed mattercontext, this breaking is a natural feature of any modeldescription of Dirac and Weyl semimetals [41–43]. InSection III we demonstrate that the effect of this explicitbreaking is, however, not very large (see also [44]). There-fore we hope that our results should be also at least qual-itatively relevant in the context of high-energy physics,where the chiral symmetry is exact at the level of theLagrangian, or tends to be exact at sufficiently high en-ergies.In order to introduce the initial chirality imbalance, wehave started the simulations with a state in which moreright-handed eigenstates and less left-handed eigenstatesare filled, as depicted on the right panel of Fig. 1. Such astate is an excited state of the many-body Dirac Hamil-tonian, even in the absence of electromagnetic fields.It is an idealized description of the result of “chiralitypumping” process in parallel electric and magnetic fields[10, 45] or in intense circularly polarized laser beams.Alternatively, we have also considered the introductionof the chiral chemical potential into the single-particleDirac Hamiltonian, which changes the energies of theright- and the left-handed Dirac points (see left panel of Fig. 1). Such initial state has nonzero axial charge butis still the ground state of the many-body Hamiltonianin the absence of interactions with electromagnetic fields.In simulations which started from this state we have notfound any signatures of instability or the transfer of he-licity from fermions to electromagnetic fields. The ax-ial charge density exhibited only small fluctuations ontop of the large mean value. Presumably, the reason forsuch a behavior is that nonzero chiral chemical potentialcorresponds to the physical situation in which our sys-tem is connected to an infinite reservoir of axial charge,which is capable of maintaining its initial value at a con-stant level. Since the anomaly equation (5) holds also atnonzero chiral chemical potential, this also implies thatthe magnetic helicity can only exhibit small fluctuations,possibly related to the violation of the anomaly equationdue to lattice artifacts. Since such behavior is not reallyinteresting, we do not discuss this setup in what follows.The structure of this paper is the following: in Sec-tion II we start with a brief summary of the details ofour numerical CSFT algorithm. In Appendix A we pro-vide a more detailed derivation of this algorithm with abias towards non-relativistic field theories and condensedmatter systems, which might hopefully complement theexisting literature on CSFT (see e.g. [33, 38] for deriva-tions which are more in the spirit of relativistic quan-tum field theory). In this Appendix we also demonstrateexplicitly the absence of any nontrivial Jacobian in theintegration measure in the CSFT algorithm, and discusssome practical aspects of our CSFT simulations on paral-lel computers. In Section III we present the results of thesimulations of the chirality pumping process. First, weconsider chirality pumping in external parallel electricand magnetic fields in the absence of backreaction andverify the validity of the anomaly equation (5) in our nu-merical setup. After that we consider the effect of back-reaction of fermionic current on the chirality pumpingprocess and demonstrate that the dynamical screeningof the external electric field prevents the system from ac-quiring large axial charge density at late evolution times.In Section IV, we consider the decay of the initial chi-ral imbalance and the generation of electromagnetic fieldswith nonzero helicity. In order to trigger the decay, westart simulations with several initially excited modes ofelectromagnetic field. Following the energies of helicalmodes in momentum space, we demonstrate that onlylong-wavelength modes of definite helicity grow and allother modes decay with time. This is a direct numericalevidence of the inverse cascade phenomenon [3, 6, 13, 46]due to chiral plasma instability. We find, however, thatthe dependence of the strength of the inverse cascadeon the initial conditions and parameters of the simu-lations is significantly more complex than predicted bythe anomalous Maxwell equations (2). In particular, inthe simulations which exhibit most rapid growth of he-lical magnetic fields the axial charge does not decay atall. Correspondingly, the mechanism which stops the in-verse cascade in our simulations is not related to the axialcharge decay, again in contrast to the expectations basedon the equations (2) [3, 4, 6, 13–15]. Rather, we observethat in most simulations which do exhibit axial chargedecay the inverse cascade emerges for the electric ratherthan for magnetic fields. In Section V we conclude witha general discussion of our results and an outlook. II. CLASSICAL STATISTICAL FIELD THEORYAPPROXIMATION TO REAL-TIMEEVOLUTION
We consider the many-body fermionic Hamiltoniancoupled to dynamical non-compact electromagnetic fieldson the lattice, so that the full Hamiltonian ˆ H of our sys-tem is ˆ H = ˆ H F + ˆ H EM . The fermionic Hamiltonian ˆ H F reads ˆ H F = (cid:88) x,y ˆ ψ † x h x,y ˆ ψ y , (6)where the labels x , y denote the sites of the three-dimensional cubic lattice, ˆ ψ † x , ˆ ψ x are the spinor-valuedfermionic creation and annihilation operators which sat-isfy the anti-commutation relation (cid:110) ˆ ψ † x , ˆ ψ y (cid:111) = δ xy and h x,y is the massless single-particle Wilson-Dirac Hamil-tonian with the Wilson coefficient r = 1: h x,y = 3 v F βδ x,y + i v F (cid:88) i =1 ( iβ + α i ) e igA x,i δ y,x + e i ++ i v F (cid:88) i =1 ( iβ − α i ) e − igA x − ei,i δ y,x − e i . (7)Here A x,i is the vector potential of the lattice gauge field, e i denotes the unit lattice vector in the direction i , v F is the Fermi velocity, β and α i are the Dirac β - and α -matrices and γ is the generator of chiral rotations: β = (cid:18) (cid:19) , α i = (cid:18) σ i − σ i (cid:19) , γ = (cid:18) − (cid:19) , (8)where σ i are the Pauli matrices. In (7) we have assumedthat the lattice spacing is unity. Thus in what followsall dimensionful quantities are expressed in units of thelattice spacing.The lattice Hamiltonian ˆ H EM of electromagnetic fieldis the straightforward lattice discretization of the corre-sponding continuum Hamiltonian:ˆ H EM = (cid:88) x (cid:88) i =1 ˆ E x,i (cid:88) j = i ˆ F x,ij A x,i J x,i ( t ) , (9)where J x,i ( t ) is the external current (which is required,e.g., to switch the external electric and magnetic fields onand off) and the operator of the magnetic field strengthtensor ˆ F x,ij is defined in terms of the finite differences ofthe vector potential operator ˆ A x,i asˆ F x,ij = ˆ A x,i + ˆ A x + e i ,j − ˆ A x + e j ,i − ˆ A x,j . (10) The operators of the electric field ˆ E x,i and the vectorpotential ˆ A x,i are canonically conjugate and satisfy thecommutation relations (cid:104) ˆ E x,i , ˆ A y,j (cid:105) = − iδ xy δ ij . We im-pose periodic boundary conditions in all spatial direc-tions both for the gauge and the fermionic fields.In our CSFT algorithm, described in detail in Ap-pendix A, we numerically solve the classical equations ofmotion of the electromagnetic field with the Hamiltonian(9) ∂ t A x,i ( t ) = −J x,i ( t ) − (cid:104) j x,i ( t ) (cid:105) −− (cid:88) j (cid:0) F x,ij ( t ) − F x − e j ,ij ( t ) (cid:1) , (11)where the initial value of the time derivative ∂ t A x,i ( t ) | t =0 is the initial value of the electric field E x,i (0) and (cid:104) j x,i ( t ) (cid:105) is the vacuum expectation value ofthe fermionic electric current, which can be calculatedas (cid:104) j x,i ( t ) (cid:105) = Tr (cid:0) ρ u (0 , t ) j x,i u † (0 , t ) (cid:1) , (12)where j x,i = ∂h∂A x,i is the single-particle operator of elec-tric current, u (0 , t ) is the quantum evolution operatordefined by the single-particle Schr¨odinger equation ∂ t u (0 , t ) = ih [ A x,i ( t )] u (0 , t ) , u (0 ,
0) = I, (13)and ρ is the initial density matrix which characterizesthe initial occupation numbers n a of single-particle states | ψ a (cid:105) : ρ = (cid:88) a | ψ a (cid:105) n a (cid:104) ψ a | . (14)In our case, | ψ a (cid:105) are the eigenstates of the single-particleHamiltonian (7). If some occupation numbers are exactlyzero (which can be the case at zero temperature), somecomponents of the quantum evolution operator u com-pletely decouple and can be discarded in the solution ofthe equation (13). This can be used to speed up thealgorithm, typically by a factor of two (for a standardzero-temperature Fermi distribution). The expectationvalue (cid:104) j x,i ( t ) (cid:105) in (11) describes the effect of backreac-tion of fermions on the electromagnetic fields.We thus have a closed set of equations (11), (12) and(13), which allows to evolve the fermionic quantum statesand the classical electromagnetic fields in a self-consistentway. One can also explicitly check that this evolutionconserves the total energy H EM + (cid:104) ˆ H F (cid:105) of electromag-netic field and fermions up to the work done by the ex-ternal current J x,i ( t ): ∂ t (cid:16) H EM + (cid:104) ˆ H F (cid:105) (cid:17) = − (cid:88) x,i J x,i ( t ) E x,i ( t ) ,H EM = 12 (cid:88) x,i ( ∂ t A x,i ) + (cid:88) j F x,ij , (cid:104) ˆ H F (cid:105) = Tr (cid:0) ρ u (0 , t ) h u † (0 , t ) (cid:1) . (15)We have solved the evolution equations (11) and (13)using the leap-frog integrator, which slightly violates theconservation of energy (15). At sufficiently small timestep this violation is completely under numerical control,see Fig. 14 in the Appendix A.In the CSFT approach, one can also partially takeinto account the quantum fluctuations of the electromag-netic fields, encoded in the nontrivial Wigner transform¯ ρ EM ( A , E ) of the initial density matrix ρ EM , where A ≡ A ( t = 0) and E ≡ E ( t = 0) are the initial valuesof electric and magnetic fields. To this end one shouldadditionally average all observables over A and E , sam-pled with the probability ¯ ρ EM ( A , E )[58]. However, inour work we have not taken these initial quantum fluc-tuations into account for the following reasons. First, inthe case of anomaly equation (5) with massless fermionswe have found that the effect of initial fluctuations ofelectromagnetic fields is much more significant than inthe case of e.g. Schwinger pair production [38][59], andhence much more samples of the initial fields are requiredto reach acceptable statistical errors. Partially this canbe explained by the large value of the electromagneticcoupling constant, which was g = 1 . g , say, g = 0 .
1. Inthe latter case, however, the characteristic time scale ofthe chiral plasma instability increases significantly aboveour current simulation times.In addition, taking into account initial quantum fluctu-ations of electromagnetic fields makes it impossible to as-sume spatial homogeneity of electromagnetic fields alongsome of the lattice directions, which is essential to speedup the CSFT simulations at large lattice sizes.Thus while the effect of quantum fluctuations on thechiral plasma instability might be potentially very sig-nificant and interesting, we cannot study it with ourpresently available computational resources and leave itfor the future work. In this work, we avoid the statis-tical averaging over the initial values of the fields E , A by using the very simple form of the initial densitymatrix ¯ ρ EM ( A , E ) which is just a delta-function onsome particular, specifically chosen initial values. Wethus completely neglect the quantum fluctuations of theelectromagnetic fields. Nevertheless, this approximationis still certainly wider than the chiral kinetic theory orhydrodynamical approximation. III. CHIRALITY PUMPING IN PARALLELELECTRIC AND MAGNETIC FIELDS
In this Section we study the real-time evolution of theaxial charge Q A in the background of constant parallelexternal electric and magnetic fields. In the absence ofbackreaction, such setup provides a direct check of howwell the anomaly equation (5) holds for the Wilson-DiracHamiltonian with inexact chiral symmetry [44], which wefurther use to study the chiral plasma instability in Sec- tion IV. The effect of backreaction is also interesting sincethe anomaly equation (5) is known to receive nontriv-ial corrections if the electromagnetic fields are dynamical[47–49].In order to induce the constant external electric field (cid:126)E = E (cid:126)e , we switch on the external current of the form J x,i ( t ) = δ i, E t . Constant external magnetic field isinduced by the static circular external current flowingaround the plaquettes with x = L − x = L − x = 0 . . . L −
1, where L , L and L are latticesizes. This static current is like a thin solenoid piercing astack of lattice plaquettes, with the field strength beingequal to B outside of solenoid and B − BL L inside it.Lattice fermions, however, acquire only the Aharonov-Bohm phase e igB when encircling such a solenoid if oneimposes the flux quantization condition gBL L = 2 π Φ , Φ ∈ Z . (16)The total external current which we insert in the equa-tions (11) is the sum of the two currents which createconstant electric and magnetic fields.For the initial state of the fermionic fields, we usethe eigenstates of the Wilson-Dirac Hamiltonian h [ A ],where A is the initial gauge field configuration with con-stant magnetic field B as described above. In this workwe consider only the limit of zero temperature, corre-spondingly, only the states with negative energies areinitially occupied.The Wilson-Dirac Hamiltonian which we use in oursimulations does not have exact chiral symmetry, andthere is no uniquely defined axial charge operator whichwould exactly satisfy the anomaly equation (5) and com-mute with the Hamiltonian. Rather, the anomaly equa-tion (5) can only hold approximately, in the limit of largelattice volume and sufficiently smooth, slowly changingand small gauge fields [50, 51]. We thus use the simplestpossible definitions of the operators of the axial chargedensity q A x and the total axial charge Q A :ˆ q A x = ˆ ψ † x γ ˆ ψ x , ˆ Q A = (cid:88) x ˆ q A x . (17)The time-dependent expectation value of the axial chargedensity is calculated similarly to the expectation value ofelectric current in (12): (cid:104) q A x ( t ) (cid:105) = Tr (cid:0) ρ u (0 , t ) γ P x u † (0 , t ) (cid:1) , (18)where P x is the single-particle projector on a single latticesite x : [ P x ] x x = δ x x δ xx .First we neglect the backreaction of fermionic electriccurrent on the electromagnetic field and measure the timedependence of the axial charge in constant parallel exter-nal electric and magnetic fields. The results are shownon the left panel of Fig. 2 for the 10 × ×
32 lattice withΦ = 1 quantum of magnetic field flux. One can see that Q A grows linearly with time until it reaches some maxi-mal value Q A /V ≈ . V = L L L is the totalnumber of lattice sites (lattice volume). After that, Q A -2 0 2 4 6 8 0 20 40 60 80 100 120 140 Q A / V × t E=0.01E=0.02E=0.03E=0.04E=0.05 -16-14-12-10-8-6-4-2 0 0 0.01 0.02 0.03 0.04 0.05 0.06 α ( E ) × E FIG. 2: Chirality pumping without backreaction for various external electric fields and external magnetic field with flux Φ = 1on the 10 × ×
32 lattice. On the left: time dependence of the axial charge Q A and its linear fits at early times. On theright: dependence of the slope of these fits on the external electric field and the linear fit of this dependence. Here and in whatfollows time is expressed in units of lattice spacing. decreases again. This decrease is a lattice artifact relatedto the fact that the characteristic momentum p ∼ Et offermions accelerated by an electric field E approachesthe UV cutoff set by the compact size of the lattice mo-mentum space k µ ∈ [ − π . . . π ]. Due to the periodicity oflattice momentum space, at large time scales the behav-ior of the axial charge (in the absence of backreaction) iswell described by Q A ∼ sin ( Et/ Q A ( t ) /V = α ( E ) t in the range t ∈ [0 . . .
50] for E = 0 .
01 and E = 0 .
02 and in the range t ∈ [0 . . .
30] for other val-ues of E . The dependence of the coefficient α ( E ) onthe electric field is shown on the right panel of Fig. 2.Again, this dependence is linear with a good precision,and we perform another linear fit α ( E ) = CE , where C corresponds to the anomaly coefficient relating ∂ t Q A and (cid:82) d (cid:126)x (cid:126)E · (cid:126)B in (5). On Fig. 3 we show the dependence of C on the size of the lattice (in the directions perpendicularto the magnetic field). One can see how C approachesthe value C = π in the limit of large lattices, in agree-ment with the anomaly equation (5). Let us also notethat for larger number of flux quanta one can performa similar fitting procedure. However, finite-volume ar-tifacts in C are significantly larger for larger magneticfluxes. For this reason, in this work we have only usedexternal magnetic field with one flux quantum.It is also interesting to check how the axial charge de-pends on time after the external electric field is switchedoff and the external magnetic field remains constant(which we believe to be a more realistic experimentalsetup than the simultaneous switching off of all fields).The time dependence of the axial charge for such a sit-uation is shown on the left plot of Fig. 4 (green line). π C Lattice size
FIG. 3: Dependence of the lattice anomaly coefficient( ∂ t Q A = C (cid:82) d (cid:126)x (cid:126)E · (cid:126)B ) on the transverse lattice size. Lat-tice size in the direction parallel to the magnetic field is fixedto L = 32. External electric field is switched off at the time t = 50.One can see that starting from this moment of time theaxial charge exhibits only some small-scale fluctuationsaround the nonzero mean value. This demonstrates thatthe effect of explicit chiral symmetry breaking due to theWilson term in the Hamiltonian (7) is rather small forsuch simulation parameters, and the total axial charge isalmost a conserved quantity.After establishing the validity of the anomaly equa-tion (5) in our simulation setup, we study the effect ofbackreaction of dynamical electromagnetic fields on thechirality pumping process. Technically, the backreactionis taken into account by inserting the expectation valueof the fermionic electric current (cid:104) j x,i (cid:105) into the Maxwellequations for the electromagnetic field. We now consider -0.5 0 0.5 1 1.5 2 2.5 0 50 100 150 200 Q A / V × t no backreaction, E field is switched off at t=50with backreactionno backreaction -4-2 0 2 4 6 8 10 0 50 100 150 200 E × t FIG. 4: A comparison of chirality pumping processes with and without backreaction of the fermionic electric current on theelectromagnetic field. On the left: time dependence of the axial charge. On the right: time dependence of the component ofthe volume-averaged electric field parallel to the magnetic field. Lattice size is 10 × ×
32, the flux of external magnetic fieldis Φ = 1, external electric field is 0 . the situation in which the external electric and magneticfields are switched on permanently. As we will see, insimulations with backreaction switching off the electricfield at sufficiently late times anyway does not affect theevolution significantly due to screening by the dynam-ically generated electric field. Time dependence of theaxial charge Q A ( t ) for simulation with backreaction isshown on Fig. 4. For comparison, on the same Figure wealso show Q A ( t ) for simulations without backreaction,where the electric field is permanent or switched off at t = 50.One can see that while at t (cid:46) Q A ( t ) grows approx-imately linearly with t both with and without backreac-tion, at later times backreaction leads to a rapid decay of Q A with subsequent fluctuations around zero. In orderto understand the origin of this effect, remember that ax-ial anomaly can be also regarded as the Schwinger pairproduction in the effective 1 + 1-dimensional theory offermions on the lowest Landau level. It is thus natural toexpect that particle-anti-particle pairs produced by theexternal electric field will tend to screen this field, justas in the case of Schwinger effect in (3 + 1) dimensions.To check this conjecture, on the right panel of Fig. 4 weplot the volume-averaged electric field projected on thedirection of the magnetic field. One can see that indeedit quite quickly decreases from the initial value E = 0 . t ≈
30 - exactly at the time atwhich the growth of the axial charge stops (see left panelof the same Figure). After that the electric field exhibitssome fluctuations around zero with the amplitude whichis approximately five times smaller than the initial fieldvalue. We thus conclude that the effect of backreactionon the chirality pumping is to stop the growth of the ax-ial charge by screening the external electric field down tozero.
IV. CHIRAL PLASMA INSTABILITY ANDDECAY OF AXIAL CHARGE
In this Section, we consider a situation in which someinitial chiral imbalance is already created e.g. by chiral-ity pumping, and the parallel electric and magnetic fieldsare adiabatically switched off while keeping nonzero valueof the total axial charge Q A and hence the chiral chem-ical potential µ A . In this setup we would like to studythe existence and the late-time evolution of the exponen-tially growing solutions (4) of the anomalous Maxwellequations (2), as well as the associated inverse cascade ofenergy of helical electromagnetic fields.In order to implement the initial chirality imbalance asdiscussed in the introductory Section I (see right panelof Fig. 1), we divide all the eigenstates of the single-particle Wilson-Dirac Hamiltonian h [ A ], where A isthe initial value of the vector potential, into the positivechirality states with (cid:104) ψ a | γ | ψ a (cid:105) > (cid:104) ψ a | γ | ψ a (cid:105) <
0. For positive chiralitystates we fill all the levels with (cid:15) a < µ A , and for negativechirality states - all the levels with (cid:15) a < − µ A . Whilethe eigenstates of the Wilson-Dirac Hamiltonian are notin general the eigenstates of the γ operator, for eigen-states with sufficiently small momenta our definition ismaximally close to the notion of distinct Fermi levels ofleft- and right-handed continuum massless fermions. Inpractice, this definition is unambiguous as long as all theenergy levels (cid:15) a are non-degenerate. In the case of de-generate energy levels (as e.g. in the case of the Wilson-Dirac Hamiltonian with zero gauge fields or in the back-ground of a single plane wave), one can additionally ro-tate the eigenstates within the degenerate subspaces inorder to maximize the absolute values of matrix elements (cid:104) ψ a | γ | ψ a (cid:105) .With lattice discretizations of the Dirac Hamiltonianit is not possible to have very large values of the chiralchemical potential µ A , since the dispersion relation atthe Fermi energy µ A (cid:38) (cid:15) (cid:16) (cid:126)k (cid:17) = v F | (cid:126)k | due to lattice artifacts. At µ A = 2,the Fermi energy touches the lowest van Hove singularity(saddle point) of the dispersion relation, and the excita-tions around the Fermi surface no longer correspond toDirac fermions. On the other hand, the solution (4) ofthe anomalous Maxwell equations (2) with the conven-tional value σ CME = µ A π of the chiral magnetic conduc-tivity suggests that the wave vectors at which the chiralplasma instability can occur are bounded by | (cid:126)k | < µ A π .On a finite spatial lattice of size L with periodic boundaryconditions, the smallest nonzero value of | (cid:126)k | is | (cid:126)k | = πL ,which dictates the lower bound on the size of the latticewhere the instability can be observed: L > π µ A . (19)Thus it is advantageous to use large values of µ A in or-der to reduce the lattice size used for simulations. Takingthe moderate value µ A = 0 .
75, at which the dispersionrelation is still linear with a good precision, we obtain
L > L in the direction x of electromagnetic wavepropagation is much larger than the sizes L = L ≡ L s in the transverse directions x and x . In addition, wehave assumed that electromagnetic fields do not dependon the transverse coordinates x and x . This allows usto represent the single-particle evolution operator u (0 , t )in the block-diagonal form in the basis of plane wavespropagating along x and x , which greatly reduces thedimensionality of the linear space on which the single-particle Schr¨odinger equation (13) should be solved. Bycomparing the results of simulations with L s = 20 and L s = 40 at fixed L = 200 (see Table I and Figs. 5,6, 10, 11 and 12) we have checked that the dependenceon the transverse lattice size is rather weak. Let us alsonote that one of the reasons for not using the final state ofchirality pumping process described in Section III for thestudy of chiral plasma instability is that in this case it isnot possible to assume spatial homogeneity in transversedirections due to the breaking of translational invarianceby external magnetic field [52].While the fermionic initial state described above is anexcited state which can spontaneously decay due to chiralplasma instability, in numerical simulations one alwaysneeds some small “seed” perturbation to start the decayprocess in a controllable way. For this reason we havestarted our simulations with a state in which also some fi-nite number n of electromagnetic field modes are excited.All of them are plane waves propagating along the latticedirection x with the largest size L , with a few small-est nonzero wave numbers k m = πmL , m = 1 . . . n and Set.No. L L s µ A n f v F Q A ↓ I Bk ↑ I Ek ↑ (cid:51) (cid:55) (cid:51) (cid:51) (cid:55) (cid:51) (cid:51) (cid:51) (cid:51) (cid:55) (cid:51) (cid:55) (cid:51) (cid:51) ? (cid:55) (cid:51) ? (cid:55) (cid:51) ? (cid:51) (cid:55) (cid:51) (cid:51) (cid:55) (cid:55) TABLE I: Summary of parameters and results of our simu-lations of chiral plasma instability. The column Q A ↓ sum-marizes the decay of the axial charge and the columns I E,Bk ↑ summarize the growth of the energies of long-wavelength elec-tric and magnetic fields. The symbols (cid:51) , (cid:55) and ? denote,respectively, the clearly visible growth, clearly visible absenceof growth and intermediate situations for which it is difficultto make any conclusion within a finite simulation time. random linear polarizations. In order to facilitate the de-tection of the inverse cascade, we choose the amplitudesof all modes in such a way that their contributions to thetotal energy of electromagnetic field are equal. Thus theexplicit form of our initial electromagnetic field configu-ration is A x,i ( t = 0) = n (cid:88) m =1 fw ( k m ) n m i cos ( k m x + φ m ) ,E x,i ( t = 0) ≡ ∂ t A x,i ( t ) | t =0 == n (cid:88) m =1 f n m i sin ( k m x + φ m ) , (20)where n m i are the random unit transverse polarizationvectors which are chosen to coincide with one of the basisvectors (cid:126)e , (cid:126)e with equal probability, φ m ∈ [0 , π ] are therandom phases and w ( k m ) = (cid:113) (cid:0) k m (cid:1) correspondsto the lattice dispersion relation for free massless fieldson the lattice.In order to understand how the evolution process de-pends on various lattice parameters, we have performedsimulations with 9 different parameter sets, which aresummarized in Table I. We have varied both the trans-verse and the longitudinal lattice sizes, the initial elec-tromagnetic field amplitude f , the number n of initiallyexcited electromagnetic field modes, the initial value ofaxial charge and the Fermi velocity v F . The parameterset No. 1 with L = 200, L s = 20, µ A = 0 . n = 10, f = 0 . v F = 1 is the “default” parameter set, andall other sets differ from it by a change in a few pa-rameters. Correspondingly, in what follows we label thedata points on the plots which combine the results fromseveral simulations by the number of parameter set (pre-ceded by the hash symbol -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0.1 1 10 100 1000 10000 100000 Q A / V t s =40) F =0.75) =20,n=1, µ A =1.0) 3.5 4 4.5 5 5.5 6 6.5 0.1 1 10 100 1000 10000 Q A / V t µ A =1.5) µ A =1.5,f=0.05) FIG. 5: Time dependence of the total axial charge for different simulation parameters. On the right plot we compare simulationswith different µ A and linearly rescale q A ( t ) → cq A ( t ) so that the initial values of cq A (0) agree for all simulations. In the plotlabels, numbers preceded by the hash symbol L = 200, L s = 20, µ A = 0 . n = 10, f = 0 . v F = 1, parameter set No. 1). On Fig. 5 we show the time dependence of axial chargein simulations with parameters summarized in Table I.In simulations with initial amplitude of electromagneticfields being equal to f = 0 . Q A de-cays with time. Interestingly, simulations with the small-est lattice size (parameter set No. 9) exhibit the fastestdecay of Q A . On the other hand, with the initial am-plitude f = 0 .
05 the axial charge density exhibits onlya rather small decrease at intermediate evolution times,subsequently followed by a slight increase. This non-trivial dependence on the electromagnetic field strengthsuggests that the dynamics of the decay process is morecomplicated than suggested by the anomalous Maxwellequations (2). It is interesting that the evolution of theaxial charge seems to depend only weakly on simulationparameters other than the initial amplitude f and thelongitudinal lattice size L (through the value of the low-est wave number πL ). Even the dependence on the chiralchemical potential µ A appears to be rather weak (aftera trivial rescaling with respect to the initial value), seeright plot on Fig. 5. The characteristic time scale forthe evolution of the axial charge appear to be essentiallylarger than in the chirality pumping simulations in theprevious Section. This difference can be qualitativelyexplained by much weaker field strengths in the simula-tions described in this Section. We also note that theinitial values of the axial charge are roughly consistentwith the continuum formula Q A /V = µ A / (cid:0) π (cid:1) , where V is the lattice volume. Deviations from this value canbe explained, first, by the inclusion of the initial vectorpotential A in the initial Hamiltonian, and second, bythe smaller value of chirality |(cid:104) ψ a | γ | ψ a (cid:105)| < | ψ a (cid:105) of the Wilson-Dirac Hamiltonian(7).A scaling analysis of the anomalous Maxwell equationssuggests that at late evolution times the time dependence Q A ( ) / Q A ( t ) t/10 s =40) µ A =1.5) F =0.75) =20,n=1, µ A =1.0) FIG. 6: Universal late-time scaling Q A ( t ) ∼ / √ t of axialcharge density in simulations with f = 0 .
2. Time depen-dence of ( Q A (0) /Q A ( t )) in the second half of the total evo-lution time is fitted, where appropriate, by linear functions( Q A (0) /Q A ( t )) = A + Bt . The fits are shown with dashedblack lines. of the axial charge density approaches the simple powerlaw [13, 46] Q A ( t ) ∼ / √ t. (21)In order to check this scaling, on Fig. 6 we plot thetime dependence of the inverse square of the axial charge,which should approach the linear function 1 /Q A ( t ) ∼ t according to (21). This asymptotic behavior indeedseems to emerge at late evolution times for simulationswith L s = 20, n = 10, f = 0 . µ A = 0 .
75, both with v F = 1 and v F = 0 .
75 (parameter sets No. 1 and 8). Thelinear fits of ( Q A (0) /Q A ( t )) for these simulations areshown on Fig. 6 with dashed black lines.0We now check whether the decay of the axial charge isaccompanied by the growth of the long-wavelength modesof electromagnetic field, as predicted by the anomalousMaxwell equations (2). To this end we perform theFourier transforms of the transverse electric and mag-netic fields (taking into account that they depend onlyon the x coordinate) E k,i ( t ) = 1 √ L (cid:88) x e ikx E x,i ( t ) ,B k,i ( t ) = 1 √ L (cid:88) x e ikx B x,i ( t ) , (22)where i = 1 ,
2, and further decompose the Fourier-transformed fields into the helical components E k,R/L ( t )and B k,R/L with right- and left-handed helicities: B k,R ( t ) = 12 ( B k, ( t ) + B − k, ( t )) ++ 12 i ( B k, ( t ) − B − k, ( t )) ,B k,L ( t ) = 12 i ( B k, ( t ) − B − k, ( t )) ++ 12 ( B k, ( t ) + B − k, ( t )) . (23)For electric fields, the definition of helical componentsis exactly the same. Again, here the term “helicity”refers to the direction of rotation of transverse electricand magnetic fields along the spatial direction of wavepropagation (the x axis in our setup).After such a decomposition, we calculate the energiesof left- and right-handed helical electric and magneticfields with a given wave number k as I Bk,R/L ( t ) = (cid:12)(cid:12) B k,R/L ( t ) (cid:12)(cid:12) / (cid:12)(cid:12) B − k,R/L ( t ) (cid:12)(cid:12) / ,I Ek,R/L ( t ) = (cid:12)(cid:12) E k,R/L ( t ) (cid:12)(cid:12) / (cid:12)(cid:12) E − k,R/L ( t ) (cid:12)(cid:12) / , (24)where k = πmL and m = 0 . . . (cid:98) L / (cid:99) now spans only thehalf of the discrete lattice momenta. Since the initial con-figuration of electromagnetic fields contains plane waveswith (random) linear polarizations and equal energies, at t = 0 the energies I k,R/L ( t ) of all left- and right-handedelectromagnetic modes with 0 < k ≤ πnL are equal.We have found that for all simulations the energies I E,Bk,R/L ( t ) exhibit quite large short-scale fluctuations withperiod of order ∆ t ∼ . . . I Ek,R ( t ) within a short initial period of timefor simulation with parameter set No. 1. These oscil-lations indicate that the helical magnetic and electricfields represented by the basis (23) are not the eigen-states of the evolution process, which is in sharp contrastto the solution (4) of the anomalous Maxwell equations(2). We have explicitly checked that if the backreactionof fermions on the electromagnetic fields is neglected, FIG. 7: Time dependence of the energies I Ek,R/L ( t ) of right-handed components of electric field on a short time interval atthe beginning of evolution for parameter set No. 1 ( L s = 200, n = 10, µ A = 0 . f = 0 . k = πL (largest wavelength) to pure red for k = πnL . these oscillations disappear and the energies I E,Bk,R/L ( t )are constant in time for all values of k and for all polar-izations. This observation suggests that the short-scaleoscillations might originate from the nontrivial depen-dence of fermionic current on the frequency, wave numberand amplitude of electromagnetic field, which turns thesolutions of the anomalous Maxwell equations (2) intowaves with generic elliptic polarizations.Despite the short-scale fluctuations, we still find it use-ful to decompose our fields in the basis (23), since thecorresponding electromagnetic modes carry definite he-licity and thus the energies of helical modes can be usedto define, at least approximately, the helicity on the lat-tice. This definition is advantageous since direct latticediscretizations of the continuum formula H ∼ (cid:82) d x (cid:126)A · (cid:126)B are in general flawed by lattice artifacts. In order to ab-stract ourselves from the short-scale fluctuations, we de-fine the energies ¯ I k,R/L ( t ) which are averaged over somefinite time interval T :¯ I E,Bk,R/L ( t ) = 1 T t + T/ (cid:90) t − T/ dt (cid:48) I E,Bk,R/L ( t (cid:48) ) . (25)We have used the value T = 25, which is sufficient toremove practically all short-scale oscillations.On Figures 8 and 9 we separately illustrate the timedependence of the energies of the left-handed (on the left)and right-handed (on the right) helical magnetic and elec-tric fields with wave numbers k ≤ πnL for several mostcharacteristic sets of simulation parameters. The wavenumbers are coded in color, from pure blue for the small-est nonzero value k = πL (largest wavelength) to pure redfor k = πnL . Semi-transparent colored regions show therange of short-scale oscillations of I E,Bk,R/L ( t ), and thicksolid lines show the time dependence of the time-smeared1 FIG. 8: Time dependence of the energies of helical magnetic fields for several selected sets of simulation parameters. Thewave numbers are coded in color, from pure blue for the smallest nonzero value k = πL (largest wavelength) to pure red for k = πnL . Semi-transparent colored regions show the range of short-scale oscillations of I Bk,R/L ( t ), and thick solid lines showthe time dependence of the time-smeared energies ¯ I Bk,R/L ( t ) defined in (25). In black-and-white version pure blue and pure redcorrespond to black and light grey, respectively. Horizontal green (light grey) lines show the initial energies which are equal forall modes. Left-handed and right-handed modes are in the left and in the right columns, respectively. FIG. 9: Time dependence of the energies of helical electric fields for several selected sets of simulation parameters. Thewave numbers are coded in color, from pure blue for the smallest nonzero value k = πL (largest wavelength) to pure red for k = πnL . Semi-transparent colored regions show the range of short-scale oscillations of I Ek,R/L ( t ), and thick solid lines showthe time dependence of the time-smeared energies ¯ I Ek,R/L ( t ) defined in (25). In black-and-white version pure blue and pure redcorrespond to black and light grey, respectively. Horizontal green (light grey) lines show the initial energies which are equal forall modes. Left-handed and right-handed modes are in the left and in the right columns, respectively. I k,R/L ( t ) defined in (25). Horizontal green linesshow the initial energies which are equal for all helicalcomponents of electric and magnetic fields.From Fig. 8 we see that in some simulations (param-eter sets No. 3 – 7) the energies of the helical compo-nents of magnetic field exhibit the expected signaturesof the inverse cascade due to chiral plasma instability[3, 6, 13, 46]. Namely, the energy of a single longest-wavelength right-handed helical mode rapidly grows atearly times and reaches some saturation limit at lateevolution time, whereas the energies of all other modesdecrease with time. As expected from the anomalousMaxwell equations (2) with σ CME = µ A / (cid:0) π (cid:1) , increas-ing µ A by a factor of two (to µ A = 1 .
5, parameter setNo. 3, second row in Fig. 8) results in the growth oftwo right-handed modes. Comparing the data on Fig. 8and Fig. 5, we conclude that the growth of helical mag-netic fields is not necessarily accompanied by the decayof the axial charge, and vice versa (see also Table I fora summary of all simulations). Yet another observationwhich supports this conclusion is that in simulations onthe smallest lattice, for which the axial charge exhibitsmost rapid decay, we have not found any signatures ofthe growing electromagnetic fields. Interestingly, increas-ing the value of µ A and/or the number of initially ex-cited electromagnetic field modes also does not necessar-ily speed up the inverse cascade and the decay of Q A .An even more interesting picture emerges if we alsoconsider the energies of the helical components of elec-tric fields, shown on Fig. 9. It turns out that for somesimulation parameters the long-wavelength helical com-ponents of the electric field, rather than the magneticfield, are enhanced during the evolution (parameter setsNo. 1, 2, 8). For parameter set No. 3, both magnetic andelectric fields grow in time. It is remarkable that pre-cisely for these parameter sets the axial charge exhibitsmost rapid decay. It seems that both the growth of thehelical electric fields and the decay of the axial chargeare triggered by sufficiently large initial amplitudes ofelectromagnetic fields. Thus it seems that the roles ofelectric and magnetic fields in the chiral plasma insta-bility scenario are essentially different, in contrast to thesimple solution (4) of the anomalous Maxwell equations.It is also interesting to note that we observe the maximalgrowth of long-wavelength helical electric fields in simu-lations with a smaller value of Fermi velocity v F = 0 . ξ B ( t ) and ξ E ( t ) as ξ E,B ( t ) = (cid:80) k πk I E,Bk ( t ) (cid:80) k I E,Bk ( t ) , (26)where I E,Bk ( t ) = I E,Bk,R ( t ) + I E,Bk,L ( t ). The time depen-dence of ξ E ( t ) and ξ B ( t ), shown on Fig. 10, quantifies the direction of the transfer of energy between short- andlong-wavelength modes. ξ B ( t ) and ξ E ( t ) can be alsothought of as the average wavelengths of magnetic andelectric fields at a given moment of time. The data shownon Fig. 10 indicates that ξ E and ξ B on average increasewith time practically for all our simulations, thus provid-ing a more quantitative evidence for the inverse cascade.The growth is somewhat more pronounced for the elec-tric correlation length ξ E , especially in simulations withlarger initial amplitude f = 0 .
2. At late evolution times, ξ E saturates at its upper bound ξ E = L equal to the lat-tice size. In contrast, the magnetic correlation length ξ B exhibits rapid growth only at early times, and later seemsto saturate at values smaller than L . On Fig. 10 we donot show the data for the parameter set No. 9, since inthis case we have found that only a single initially ex-cited mode strongly dominates the spectrum throughoutthe whole evolution process, and the quantities ξ E and ξ B are trivially equal to L up to some very small cor-rections.The effect of saturation of ξ E,B at late evolution timesprevents us from checking the universal late-time behav-ior ξ E,B ∼ √ t which follows from the scaling analysis ofthe anomalous Maxwell equations [13, 46], similarly to(21). It seems, however, that the late-time behavior of ξ E is more similar to √ t than that of ξ B .So far almost all theoretical studies of the chiral plasmainstability assume that the axial charge is distributedhomogeneously in space and can be described by acoordinate-independent chiral chemical potential µ A atall evolution times. The extension of the anomalousMaxwell equations (2) which allows to consider spatiallyinhomogeneous distributions of axial charge density hasbeen constructed only recently in [27]. It is thus inter-esting to check how well the assumption of spatial ho-mogeneity of the axial charge density q A x holds in oursimulations. In order to quantify the spatial inhomo-geneity of q A x , we consider the space-averaged squareddeviation of q A x from its space-averaged value Q A /V : σ [ q A ] = (cid:115)(cid:88) x ( q A x − Q A /V ) . (27)The time dependence of the ratio σ [ q A ] / ( Q A /V ) isshown on Fig. 11 for those sets of simulation parame-ters which exhibit axial charge decay (in particular, thisfixes f = 0 . A of vec-tor potential, even at the start of the evolution the axialcharge density is slightly inhomogeneous, with deviationsfrom mean value being of order of 5 %. As one can seefrom Fig. 11, at late evolution times the inhomogeneityof q x A tends to slightly increase, however, this increaseis not dramatic and does not exceed 20%. This suggeststhat the approximation of spatially homogeneous axialcharge distribution is not unreasonable even when thelong-wavelength modes are strongly enhanced and domi-nate the evolution. For simulations with f = 0 .
05 which4 ξ B / L z t/10 s =40) µ A =1.5) µ A =1.5,f=0.05) F =0.75) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 ξ E / L z t/10 s =40) µ A =1.5) µ A =1.5,f=0.05) F =0.75) FIG. 10: Magnetic and electric correlation lengths ξ B and ξ E (defined in (26)) as functions of time for different parameter sets.A smearing procedure similar to (25) was applied in order to suppress minor short-scale fluctuations in the data. σ [ q A ]/ ( Q A / V ) t/10 =20,n=1, µ A =1.0) s =40) µ A =1.5) F =0.75) FIG. 11: Time dependence of the standard deviation of theaxial charge density σ [ q A ] from its volume-averaged value Q A /V in simulations with f = 0 . do not exhibit the decay of the axial charge, the inho-mogeneities of the axial charge density remain approxi-mately constant or even tend to decrease.An interesting question is also the net transfer of en-ergy between fermions and electromagnetic fields. As dis-cussed in Section II, in classical statistical field theory al-gorithm the total energy of fermions and electromagneticfields is conserved up to the work performed by the ex-ternal current (see Fig. 14 in Appendix A for a numericaldemonstration). Since in the simulations of chiral plasmainstability discussed in this Section the external currentsare absent, the transfer of energy can be characterizedby the time dependence of the energy of electromagneticfield alone, which is illustrated on Fig. 12. We see that inalmost all simulations the energy of electromagnetic fielddecreases or stays constant. The only exception is thesimulation with parameter set No. 7 ( n = 10, f = 0 . µ A = 1 . t (cid:38) · .Analysis of power spectra suggests that this increase canbe at least partly attributed to the enhancement of helicallong-wavelength electric fields (similarly to the one ob-served for parameter set No. 3, see second row on Fig. 9).The decrease of electromagnetic field energy in all othersimulations indicates that it might be not completely cor-rect to think of chiral instability as of a “discharge” ofon excited state of Dirac sea into electromagnetic waves. V. CONCLUSIONS
In this work, we have studied the real-time quan-tum evolution of chirally imbalanced Wilson-Dirac latticefermions coupled to dynamical classical electromagneticfield within the classical statistical field theory approach.The quantum evolution of fermions was simulated ex-actly (up to small fully controlled errors originating fromdiscretization of time). Our simulations of the chiral-ity pumping process, described by the volume-integratedanomaly equation (5), suggest that the effect of explicitchiral symmetry breaking due to the Wilson term in thelattice Dirac Hamiltonian is not very large. We hopetherefore that our results can be confronted at least atthe qualitative level with the theoretical predictions forcontinuum chiral fermions.We have considered both the generation of chirality im-balance in parallel electric and magnetic fields and the de-cay of initially present chirality imbalance at the expenseof generating electromagnetic fields with nonzero helicity.We have observed that in general the backreaction of dy-namical electromagnetic fields prevents fermions from ac-quiring large chirality imbalance - either by suppressionof the chirality pumping or by accelerating the decay ofinitially present chirality imbalance. The suppression ofthe chirality pumping process can be understood as thedynamical screening of the external electric field, sim-ilarly to what happens in the Schwinger pair creation5 H E M t/10 s =40) µ A =1.5) F =0.75) =20,n=1, µ A =1.0) 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 H E M t/10 µ A =1.5,f=0.05) FIG. 12: Time dependence of the total energy of electromagnetic field in simulations with f = 0 . f = 0 . process [36, 38].In simulations with nonzero initial axial charge Q A wehave also found a numerical evidence of the inverse cas-cade phenomenon due to the chiral plasma instability -that is, rapid growth of long-wavelength magnetic fieldsof definite helicity at early evolution times and the decayof all other magnetic field components. In some cases, he-lical electric fields were found to grow, even when mag-netic fields did not exhibit any enhancement. A sum-mary of our simulations given in Table I suggests thatthe growth (or at least the absence of decay) of long-wavelength helical electric fields is a necessary conditionfor the dynamical decay of the axial charge. The factthat the enhancement of helical electric fields is switchedon only for sufficiently large initial amplitude of electro-magnetic field indicates that nonlinear responses such asthe dynamical refringence [53] might be important forthe evolution of chirally imbalanced plasma.We have observed the mechanism which eventuallystops the growth of long-wavelength modes in our sim-ulations is not directly related to the decay of the ax-ial charge. This observation, together with quite differ-ent roles of electric and magnetic fields in the evolutionprocess, suggests that the nontrivial momentum and fre-quency dependence of both the electric conductivity andthe chiral magnetic conductivity might be important forthe quantitative description of chiral plasma instability.On the other hand, our simulations also indicate that theapproximation of spatially homogeneous axial charge dis-tribution, assumed in most theoretical considerations ofanomalous Maxwell equations, is reasonably good evenat late evolution times, when the instability has fully de-veloped and the growth of long-wavelength helical elec-tromagnetic fields has saturated.An interesting further development of our work wouldbe to use chiral lattice fermions in the CSFT algorithm,with the possible choice of overlap Hamiltonian [54]. Inthis case, axial charge is conserved in the absence of elec-tromagnetic fields, and the effects of explicit chiral sym- metry breaking at high momenta should be absent. Suchsetup should be more relevant in the context of high-energy physics, where chiral symmetry tends to be exactat sufficiently high energies (at least at the level of thebare Lagrangian). Yet another interesting open questionis the effect of the quantum fluctuations of the electro-magnetic field, which are encoded in the nontrivial initialdensity matrix. Acknowledgments
We thank D. Kharzeev, A. Sadofyev and N. Yamamotofor interesting and stimulating discussions of the physicsof chiral media, and F. Hebenstreit and D. Gelfand foruseful discussions of the CSFT algorithm. This workwas supported by the S. Kowalevskaja award from theAlexander von Humboldt Foundation. We are also in-debted to S. Valgushev for his help with code paralleliza-tion at the initial stage of this work as well as for thecareful reading of this manuscript. The authors are grate-ful to FAIR-ITEP supercomputer center where a part ofthese numerical calculations was performed.
Appendix A: Classical statistical field theoryalgorithm
The starting point of our derivation of the CSFT al-gorithm is the general expression for the time-dependentexpectation value of some quantum operator ˆ O : (cid:104) ˆ O ( t ) (cid:105) = Tr (cid:16) ˆ ρ ˆ U ( t , t ) ˆ O ˆ U † ( t , t ) (cid:17) , (A1)where ˆ ρ is the initial density matrix and the evolutionoperator ˆ U ( t , t ) is the time-ordered exponentˆ U ( t , t ) = T exp − i t (cid:90) t dt (cid:48) ˆ H ( t (cid:48) ) , (A2)6where the Planck constant is set to one by an appropri-ate choice of units. The Hamiltonian operator is definedby equations (6), (7) and (9) in Section II. We have al-lowed for an explicit time dependence of the Hamiltonian,for example, due to the time dependence of the externalcurrent J x,i ( t ). The evolution operator (A2) can be ex-panded into a product of elementary evolution operatorsfor small time step δ = t − t N , where δ − should be muchlarger than any relevant energy scale in the system:ˆ U ( t, t ) == lim δ → (cid:16) e − i ˆ H ( t ) δ e − i ˆ H ( t + δ ) δ . . . e − i ˆ H ( t ) δ (cid:17) . (A3)Let us now insert the decompositions of the identity op-erator ˆ I = (cid:81) x,i (cid:82) dA x,i | A x,i (cid:105)(cid:104) A x,i | in the Hilbert space ofelectromagnetic field between the infinitesimal factors aswell as at the beginning and at the end of the product in(A3) in order to arrive at the path integral representationof the evolution operator (A2). We do this both for theforward and the backward evolution operators ˆ U ( t , t )and ˆ U † ( t , t ) in (A1). It is convenient to enumeratethe gauge fields which enter identity decompositions inthe forward evolution operators with the discrete latticetime variable τ = 0 . . . N , and in the backward branch -with τ = N + 1 . . . N + 1. The variable τ is a discreteparametrization of the Keldysh contour going from t to t and back (see Fig. 13 for an illustration). Now we have toexpress the matrix elements (cid:104) A τ | e ∓ i ˆ H ( τ ) δ | A τ +1 (cid:105) of the elementary evolution operators in terms of the fields A τ and A τ +1 . In the derivation of the CSFT algorithm, it ismost convenient to use the approximate expression (cid:104) A τ | e − i ˆ H ( τ ) δ | A τ +1 (cid:105) ≈ e + i δ (cid:80) x,i ( A τ +1 x,i − A τx,i ) ×× e − iδ (cid:80) x,i,j ( F τx,i,j ) − iδ ˆ H F [ A τ ] − iδ (cid:80) x,i A τx,i J τx,i (A4)for the forward evolution operators, and the different ap-proximate expression (cid:104) A τ | e + i ˆ H ( τ ) δ | A τ +1 (cid:105) ≈ e − i δ (cid:80) x,i ( A τ +1 x,i − A τx,i ) ×× e + iδ (cid:80) x,i,j ( F τ +1 x,ij ) + iδ ˆ H F [ A τ +1 ] + iδ (cid:80) x,i A τ +1 x,i J τ +1 x,i (A5)for the backward evolution operators. In the first expres-sion (A4), we order the electromagnetic field operatorsˆ E x,i and ˆ A x,i in such a way that all the operators inthe exponent containing ˆ A x,i act on the vector (cid:104) A τ | . Inthe second expression (A5), these operators act on thevector | A τ +1 (cid:105) . These approximations are both valid toorder O ( δ ) and differ only in the terms of order O ( δ ),hence being equivalent in the limit δ → (cid:104) ˆ O ( t ) (cid:105) , in whichwe integrate over the gauge fields living on the discretizedKeldysh contour: (cid:104) ˆ O ( t ) (cid:105) = (cid:90) dA . . . dA N +1 ρ EM (cid:2) A , A N +1 (cid:3) ×× Tr (cid:32) ˆ ρ F e + i δ (cid:80) x,i ( A x,i − A x,i ) e − iδ ˆ H F [ A ] e − iδ (cid:80) x,i,j ( F x,ij ) − iδ (cid:80) x,i A x,i J x,i × . . .. . . × e + i δ (cid:80) x,i ( A Nx,i − A N − x,i ) e − iδ ˆ H F [ A N − ] e − iδ (cid:80) x,i,j ( F N − x,ij ) − iδ (cid:80) x,i A N − x,i J N − x,i ˆ O (cid:2) A N , A N +1 (cid:3) ×× e − i δ (cid:80) x,i ( A N +2 x,i − A N +1 x,i ) e iδ ˆ H F [ A N +2 ] e iδ (cid:80) x,i,j ( F N +2 x,ij ) + iδ (cid:80) x,i A N +2 x,i J N − x,i × . . .. . . × e − i δ (cid:80) x,i ( A N +1 x,i − A Nx,i ) e iδ ˆ H F [ A N +1 ] e iδ (cid:80) x,i,j ( F N +1 x,ij ) + iδ (cid:80) x,i A N +1 x,i J x,i (cid:33) . (A6)It is important to stress that at this point we have usedthe path integral representation only for the bosonicfields, and the exponential factors e ± iδ ˆ H F [ A τ ] in (A6)are still operators in the fermionic many-body Hilbertspace. Correspondingly, the trace in (A6) is takenover this Hilbert space. The operator of the observ-able ˆ O (cid:2) A N , A N +1 (cid:3) = (cid:104) A N | ˆ O | A N +1 (cid:105) is also an oper-ator on the fermionic Hilbert space which depends onfields A N and A N +1 . In deriving the above path inte-gral representation, we have made a simplifying assump- FIG. 13: An illustration of the Schwinger-Keldysh contourwith the discrete lattice time τ . ρ factorizes intothe direct product of the fermionic density matrix ˆ ρ F (which might in general be correlated with the initialstate of the electromagnetic field) and the density matrixˆ ρ EM of the electromagnetic field with matrix elements ρ EM (cid:0) A , A N +1 (cid:1) = (cid:104) A | ˆ ρ EM | A N +1 (cid:105) . While this as-sumption is certainly not valid, say, for the density ma-trix ˆ ρ = e − ˆ H/T describing the thermal equilibrium stateof the full Hamiltonian ˆ H EM + ˆ H F , it is still justifiable inthe case of almost classical dynamics of electromagneticfields.At this point let us assume that the observable oper-ator ˆ O (cid:2) A N , A N +1 (cid:3) can be represented as a sum of theidentity operator in fermionic Hilbert space (this sum-mand corresponds to purely bosonic observables) and ofall possible fermionic bilinear operators:ˆ O (cid:0) A N , A N +1 (cid:1) = O B (cid:0) A N , A N +1 (cid:1) ˆ I ++ (cid:88) x,y ˆ ψ † x (cid:2) O F (cid:0) A N , A N +1 (cid:1)(cid:3) x,y ˆ ψ y . (A7)This form is sufficiently general to describe all the ob-servables which we consider in this work. Furthermore,let us assume that the fermionic density matrix ˆ ρ F canbe represented as an exponent of some fermionic bilinearoperator ˆ H = (cid:80) x,y ˆ ψ † x [ h ] x,y ˆ ψ y :ˆ ρ F = Z − exp (cid:16) − ˆ H /T (cid:17) , (A8)where T is some (perhaps fictitious) temperature. Notethat in the case of evolution which starts from non-equilibrium state the operator ˆ H can be different fromthe Wilson-Dirac Hamiltonian ˆ H F which governs thequantum evolution. For instance, the excited state withinitial chiral imbalance considered in Section IV corre-sponds to the following form of h in the limit T → h = h (cid:2) A (cid:3) + µ A (cid:88) a | ψ a (cid:105) sign ( (cid:104) ψ a | γ | ψ a (cid:105) ) (cid:104) ψ a | , (A9)where | ψ a (cid:105) are the eigenstates of the Wilson-DiracHamiltonian (7) with the initial gauge field A . It is obvious that for exactly chiral Dirac Hamiltonian with (cid:104) ψ a | γ | ψ (cid:105) a = ± h = h + µ A γ .Now we are in the position to further simplify the traceover the many-body fermionic Hilbert space in (A6). Tothis end we use the identitiesTr (cid:16) e ˆ B . . . e ˆ B n (cid:17) = det (cid:0) e B . . . e B n (cid:1) , (A10)Tr (cid:16) e ˆ B . . . e ˆ B n ˆ O F (cid:17) = det (cid:0) e B . . . e B n (cid:1) ×× Tr (cid:16)(cid:0) e − B n . . . e − B (cid:1) − O F (cid:17) , (A11)where the operators ˆ B i = (cid:80) x,y ˆ ψ † x [ B i ] x,y ˆ ψ y and ˆ O F = (cid:80) x,y ˆ ψ † x O x,y ˆ ψ y are the fermionic bilinear operators, andthe corresponding symbols without hats denote opera-tors on the single-particle fermionic Hilbert space withmatrix elements [ B i ] x,y and O F x,y . Correspondingly, onthe left hand side the traces are over the many-bodyfermionic Hilbert space, and the determinants and traceson the right hand side are on the single-particle fermionicHilbert space.As yet another preliminary step in the derivation ofthe CSFT algorithm, let us also decompose the gaugefields on the forward and the backward branches of theKeldysh contour into the “classical” gauge field ¯ A τx,i andthe “quantum” gauge field ˜ A τx,i as A τx,i = ¯ A τx,i + 12 ˜ A τx,i ,A N +1 − τx,i = ¯ A τx,i −
12 ˜ A τx,i , τ = 0 . . . N. (A12)Relying on the assumptions (A7) and (A8) and us-ing the identities (A10) and (A11), one can rewrite theexpression (A6) in terms of the operators on the single-particle fermionic Hilbert space and the variables ¯ A τx,i and ˜ A τx,i : (cid:104) ˆ O ( t ) (cid:105) = Z − (cid:90) d ¯ A . . . d ¯ A N (cid:90) d ˜ A ...d ˜ A N ρ EM (cid:32) ¯ A + ˜ A , ¯ A − ˜ A (cid:33) exp (cid:16) Tr ln (cid:16) u − e − h /T u + (cid:17)(cid:17) ×× exp iδ N − (cid:88) τ =0 (cid:88) x,i (cid:16) ˜ A τ +1 x,i − ˜ A τx,i (cid:17) (cid:0) ¯ A τ +1 x,i − ¯ A τx,i (cid:1) − iδ N − (cid:88) τ =0 (cid:88) x,i ˜ A τx,i J τx,i + (cid:88) j ¯ F τx,i,j − ¯ F τx − ˆ j,i,j ×× (cid:32) O (cid:32) ¯ A N + ˜ A N , ¯ A N − ˜ A N (cid:33) + Tr (cid:32)(cid:16) u − e + h /T u − − (cid:17) − O (cid:32) ¯ A N + ˜ A N , ¯ A N − ˜ A N (cid:33)(cid:33)(cid:33) , (A13)where the field strength tensor ¯ F τx,ij is constructed from the “classical” component of the gauge field ¯ A τx,i exactly8in the same way as in (10) and we have introduced theunitary single-particle forward and backward evolutionoperators u + = e − iδh (cid:104) ¯ A + ˜ A (cid:105) . . . e − iδh (cid:104) ¯ A N − + ˜ AN − (cid:105) ,u − = e + iδh (cid:104) ¯ A N − − ˜ AN − (cid:105) . . . e + iδh (cid:104) ¯ A − ˜ A (cid:105) . (A14)The path integral representation (A6) is exact in thelimit N → ∞ , δ → t = N δ (up to thesimplifying assumptions on the form of the observableoperator ˆ O and the initial density matrix ˆ ρ ), but isnot suitable for numerical analysis. The key step inthe derivation of the CSFT algorithm is to expand thefermion-induced effective action of electromagnetic field S F = Tr ln (cid:0) u − e − h /T u + (cid:1) in the first line of (A13) tothe linear order in the “quantum” electromagnetic field˜ A x,i : S F ≈ S F | ˜ A τx,i =0 + N (cid:88) τ =0 (cid:88) x,i ˜ A τx,i ∂∂ ˜ A τx,i S F (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ A τx,i =0 . (A15)In order to calculate the first derivative of S F over ˜ A τx,i ,we use the identities ∂∂ ˜ A τx,i u + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ A τx,i =0 = − iδ u (0 , τ ) j (cid:2) ¯ A τ (cid:3) u ( τ, N ) ,∂∂ ˜ A τx,i u − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ A τx,i =0 = − iδ u † ( τ, N ) j (cid:2) ¯ A τ (cid:3) u † (0 , τ ) , (A16)where we have introduced the single-particle operator ofthe conserved electric current j x,i [ A ] = ∂h [ A ] ∂A x,i (A17)as well as the single-particle evolution operator in thebackground of the “classical” electromagnetic field ¯ A τx,i : u ( τ , τ ) = e − iδh [ ¯ A τ ] . . . e − iδh [ ¯ A τ − ] , τ ≥ τ . (A18)The identities (A16) are exact up to the order O (cid:0) δ (cid:1) ,since in the derivatives of the forward and backward evo-lution operators we have used different orderings of theelementary evolution operator e − iδh [ ¯ A τ ] and the currentoperator j (cid:2) ¯ A τ (cid:3) .Using (A16) and the relations u + (0 , N ) | ˜ A =0 = u (0 , N ) and u − (0 , N ) | ˜ A =0 = u † (0 , N ), after some sim-ple algebraic manipulations we can rewrite the derivativeover ˜ A τx,i in (A15) as ∂∂ ˜ A τx,i S F (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ A τx,i =0 ≡ (cid:104) j τx,i (cid:105) == − iδ Tr (cid:18)
11 + e h /T u (0 , τ ) j x,i (cid:2) ¯ A τ (cid:3) u † (0 , τ ) (cid:19) (A19) Now that our action (A15) is assumed to be lin-ear in the “quantum” field ˜ A τx,i for τ = 1 . . . N −
1, thequantum field ˜ A τx,i can be integrated out in a straight-forward way in the case of purely bosonic observableswith O F ≡
0. The case of fermionic observableswith nontrivial O F (cid:16) ¯ A N + ˜ A N , ¯ A N − ˜ A N (cid:17) is more sub-tle, since in the path integral representation (A13) thefermionic observable itself depends on ˜ A τx,i (via the fac-tor (cid:0) u − e + h /T u − − (cid:1) − under the fermionic trace inthe last line of (A13)). It is a common assumption inthe derivation of the CSFT algorithm to neglect the˜ A τx,i dependence of the fermionic observables (see e.g.[38]), which can be justified, e.g., if the relevant physi-cal processes involve large number of virtual fermionicparticles. In this case one can argue that the expo-nent of the effective action S F has a much stronger de-pendence on ˜ A τx,i than the observable. A heuristic ar-gument in favor of such assumption is that if one ne-glects the ˜ A τx,i dependence of the fermionic observables,the expectation values of all fermionic bilinear opera-tors take exactly the same form as the expectation val-ues of the electric current (A19) and the fermionic en-ergy (see equation (15) below). Since these quantitiesare related to the observables characterizing the classi-cal electromagnetic field via the inhomogeneous Maxwellequations and the energy conservation law, they arecertainly also physical observables. While the ˜ A τ de-pendence of the observable operator might still encodesome interesting effects of the backreaction of measure-ments on the quantum evolution, taking it into ac-count would presumably lead to a significant compli-cation of the CSFT algorithm. For all these reasons,we also assume that the factor (cid:0) u − e + h /T u − − (cid:1) − in(A13) depends negligibly weakly on ˜ A τ and replace it by (cid:16) u − (0 , N ) e + h /T u (0 , N ) †− (cid:17) − .In order to integrate out the fields ˜ A x,i and ˜ A Nx,i atthe endpoints of the Keldysh contour, it is convenientto introduce the Wigner transforms ¯ ρ EM (cid:0) ¯ A x,i , ¯ E x,i (cid:1) and¯ O , (cid:0) ¯ A Nx,i , ¯ E Nx,i (cid:1) of the initial density matrix ˆ ρ EM and theoperators O F,B in (A7): ρ EM (cid:32) ¯ A + ˜ A , ¯ A − ˜ A (cid:33) == (cid:90) d ¯ E x,i ¯ ρ EM (cid:0) ¯ A x,i , ¯ E x,i (cid:1) e i (cid:80) x,i ¯ E x,i ˜ A x,i (A20) O F,B (cid:32) ¯ A N + ˜ A N , ¯ A N − ˜ A N (cid:33) == (cid:90) d ¯ E Nx,i ¯ O F,B (cid:0) ¯ A Nx,i , ¯ E Nx,i (cid:1) e − i (cid:80) x,i ¯ E Nx,i ˜ A Nx,i , (A21)where ¯ E τx,i is the “classical” electric field. We also notethat the first sum over τ in the second line of (A13) can9be rewritten as N − (cid:88) τ =0 (cid:88) x,i (cid:16) ˜ A τ +1 x,i − ˜ A τx,i (cid:17) (cid:0) ¯ A τ +1 x,i − ¯ A τx,i (cid:1) == − N − (cid:88) τ =1 (cid:88) x,i ˜ A τx,i (cid:0) ¯ A τ +1 x,i + ¯ A τ − x,i − A τx,i (cid:1) + + (cid:88) x,i (cid:16) ˜ A Nx,i (cid:0) ¯ A Nx,i − ¯ A N − x,i (cid:1) − ˜ A x,i (cid:0) ¯ A x,i − ¯ A x,i (cid:1)(cid:17) . (A22)Finally, we are ready to integrate out the “quantum”electromagnetic field ˜ A τx,i , which leads to the followingexpression for the expectation value (cid:104) ˆ O ( t ) (cid:105) : (cid:104) ˆ O ( t ) (cid:105) = (cid:90) d ¯ E d ¯ E N (cid:90) d ¯ A . . . d ¯ A N ˜ ρ EM (cid:0) ¯ A , ¯ E (cid:1) ×× δ (cid:20) ¯ E − ¯ A − ¯ A δ − δ R (cid:21) N − (cid:89) τ =1 δ (cid:20) ¯ A τ +1 + ¯ A τ − − A τ δ + δ R τ (cid:21) δ (cid:20) ¯ E N − ¯ A N − ¯ A N − δ (cid:21) ×× (cid:18) ¯ O (cid:0) ¯ A N , ¯ E N (cid:1) + Tr (cid:18)(cid:16) e h /T (cid:17) − u (0 , N ) ¯ O (cid:0) ¯ A N , ¯ E N (cid:1) u † (0 , N ) (cid:19)(cid:19) , (A23)where R τx,i = J τx,i + (cid:104) ˆ j τx,i (cid:105) + (cid:88) j (cid:16) ¯ F τx,ij − ¯ F τx − ˆ j,ij (cid:17) , (A24)and (cid:104) j τx,i (cid:105) is the expectation value of the electric cur-rent defined as in (A19). We note that the normalizationfactor Z − in (A8) and (A13) is cancelled by the zeroth-order term of the expansion (A15).From the explicit expression (A19) for the fermionicelectric current (cid:104) j τx,i (cid:105) one can immediately see that itdepends only on the classical electromagnetic field ¯ A τ (cid:48) x,i with τ (cid:48) < τ . Therefore the delta-functions in the in-tegral (A23) can be regarded as the constraints on thedeterministic evolution of the classical electromagneticfield ¯ A τx,i interacting with the quantum fermionic field.To make this more obvious, we can rewrite the chain of δ -functions in (A23) as δ (cid:2) ¯ A − A (cid:2) ¯ A , ¯ E (cid:3)(cid:3) ×× N (cid:89) τ =2 δ (cid:2) ¯ A τ − A τ (cid:2) ¯ A τ − , ¯ A τ − , (cid:104) j τ − (cid:105) (cid:3)(cid:3) ×× δ (cid:20) ¯ E N − ¯ A N − ¯ A N − δ (cid:21) , A (cid:2) ¯ A , ¯ E (cid:3) = ¯ A + δ (cid:0) ¯ E − δ R (cid:1) , A τ (cid:2) ¯ A τ − , ¯ A τ − , (cid:104) j τ − (cid:105) (cid:3) == 2 ¯ A τ − − ¯ A τ − − δ R τ − , (A25)where the last definition is for τ = 2 . . . N . From thisexpression one can see that one can sequentially inte-grate out the fields ¯ A τ with τ = 1 . . . N − A N , ¯ E N in terms of the initial values ¯ A ,¯ E . Namely, integrating out the field ¯ A first, we re-move the first delta-function in the product in (A25)and replace ¯ A by A (cid:2) ¯ A , ¯ E (cid:3) in the arguments of all the other delta functions. Integrating out ¯ A , weremove the second delta-functions and replace ¯ A by A (cid:2) ¯ A , ¯ E (cid:3) ≡ A (cid:0) A (cid:2) ¯ A , ¯ E (cid:3) , ¯ A , (cid:104) j (cid:105) (cid:1) . We can re-peat this process for all τ up to N −
1, each time ex-pressing A τ in terms of the initial values ¯ A and ¯ E andthe functionals A τ (cid:48) with τ (cid:48) < τ . It is important that insuch a sequential integration, the integrand ¯ A τ alwaysenters the argument of the delta-function being removedlinearly. Therefore despite the nonlinearity of the chainof evolution equations, such intermediate integrations donot produce any nontrivial Jacobian. To our knowledge,the absence of Jacobian in the integration measure in theCSFT algorithm so far has only been demonstrated forscalar field theory [55, 56]. It is nice to see here explicitlyits absence for lattice gauge theory coupled to fermions.After integrating out all the intermediate fields ¯ A τ with τ = 1 . . . N −
1, we are left with the following form ofequation (A23) (cid:104) ˆ O ( t ) (cid:105) = (cid:90) d ¯ A d ¯ E (cid:90) d ¯ A N d ¯ E N ¯ ρ EM (cid:0) ¯ A , ¯ E (cid:1) ×× δ (cid:2) ¯ A N − A N (cid:2) ¯ A , ¯ E (cid:3)(cid:3) δ (cid:2) ¯ E N − E N (cid:2) ¯ A , ¯ E (cid:3)(cid:3) ×× (cid:0) ¯ O (cid:0) ¯ A N , ¯ E N (cid:1) ++Tr (cid:18)
11 + e h /T u (0 , N ) ¯ O (cid:0) ¯ A N , ¯ E N (cid:1) u † (0 , N ) (cid:19)(cid:19) , (A26)where E N (cid:2) ¯ A , ¯ E (cid:3) = A N (cid:2) ¯ A , ¯ E (cid:3) − A N − (cid:2) ¯ A , ¯ E (cid:3) δ (A27)and we have expressed the functionals A N and A N − interms of the initial values ¯ A , ¯ E of the vector gaugepotential and the electric field. In this expression, itis straightforward to integrate out ¯ A N and ¯ E N , whichamounts to substituting A N (cid:2) ¯ A , ¯ E (cid:3) and E N (cid:2) ¯ A , ¯ E (cid:3) in place of ¯ A N and ¯ E N in the Wigner transforms of theobservable operators O B and O F .0To summarize, the CSFT algorithm amounts to a si-multaneous time evolution of the classical electromag-netic fields, described by the vector potential ¯ A τx,i andthe electric field ¯ E τx,i , and the quantum fermionic fields,described by single-particle evolution operator u (0 , τ ) in(A18). The discrete equations which govern this evolu-tions (arguments of the delta functions in (A23)) have awell-defined continuum limit δ → N → ∞ with fixed t = N δ , which is given by the equations (11), (12) and(13) in the main text. To simplify the notation, in themain part of the text we omit the bar over the classi-cal components of the gauge field and denote them as A x,i ( t ) ≡ ¯ A τx,i , E x,i ( t ) ≡ ¯ E τx,i , F x,ij ( t ) ≡ ¯ F τx,ij .In practice, however, the numerical solution of the con-tinuum equations (11), (12) and (13) should necessar-ily involve some discretization of time. While the sim-ple discretization of the Keldysh contour used in theabove derivation can be in principle used for numericalsolution at sufficiently small δ , for a given finite δ onecan construct different, more advanced discretizationswhich would reduce discretization errors, thus improvingthe conservation of energy (15) and making the single-particle evolution operator u (0 , τ ) numerically closer toa unitary matrix.In this work we follow [35, 38] and use the leapfrogevolution scheme for the single-particle evolution opera-tor u τ ≡ u (0 , τ ), which significantly improves the con-servation of the unitarity condition u (0 , τ ) u † (0 , τ ) = 1at finite discrete time step δ : u τ +1 = u τ − − iδh (cid:2) ¯ A τ (cid:3) u τ , τ = 1 . . . N − u = u − iδh (cid:2) ¯ A (cid:3) u , u = 1 . (A28)In practice it is convenient to work with the componentsof u τ in the basis of eigenstates of the initial single-particle Hamiltonian h (cid:2) ¯ A (cid:3) . In particular, if transla-tional invariance along some of the space directions is pre-served during the evolution, u τ remains block-diagonalin the basis of plane waves propagating along these di-rections. This block-diagonal structure can be used togreatly reduce the number of independent componentsof u τ which enter the equations (A28). We have usedtranslational invariance in two out of three spatial lat-tice directions to speed up the evolution algorithm onlarge lattices with sizes up to 200 × ×
40, assumingtranslational invariance in two out of three spatial direc-tions.For the evolution of electromagnetic field we use theequations which directly follow from (A23):¯ E τ +1 x,i − ¯ E τx,i δ == −J τx,i − (cid:104) ˆ j τx,i (cid:105) − (cid:88) j (cid:16) ¯ F τx,ij − ¯ F τx − ˆ j,ij (cid:17) , ¯ A τ +1 x,i − ¯ A τx,i δ = ¯ E τ +1 x,i , ¯ A x,i − ¯ A x,i δ = ¯ E x,i − E t/10 H EM
20 lattice with n = 10, f = 0 . µ A = 0 . v F = 1(parameter set No. 1 in Table I). − δ J x,i + (cid:104) ˆ j x,i (cid:105) + (cid:88) j (cid:16) ¯ F x,ij − ¯ F x − ˆ j,ij (cid:17) . (A29)In our simulations, we use the value δ = 0 .
05. We havechecked that decreasing δ down to 0 .
02 does not changeour results up to some small unimportant fluctuations.In principle, leapfrog-type time discretization (A28) al-lows the existence of fermionic doublers in time direc-tion - that is, the symmetric finite differences in (A28)are zero if the mode functions oscillate as ( − τ . Suchdoubler modes correspond to another flavour of Diracfermions with an opposite signature of the γ matrix.Thus if such modes are excited, they can also contributeto the anomaly equation (5) and effectively decrease theanomaly coefficient, or lead to the decay of the initialvalue of the axial charge [33]. In order to check whetherfermionic modes with such high frequencies are excitedwe have calculated the average norm of forward finitedifferences of u τ as V Tr (cid:16)(cid:0) u τ +1 − u τ (cid:1) † (cid:0) u τ +1 − u τ (cid:1)(cid:17) .Since the size of the single-particle Hilbert space is equalto 4 V this quantity should be of order of δ if u τ aresmooth functions of τ . On the other hand, doubler modeswith u τ ∼ ( − τ yield the contribution of order of unity.In our simulations we have checked that the above normremains of order of 10 − for all evolution times and doesnot exhibit any tendency to grow. This suggests that thedoubler modes remain practically unexcited during theevolution.Another important characteristic of the discretizationof the evolution equations (11) and (13) is the precisionwith which the conservation of energy (15) holds. Forthe leap-frog equations (A28) and (A29) the total energyof electromagnetic fields and fermions is conserved onlyapproximately, up to the terms of order of δ . In order toillustrate the conservation of energy in our simulations,on Fig. 14 we show the time dependence of the fermionic1energy (cid:104) ˆ H F (cid:105) , the energy H EM of electromagnetic fieldsand their total. One can see that while both (cid:104) ˆ H F (cid:105) and H EM change quite strongly during the evolution, theirsum is conserved with a very good precision, which againshows that the time step δ = 0 .
05 is small enough.The discrete evolution equations (A28) and (A29) areideally suited for parallelization on multi-node comput-ers. Indeed, the largest amount of computer time isrequired to solve the evolution equation (A28) for the(4 V ) × (4 V ) matrix u τ . Taking into account that onlyhalf of the single-particle fermionic states are filled atzero temperature, this size can be reduced by a factorof two down to (2 V ) × (4 V ). The evolution of elec-tromagnetic field (A29) requires only the total electriccurrent summed over all fermionic modes, and is com-putationally very cheap. Thus it is natural to distributethe rows of the u τ matrix over multiple nodes. On eachnode, one performs the elementary evolution step (A28)for the rows attributed to this node and calculates thepartial traces of the electric current and other fermionicbilinear operators over these rows. The results are sentto the master node, which calculates the total currentand performs the evolution of the electromagnetic field according to (A29). Since the amount of data transferredby network from each slave node to the master node is sig-nificantly smaller than the amount of data stored at eachslave mode (for realistic lattice sizes and node numbers,the number of rows of u τ per node is large), the speed ofthe algorithm scales practically linearly with the numberof slave nodes.Such parallelization also solves the problem with verylarge amount of RAM memory required to store the ma-trix u τ ( ∼
36 Gb for the 20 × ×
20 lattice when oneuses 8-byte double accuracy numbers for all fields), whichis simply split over different nodes. In order to furtherdecrease the required RAM size, we use the 4-byte floatnumbers to store u τ . We have explicitly checked that thereduction from double to float real numbers practicallydoes not affect our results.The extensive parallelization also allows us to avoid thestochastic summation over all modes [35–38], which in-troduces additional statistical noise in the results and cantherefore significantly affect potentially unstable evolu-tion which we study in this work. Instead, we perform ex-act summation over all initially occupied fermionic states. [1] M. Joyce and M. Shaposhnikov, Phys.Rev.Lett. , 1193(1997), ArXiv:astro-ph/9703005.[2] V. B. Semikoz and J. W. F. Valle, JCAP , 048 (2011),ArXiv:1104.3106.[3] A. Boyarsky, J. Froehlich, and O. Ruchayskiy,Phys.Rev.Lett. , 031301 (2012), ArXiv:1109.3350.[4] H. Tashiro, T. Vachaspati, and A. Vilenkin, Phys.Rev.D , 105033 (2012), ArXiv:1206.5549.[5] M. Giovannini, Phys.Rev.D , 063536 (2013),ArXiv:1307.2454.[6] A. Boyarsky, J. Fr¨ohlich, and O. Ruchayskiy, Phys.Rev.D , 043004 (2015), ArXiv:1504.04854.[7] G. Sigl and N. Leite, J.Cosmol.Astropart.Phys. (2015),ArXiv:1507.04983.[8] G. E. Volovik, The Universe in a Helium Droplet (Clarendon Press, 2003).[9] H. Ooguri and M. Oshikawa, Phys.Rev.Lett. , 161803(2012), ArXiv:1112.1414.[10] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosic,A. V. Fedorov, R. D. Zhong, J. A. Schneeloch, G. D.Gu, and T. Valla, Nature Phys. , 550 (2016),ArXiv:1412.6543.[11] H. J. Kim, K. S. Kim, J. F. Wang, M. Sasaki,N. Satoh, A. Ohnishi, M. Kitaura, M. Yang, and L. Li,Phys.Rev.Lett. , 246603 (2013), ArXiv:1307.6990.[12] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan,W. Wang, R. J. Cava, and N. P. Ong, Signature of thechiral anomaly in a Dirac semimetal: a current plumesteered by a magnetic field (2015), ArXiv:1503.08179.[13] Y. Hirono, D. Kharzeev, and Y. Yin, Phys.Rev.D ,125031 (2015), ArXiv:1509.07790.[14] Y. Akamatsu and N. Yamamoto, Phys.Rev.Lett. ,052002 (2013), ArXiv:1302.2125.[15] C. Manuel and J. M. Torres-Rincon, Phys.Rev.D , 074018 (2015), ArXiv:1501.07608.[16] K. Tuchin, Phys.Rev.C , 064902 (2015),ArXiv:1411.1363.[17] K. Fukushima, D. E. Kharzeev, and H. J. Warringa,Phys.Rev.D , 074033 (2008), ArXiv:0808.3382.[18] D. E. Kharzeev and H. J. Warringa, Phys.Rev.D ,034028 (2009), ArXiv:0907.5007.[19] D. Hou, H. Liu, and H. Ren, JHEP , 046 (2011),ArXiv:1103.2035.[20] D. Satow and H. Yee, Phys.Rev.D , 014027 (2014),ArXiv:1406.1150.[21] P. V. Buividovich, Nucl. Phys. A , 218 (2014),ArXiv:1312.1843.[22] I. Amado, K. Landsteiner, and F. Pena-Benitez, JHEP , 081 (2011), ArXiv:1102.4577.[23] K. Landsteiner, E. Megias, and F. Pena-Benitez, Anoma-lous transport from Kubo formulae , in Lect. Notes Phys.Strongly interacting matter in magnetic fields (Springer),edited by D. Kharzeev, K. Landsteiner, A. Schmitt, H.-U. Yee (2012), ArXiv:1207.5808.[24] P. V. Buividovich, M. Puhr, and S. N. Valgushev,Phys.Rev.B , 205122 (2015), ArXiv:1505.04582.[25] C. Manuel and J. M. Torres-Rincon, Phys.Rev.D ,096002 (2014), ArXiv:1312.1158.[26] C. Manuel and J. M. Torres-Rincon, Phys.Rev.D ,076007 (2014), ArXiv:1404.6409.[27] E. V. Gorbar, I. A. Shovkovy, S. Vilchinskii, I. Rudenok,A. Boyarsky, and O. Ruchayskiy, Phys.Rev.D , 105028(2016), ArXiv:1603.03442.[28] Z. V. Khaidukov, V. P. Kirilin, A. V. Sadofyev, and V. I.Zakharov, On magnetostatics of chiral media (2013),ArXiv:1307.0138.[29] V. P. Kirilin, A. V. Sadofyev, and V. I. Zakharov,
Anomaly and long-range forces (2013), ArXiv:1312.0895. [30] M. A. Stephanov and Y. Yin, Phys.Rev.Lett. ,162001 (2012), ArXiv:1207.0747.[31] D. T. Son and P. Surowka, Phys.Rev.Lett. , 191601(2009), ArXiv:0906.5044.[32] Y. Akamatsu, A. Rothkopf, and N. Yamamoto, JHEP , 210 (2016), ArXiv:1512.02374.[33] G. Aarts and J. Smit, Nucl. Phys. B , 355 (1999),ArXiv:hep-ph/9812413.[34] G. Aarts and J. Smit, Phys.Rev.D , 025002 (1999),ArXiv:hep-ph/9906538.[35] S. Borsanyi and M. Hindmarsh, Phys.Rev.D , 065010(2009), ArXiv:0809.4711.[36] F. Gelis and N. Tanji, Phys.Rev.D , 125035 (2013),ArXiv:1303.4633.[37] J. Berges, D. Gelfand, and D. Sexty, Phys.Rev.D ,025001 (2014), ArXiv:1308.2180.[38] V. Kasper, F. Hebenstreit, and J. Berges, Phys.Rev.D , 025016 (2014), ArXiv:1403.4849.[39] K. Fukushima, Phys.Rev.D , 054009 (2015),ArXiv:1501.01940.[40] A. M. Polyakov, Phys.Lett.B , 82 (1975).[41] M. M. Vazifeh and M. Franz, Phys.Rev.Lett. , 027201(2013), ArXiv:1303.5784.[42] P. Hosur and X. Qi, Comp.Rend.Phys. , 857 (2013),ArXiv:1309.4464.[43] A. Sekine and K. Nomura, J.Phys.Soc.Jpn. , 094710(2013), ArXiv:1309.1079.[44] N. Tanji, N. Mueller, and J. Berges, Phys.Rev.D ,074507 (2016), ArXiv:1603.03331.[45] S. A. Parameswaran, T. Grover, D. A. Abanin, D. A.Pesin, and A. Vishwanath, Phys.Rev.X , 031035 (2014),ArXiv:1306.1234.[46] N. Yamamoto, Scaling laws in chiral hydrodynamic tur-bulence (2016), ArXiv:1603.08864.[47] S. L. Adler,
Anomalies to all orders (2004), ArXiv:hep-th/0405040.[48] A. A. Anselm and A. A. Iogansen, JETP Letters , 214(1989). [49] K. Jensen, P. Kovtun, and A. Ritz, JHEP , 186(2013), ArXiv:1307.3234.[50] L. H. Karsten and J. Smith, Nucl. Phys. B , 103(1981).[51] H. J. Rothe and N. Sadooghi, Phys.Rev.D , 074502(1998), ArXiv:hep-lat/9803026.[52] M. H. Al-Hashimi and U. Wiese, Ann.Phys. , 343(2009), ArXiv:0807.0630.[53] N. Mueller, F. Hebenstreit, and J. Berges, Anomaly-induced dynamical refringence in strong-field QED (2016), ArXiv:1605.01413.[54] M. Creutz, I. Horvath, and H. Neuberger,Nucl.Phys.Proc.Suppl. , 760 (2002), ArXiv:hep-lat/0110009.[55] S. Jeon, Phys.Rev.C , 014907 (2005), ArXiv:hep-ph/0412121.[56] E. Gozzi and M. Regini, Phys.Rev.D , 067702 (2000),ArXiv:hep-th/9903136.[57] It is interesting that nonzero scalar product (cid:126)E · (cid:126)B is onlypossible for exponentially growing or decaying solutions.E.g. circularly polarized electromagnetic waves in a dis-sipationless medium always have (cid:126)E · (cid:126)B = 0.[58] This of course requires that the initial density matrixshould correspond to a sufficiently classical state, so thatthe Wigner transform of the initial density matrix is non-negative.[59] In the process of development of our algorithm we havereproduced the results of [38] on the Schwinger pair pro-duction with massive fermions and explicitly checked,that initial quantum fluctuations of electromagnetic fieldencoded in the nontrivial initial density matrix ˆ ρ EM ∼ e − ˆ H EM /T with T → (cid:126) w/ A and E0