Chiral charge dynamics in Abelian gauge theories at finite temperature

We study fermion number non-conservation (or chirality breaking) in Abelian gauge theories at finite temperature. We consider the presence of a chemical potential $\mu$ for the fermionic charge, and monitor its evolution with real-time classical lattice simulations. This method accounts for short-scale fluctuations not included in the usual effective magneto-hydrodynamics (MHD) treatment. We observe a self-similar decay of the chemical potential, accompanied by an inverse cascade process in the gauge field that leads to a production of long-range helical magnetic fields. We also study the chiral charge dynamics in the presence of an external magnetic field $B$, and extract its decay rate $\Gamma_5 \equiv -{d\log \mu\over dt}$. We provide in this way a new determination of the gauge coupling and magnetic field dependence of the chiral rate, which exhibits a best fit scaling as $\Gamma_5 \propto e^{11/2}B^2$. We confirm numerically the fluctuation-dissipation relation between $\Gamma_5$ and $\Gamma_{\rm diff}$, the Chern-Simons diffusion rate, which was obtained in a previous study. Remarkably, even though we are outside the MHD range of validity, the dynamics observed are in qualitative agreement with MHD predictions. The magnitude of the chiral/diffusion rate is however a factor $\sim 10$ times larger than expected in MHD, signaling that we are in reality exploring a different regime accounting for short scale fluctuations. This discrepancy calls for a revision of the implications of fermion number and chirality non-conservation in finite temperature Abelian gauge theories, though not definite conclusion can be made at this point until hard-thermal-loops (HTL) are included in the lattice simulations.


Introduction
Anomalous processes can be relevant in a large number of phenomena, from high energy particle physics to condensed matter. One of the most well-known applications is the explanation of the π 0 → 2γ decay in quantum electrodynamics (QED), whereas in quantum chromodynamics (QCD) they play a decisive role in the resolution of the U A (1) problem [1,2]. The rate at which anomalous processes occur is actually a relevant quantity whenever we are dealing with out-of-equilibrium processes. In highly energetic and dense matter environments like in the early universe, the fluctuations of gauge and scalar fields -sphalerons [3] -lead to rapid fermion number non-conservation in the Standard Model (SM) [4], and to chirality non-conservation in QCD [5]. The SU (2) sphaleron rate is a crucial quantity to assess the viability of electroweak baryogenesis [4], and has been extensively studied across decades, both from a purely analytical side and with the help of numerical simulations. Studies have been carried in the pure SU (2) theory in Refs. [6][7][8][9][10][11][12][13][14][15][16] and in the electroweak theory e.g. [17][18][19][20], see [21] for the latest up-to-date prediction using the measured Higgs-mass.
Anomalous U (1) processes have received some attention, especially in the cosmological context. In the electroweak theory of the SM, the anomaly in the fermionic and/or chiral current actually contains a U (1) contribution, which is associated with the hypercharge field in the Higgs unbroken phase and to the photon field of QED in the Higgs broken phase. As in Abelian gauge theories there are no large gauge transformations, nor vacuum configurations with different Chern-Simons numbers, there is no irreversible fermion (or chiral) number non-conservation, contrary to the case of non-Abelian theories. This does not prevent however the fermion/chiral number in Abelian theories to be transferred into gauge configurations carrying Chern-Simons number, and to re-appear back again due to changes in the gauge field background. These processes may have an important impact on the problems of baryogenesis [22][23][24], magnetic field generation in the early universe [25], and chiral asymmetry evolution at ∼ MeV temperatures [26]. Anomalous U (1) processes have also received a renewed interest in the quark-gluon plasma community, where the chiral magnetic effect [27,28] and its potential experimental signatures are being studied, see [29] for a review.
The above U (1) case have been studied mostly within the framework of magnetohydrodynamics (MHD), which is an effective description accounting for distance scales exceeding the mean free path of the charged particles involved in the problem, see e.g. [30,31], or [32] for recent numerical simulations. Despite the relevance of these processes, a full study beyond MHD, taking into account small scale fluctuations in detail, remains to be done. Some attempts in this direction were made in [33,34], where out-of-equilibrium techniques were implemented. The main limitations of these studies were the intrinsic numerical cost associated with a full-fledged treatment of fermions. Another approach was initiated in [35], of which this paper is a natural continuation. To explain the aim of the present research, let us set up a working model and fix notation. We are interested in the study of physics described by scalar electrodynamics coupled to a massless vector-like fermion field Ψ, so that our starting lagrangian is where F µν is the field strength tensor of the U (1) gauge field A µ , D µ = ∂ µ − ieA µ and V (φ) = m 2 |φ| 2 + λ|φ| 4 . Taking a positive squared mass m 2 > 0, we can chose its value so that the lagrangian (1.1) becomes a toy-model for the SM hypercharge sector close to the electroweak phase transition. The chiral fermionic current J µ =Ψγ µ γ 5 Ψ is not conserved at the quantum level, and satisfies the anomaly equation whereF µν = 1 2 µνρσ F ρσ is the hodge-dual of F µν , N f is the number of flavours and K µ = In the particular case of a homogeneous fermion distribution, the anomaly equation reduces to (we fix N f = 1 from now on) This allows us to write an effective description without fermions. Integrating them out, they can be represented by a (homogeneous) chemical potential µ sourced by the Chern-Simons number [36,37]. The anomaly, which represents a violation of chirality, becomes then an equation for µ, which reads The equations of motion of the scalar and gauge fields, together with the anomaly equation (1.5), can be actually derived from an effective action upon identifying ∂ 0 a = Λµ and Λ 2 = T 2 /12, with T the temperature of the system. This is nothing else but the effective action of a homogeneous shift-symmetric axion field a(t) linearly coupled to the U (1) P ontryagin densityF F , with coupling constant 1/Λ. In this paper, we want to study the evolution of the chiral chemical potential in a high-temperature environment, with the use of numerical real-time lattice simulations. We study first the evolution of µ starting from a large non-vanishing initial value. The chemical potential is expected to be unstable. Due to volume effects unavoidable in any lattice simulation, the chemical potential cannot decay completely, and rather relaxes into a lattice-dependent finite critical value µ c . By reaching large enough volumes in the lattice we have achieved nonetheless sufficiently small values of µ c , so that we can characterise well the chemical potential decay. In particular, we observe that the decay proceeds through a self-similar/power-law regime, while at the same time the gauge field long-wavelengths develop an inverse cascade spectrum.
We are also interested in the evolution of the chemical potential in the presence of a background magnetic field. In this case, the theory becomes more similar to its non-Abelian counterpart, as the vacuum becomes degenerate. In particular, µ decays exponentially all the way down to zero. We have measured the corresponding rate Γ 5 ≡ − d log µ dt as a function of the relevant parameters: the gauge coupling e and the external magnetic field B. We provide in this way a new determination of the parametric dependence of the chiral rate Γ 5 , which exhibits a best fit scaling as Γ 5 ∝ e 11/2 B 2 , with no residual volume dependence. Furthermore, we have compared Γ 5 versus the diffusion rate Γ diff of the Chern-Simons number in the absence of chemical potential, which according to a fluctuation-dissipation argument should be related as Γ 5 = 6Γ diff /T 3 , with T the temperature of the system. We confront this prediction numerically comparing the magnitude and parametric scaling of Γ 5 obtained from our lattice simulations, against direct measurements of Γ diff from (independent) simulations previously presented in Ref. [35]. The measured value of the chiral rate is in good agreement with the results of Γ diff , verifying the fluctuation-dissipation relation.
Most of what is known about anomalous U (1) dynamics comes from MHD predictions, which represent the long-wavelength effective description of our system. It is thus relevant to compare our findings with this effective theory. The phenomena we observed (e.g. self-similar decay of the chemical potential, gauge field inverse cascade dynamics, etc) are qualitatively well-modeled by MHD. This shows that they do not originate purely from the hydrodynamical regime. The expected breaking of the effective theory manifests itself when we try to measure its "matching coefficient", the electric conductivity. When confronted against our data outcome, different MHD-based observables lead to different results for the conductivity, showing that our system is away from the linear-response regime. Interestingly, the values of the effective conductivities, even if different among them, systematically appear to be smaller than the MHD value. For example, the conductivity parameter extracted from the measurement of the chiral rate is about ∼ 10 times smaller. This calls for a revision of the implications of fermion number and chirality nonconservation in finite temperature Abelian gauge theories, though admittedly no definite conclusions can be made at this point until hard-thermal-loops (HTL) are included in the lattice simulations.
On the technical side, the discretisation of this theory needs to be done with care, especially when considering the Chern-Simons number. An appropriate discretisation scheme for Abelian gauge theories, reproducing the continuum limit of the theory to quadratic order in the lattice spacing, was presented in [38]. It obeys the following properties on the lattice: i) the system is exactly gauge-invariant, and ii) shift symmetry of the axion is exact 1 . Property i) implies that physical constraints such as the Gauss law or Bianchi identities, are exactly verified on the lattice (up to machine precision). Property ii) implies that the lattice formulation naturally admits a construction of the topological number density with a total (lattice) derivative representation K ≡ F µνF µν = ∆ + µ K µ , which reproduces the continuum expression K = ∂ µ K µ ∝ E · B up to O(dx 2 µ ) corrections. Without this property, the interaction a F µνF µν cannot be interpreted as a derivative coupling, and hence is not really shift symmetric. As it is precisely the shift symmetry which justifies the functional form of the interaction in first place, it is therefore relevant to preserve exactly such symmetry at the lattice level. Hence, in the present work, we obtain results from numerical simulations based on the discretisation scheme presented in Ref. [38]. In appendix, we provide a summary of the key equations of such lattice formulation. For further technical details we refer the interested reader directly to Ref. [38]. The paper is organised as follows. In section 2, we present the outcome of our lattice simulations. Namely, we study the evolution of an initially non-vanishing large chiral chemical potential µ in the absence of an external magnetic field. This leads to a characterisation of the self-similar decay of µ, and of the gauge field inverse cascade dynamics. We then add an external magnetic field, which introduce a vacuum degeneracy, and allows for the complete decay of the chemical potential. This allows us to measure the amplitude and parametric dependence of the decay rate Γ 5 , which then we compare with Γ diff as inferred from previous simulations. Section 3 is devoted to the comparison of our results to MHD. We provide there a simplified MHD model which allows us to describe the gross feature of the system. We then characterise the inverse cascade phenomenon and see that its dynamical evolution can be well fit to an MHD-like ansatz. We also discuss the expected parametric dependence of the chiral rate. In section 4, we focus on the electric conductivity and the breakdown of MHD at the quantitative level. In section 5 we discuss our results and present future outlooks. Appendix A recaps the lattice setup and in appendix B we review the derivation of the fluctuation-dissipation relation between the chiral rate and the Chern-Simons diffusion one.

Lattice results
In this section, we report the results obtained from our lattice simulations. We select from a thermal ensemble (with µ = 0) initial configurations, which we then evolve following the (lattice version of) the classical equations of motions derived from (1.6), with initial condition µ(0) = µ 0 = 0. We follow the evolution of the chemical potential, which in the symmetric phase is expected to be unconditionally unstable, leading to the creation of long-range gauge fields [41,42]. We start by investigating the general features of the chemical potential evolution in the absence of a background magnetic field, making sure that the volume dependence of our observables is under control. Then, we study the selfsimilar behaviour of the decay for large initial chemical potential, and characterise the associated inverse cascade dynamics of the gauge field. We show that the magnetic power spectrum flows to the infrared modes. After that, we switch on a homogeneous background magnetic field. This is achieved by the use of twisted boundary conditions [43], which are already introduced in the Monte-Carlo process that generates the initial configuration (we describe this in appendix A). We observe different regimes, and in particular, we focus on the exponential decay induced by the presence of the external magnetic field. This allows us to extract the chiral magnetic decay rate Γ 5 = − d log µ dt , and characterise its parametric dependence.

Chemical potential decay
The existence of an instability can be understood by looking at the free energy of the system under consideration. In momentum space, the usual magnetic term H mag = 1 2 B 2 ∼ 1 2 k 2 A 2 competes with the chemical potential term H CS = µn CS ∼ ± e 2 8π 2 µkA 2 , where the different sign corresponds to the different polarisations of the gauge field, reflecting the chiral nature of the coupling. The chiral term may become dominant at small momenta, and for the gauge field polarisation for which H CS and H mag have opposite sign, an instability occurs for sufficiently infrared modes k < e 2 4π 2 µ = k c . 2). As expected, finite volume effects are observed for chemical potential close enough to critical. Far enough from the critical value, all lattices give the same results. Lower panel: Evolution of the chemical potential for different initial values, on the largest N = 448 lattice we have simulated. As the initial value is decreased, the size of the initial plateau increases. For large chemical potential, we observe a power law decay, which is related to some self-similar behaviour.
In particular, if all k's on the lattice are larger than the critical k, no instability can develop. This implies, as lattice momentum is discrete and has a minimum value k min , the existence of a critical chemical potential below which no instability can develop. Equation (2.1) can be rewritten to understand what is the largest chemical potential which is stable given a momentum resolution. Using the smallest momentum in a lattice k min = 2π N dx , with N the number of lattice points in one direction and dx the lattice spacing, we find For µ < µ c (e, N, dx) the instability is not captured any more for the given lattice. From now on we set T dx = 1, as any re-scaling of the lattice spacing can be translated into an inverse re-scaling of the gauge coupling constant e 2 . In this way, relevant scales are always guaranteed to be captured in the lattice (see [35] for more details). In the upper panels of figure 1 we show the behaviour of the chemical potential for different lattice sizes (in this figure and in the following, we indicate in brackets the corresponding units in powers of the temperature T ). The critical threshold predicted by (2.2) is clearly appreciated, as indicated by the dotted lines in the figure. As long as the chemical potential remains much larger than the critical value for a given lattice size N , all simulations agree, as expected, independently of the volume. On the other hand, some finite volume effects are observed around the critical value. In particular, we see in the top right panel that for the smallest volume considered N = 56, the chemical potential displays some oscillations just before approaching the critical value; such oscillations disappear for larger volumes. In the lower panel, we display the chemical potential evolution for different initial values, for the largest N = 448 lattices we simulated. For large enough initial values, we observe a power law decay, which as we will explain in section 2.2, corresponds to a self-similar behaviour. Another phenomenon we observe is that the smaller the initial chemical potential, the longer it takes for the decay to be triggered. In other words, we see the emergence of long plateaus, gradually longer the smaller the initial value of the chemical potential is. This phenomenon makes the regime of very small chemical potential difficult to be captured well on a lattice. Atop of requiring large volumes to decrease the critical chemical potential, one has to perform increasingly longer simulations.

Inverse cascade
In the previous subsection, we have seen that finite volume effects are under control for sufficiently large lattices N 50, and that long initial plateaus appear for small chemical potentials. A remarkable feature of the chemical potential evolution in this context is that for large initial values, its decay corresponds to a universal power-law. We show this phenomenon in figure 2 for different charges and initial values of the chemical potential. After some transient behaviour, all decays enter a power-law regime well described by a t −1/2 behaviour. All the dependence on the parameters comes from the charge and lies in the prefactor as we illustrate in the right-hand side part of the figure.
Actually, this phenomenon goes in pair with what is known as an inverse magnetic cascade. The decay of the chemical potential induces a transfer of magnetic energy from the ultraviolet (UV) into the infrared (IR). In other words, long-range magnetic fields are created. Let us focus our attention, for the time being, in a large volume simulation with N = 448, µ 0 = 60µ 448 c and e 2 = 1. We consider the magnetic power-spectrum where the quantity . . . |k| denotes an angular average over the spherical shell of radii [|k| − 1 2 ∆k, |k| + 1 2 ∆k), with ∆k = k min ≡ 2π N dx the binning width. This definition has the advantage of being volume independent and it is simply related to the real-space volume average as  . We see that the exponent is compatible with − 1 2 . Right: Chemical potential normalised by e 2 . The data collapse shows that µ(t) ∝ e 2 g(t) with g(t) a function independent of e 2 .
In figure 3, we plot the time evolution of the magnetic spectrum, where an energy flow from UV scales into IR scales, i.e. an inverse cascade, is clearly observed. As the chemical potential starts decaying, the magnetic energy is gradually transferred into the lower modes, so that at the end of the simulation most of it is peaked around the smallest lattice mode k min . Another way to display this information is by looking directly at the spatial distribution of the magnetic field, as we do in figure 4. There we show snapshots of the magnetic field arrows, together with some representative field lines emanating from the center. In the left panel of figure 4, we see that, immediately after thermalisation, no characteristic structures in the magnetic field distribution are observed. However, after the inverse cascade process occurs, we see long-range magnetic fields, as depicted in the right panel of figure 4, where the field lines reach out through the whole lattice. In the 3D-plot, we also show the evolution of the chemical potential (the reader is referred to figure 2 for the correct scale). As the decay proceeds, we assist to a steady transfer of energy towards the IR. In the end, most of the magnetic energy is stored in the minimal frequency. The UV part of the spectrum is related to the intrinsic UV sickness of the classical theory. and e 2 = 1. The blue 'ribbons' are field lines emanating from a small plane placed in the center of the box. They can, for instance, be drawn by following the motion of a test magnetic charge initially on the plane. At early times, no special structure is seen, since the magnetic field is homogeneous and well thermalised. At the end of the simulation, after the inverse cascade process has operated for a long while, the magnetic field has developed a long-range order, leading to structures which permeate the whole lattice. This is the spatial counterpart to the IR power displacement in Fourier space described below Eq. (3).

Decay in the presence of an external magnetic field
The dynamics of a chiral charge is of special interest in the presence of a background magnetic field B i = 1 2 ijk F jk . In this case, as detailed in [35], the vacuum of the theory is degenerate and the situation is closely related to its non-abelian counterpart. As already mentioned, such a background can be introduced on the lattice through the use of twisted boundary conditions [43]. The main difference from the case without magnetic field is that the chiral charge is now unconditionally unstable, and can decay all the way down to zero.
The early dynamics of the system depends crucially on the strength of the magnetic field. For large enough ones, it will drive the chiral charge decay from the beginning onward. For weak external magnetic fields, we expect the system to evolve initially like in the absence of magnetic field, and settle down to the critical value µ c . The chemical potential cannot stay however in the state µ = µ c in the presence of a background magnetic field, and eventually decays to zero.
These two limiting situations can be well observed in our simulations. We can also study the transition between them, see in particular figure 5. There we plot the contribution to the total energy from the chemical potential. This allows us to compare the chemical potential energy with respect to that of the external magnetic field. In the upper panel, we plot configurations with large initial values of the chemical potential, whereas the analogous plots for weaker initial chemical potential values are presented in the lower panel. As a reference, we also display the outcome from simulations with the same initial chemical potential but no magnetic field. The dashed lines represent the magnetic energy carried by the external background. Looking first at the weakest magnetic fields, we see that its influence on the chemical potential is minimal. For larger initial chemical potentials, the simulations are almost not affected while in the case of smaller chemical potentials the length of the initial plateau is slightly reduced. As expected, we see that the chemical potential relaxes first to its critical value, before undergoing a secondary decay due to the presence of the background magnetic field. For stronger background magnetic fields (top red curve, and bottom red and blue curves), we observe a transition regime where both effects are competing. In this case, the system is not sensitive any more to the finite number of infrared modes; the chemical potential does not stabilise to a critical value. In both cases, it happens when the external magnetic energy is roughly two orders of magnitude smaller than the initial chemical potential energy. For even larger magnetic energies, we enter into a different dynamical regime, where the chemical potential decays quickly to zero through damped oscillations.
It is also instructive to look at the spatial distribution of the magnetic field. In figure 6, we show the magnetic field lines corresponding to the orange and red curves in the top panel of Fig. 5, associated with large initial chemical potentials. In the left panels of Fig. 6, we show the field lines just after thermalisation, and in the right panels we show the field line distribution at the final stage of the decay. As expected, the magnetic field is initially oriented along the direction of the external magnetic field (which, without loss of generatily, was chosen to be oriented along the z-axis), and the stronger the external magnetic field the less important the thermal fluctuations are. In the final stage, we still  Figure 5: Chiral chemical potential energy evolution for a variety of external magnetic fields. Upper plot: Large initial chemical potential. We see three regimes. First, when the external magnetic field's energy is much smaller than the initial chemical potential energy, it has no effect on the initial dynamics; the system evolves as if there was no external background (orange curve). After some time spent in the critical plateau, the system eventually decays to a µ = 0 state. For a large external magnetic energy, the system does not enter the self-similar regime and the decay happens through the external magnetic field from the beginning (green curve). Between these two regimes, there is an intermediate one, where both effects contribute (pink curve). Lower plot: Small initial chemical potential. A similar discussion applies as in the upper plot, though the system exhibits now higher sensitivity to the external magnetic field, as the ratio between the chemical and magnetic energy is smaller.  see lattice-size structures. For weaker magnetic fields, the decay is essentially driven by the chiral instability of the chemical potential, so the field lines extend in all directions, as in the case without magnetic fields. Yet the effect it is still visible due to the presence of the external magnetic field, as more field lines expand along the z-direction. For a larger magnetic field background, the field lines which are not aligned to the z-axis are suppressed at the end of decay, and we essentially see only field strength lines along the z-axis.

Chiral magnetic rate
As presented in the previous section, the evolution of the system in the presence of an external magnetic field is subject to two competing dynamics. On the one hand, there is the intricate evolution related to the anomaly-induced creation of long-range gauge fields. On the other hand, the external magnetic field acts as a vacuum reservoir and induce an exponential decay of the chemical potential. In this section, we are interested in this dynamic.
To extract Γ 5 from our simulations, we proceed as follow. To disentangle the magnetic driven decay from the rest of the dynamics, we take our initial chemical potential to be the critical one. As explained in section 2.1, the anomaly-related instability does not develop for sub-critical chemical potentials. Representative evolution of the critical chemical potential in the presence of an external field are shown in the upper-left plot of 7. After an initial relaxation, we see, as expected, a clear exponential decay. This allows us to measure Γ 5 . We show its B 2 dependence for different volumes in the upper-right panel of 7. We observe a good convergence to the infinite volume limit, as the N = 448 and N = 224 cases differ only in ∼ 0.5%. We thus pursue the investigation only on the N = 448 lattices, adding a 1% systematic error to take into account any potential remaining volume dependence. We obtain that there is virtually no deviation from the expected B 2 dependence.
In the lower-left panel of Fig. 7, we look at the charge dependence of the rate. Fitting our data, we find that they are well modeled by a ∝ e 11/2 dependence. This can also be seen in the lower-right panel of Fig. 7, where we show More quantitatively, we present in Table 1 the different fits we performed to our data. All the main fits agree with each other, giving an almost exact B 2 magnetic dependence, an effective e 11/2 charge dependence and a prefactor of 10 −2.1 ≈ 0.0079. To check the robustness of the coefficients, we also to tried to fit Γ 5 e 6 with f (B) = αB β . Such a hypothesis is clearly excluded by the data as it leads to a χ 2 /dof of 155. We also present in this table a fit in Γ 5 = aB 2 e 6 ln b e to show that the charge dependence can also potentially interpreted as logarithmic corrections to a leading e 6 scaling. As we will discuss in section 3.3, this is what can be expected from theoretical predictions.
These numbers may be compared to the values of the topological charge diffusion rate, which at large time reads Upper right: Volume dependence of the measured decay. We observe convergence to in the infinite volume limit. The N = 448 results differ only in a few parts in thousands from a naive linear extrapolation using the N = 224 and N = 448 results. Lower left: Charge dependence of the decay rate. We observe a deviation from the dominant e 6 scaling. We can well describe it by an effective e 11/2 dependence. Lower right: Full scaling of the rate. We see that the quantity Γ 5 e 11/2 B 2 ext is almost constant. The remaining charge and magnetic field dependence is very weak. and it is related to Γ 5 through a standard fluctuation-dissipation argument 2 For the sake of the comparison, we present in Table 2 a re-analysis of the data of Ref. [35]. Unfortunately, the diffusion rate is a harder quantity to extract and leads 10 −3.06±0.24 · 10 −3 · B 2.05±0.08 · V 0.06±0.09 0.74 0.77 Table 2: Diffusion rate Γ diff from [35], re-analysed. We take advantage of this re-analysis to take into account the statistical fluctuations in the fit, which was not done in [35]. The errors represent two standard deviations. The relatively small χ 2 /dof shows some degree of overfitting. As we do not use it to rescale the errors, they are probably overestimated.
to a less precise quantification. As the results of [35] were obtained at different small volumes, we keep an explicit volume dependence in our fits, even if we find it to be weak. Taking as a measured value Γ 5 = 10 −2.10±0.02 B 2 e 11/2 /T 3 leads to a prediction Γ diff = 10 −2.88±0.02 B 2 e 11/2 . The prefactor and the magnetic field dependence are consistent with the measurements of Γ diff reported in [35]. At first sight, the scaling with e 2 seems to be in mild tension, but further analysis shows that the data from Ref. [35] does simply not constrain well enough the charge dependence of Γ diff . This can be seen in the fits in Table 2, where either scalings ∝ e 11/2 and ∝ e 6 are compatible with the data, exhibiting a similar χ 2 /dof . We conclude that within the errors, the fluctuation-dissipation relation Eq. (2.5) is well verified.

Comparison to MHD
Long wavelength field modes are described by hydrodynamics. A consistent way to include the electromagnetic interaction is to effectively couple it to matter through linear response theory. The equations can be further modified to take into account anomalous processes 3 [22,25]. This leads to the set of equations Despite the fact our simulations are concerned with scales outside of the MHD regime, it can be expected that some MHD features are preserved at some more fundamental scales. In this section, we will see that most of the observations we made in the previous one are also present in MHD. The fact that our approach is outside the validity range of MHD mostly manifests itself in the fact that the matching coefficient, the conductivity parameter, in this case, differs from the MHD value, see section 3. We also plot the predicted breaking time, which gives a good estimate for the end of the exponential growth. Right: Corresponding chemical potentials.

Qualitative behaviour and initial plateaus
As presented in previous sections, the chemical potential is transferred into long-range helical magnetic fields which carry a non-vanishing Chern-Simons number. To understand at a qualitative level some of the observed features, like the existence of plateaus in the initial stage, we can consider the evolution of a maximally helical field under anomalous MHD dynamics. In the presence of a background magnetic field in the z-direction, this corresponds to make the ansatz [41] which leads tof Let us first consider this toy-model without an external magnetic field (B = 0, f z (t) = 0). Then the system of equations reduces tö By plugging in the expression for µ into the first equation of (3.6), we obtain a nonlinear equation for f (t). Solving it numerically also gives us access to the evolution of µ(t). Examples of the time evolution of f (t) and µ(t) obtained with this procedure, are shown in figure 8. Despite being in qualitative agreement, this simple modeling does not capture completely the fine details, especially for large initial chemical potentials. In figure  9, we compare the modeling based on the ansatz (3.2) with the numerical outcome from our lattice simulations. The lower panel shows the chemical potential whilst the upper panel shows the Chern-Simons number, which is useful to understand the initial dynamics. Solid lines were obtained from our simulations, while dashed and dotted ones are numerical solutions to equations (3.6) for different parameters (conductivity and initial conditions). Let us first consider the case µ 0 = 8µ c , as the model gives a better description for relatively small chemical potentials. The dotted curves were obtained by fitting the ansatz model to the whole range of data. When looking at the chemical potential, we see that the model is able to describe reasonably well the initial plateau, but fails to describe well the decay. Actually, when looking at the Chern-Simons number, we see that even the initial phase is not very well described by this fit. A remedy to this is to restrict the fitting range to early times, where the system should be closer to the model. For example, the dashed curves are obtained by restricting the fitting range to times t < 10/T (delimited by a solid vertical black line). This gives a much better description of the initial phase; essentially the initial evolution is well captured by the simple ansatz. At a later stage, when a more complex dynamics has developed, such as the creation of long-range magnetic fields following a process of inverse cascade (or what can be described as a turbulent regime in MHD), the model is too simplistic to capture well the physics. The same can be said about initial chemical potentials, as in the case of µ 0 = 60µ c depicted in the figure. There, we see explicitly that the power-law decay (expected due to the generation of long-range magnetic field) is completely missed by the ansatz modeling. The red and orange vertical solid lines are a visual guide indicating the saturation of the Chern-Simons number and of the corresponding chemical potential.
We can get a better understanding of the early dynamics of f (t) and µ(t) in the ansatz modeling, by simply perturbing equation (3.6), i.e. considering f (t) = f 0 + δf (t), so that Keeping the linear order terms in δf and neglecting the terms proportional to f 2 0 (we do not start in a state with an helical field), the equation reduces to This is a driven harmonic oscillator which admits a solution as whereḟ 0 is the initial time derivative of f (t), and we have defined (3.10) Eq. (3.9) indicates, first of all, that the solution corresponds to an IR instability for the modes k < e 2 4π 2 µ 0 , which grow exponentially fast. Secondly, by estimating the range of validity of this solution, we can estimate the duration of the chemical potential plateau until the onset of its decay. Indeed, the breakdown of this approximation corresponds to the end of the exponential growth of f (t), which in turn triggers the chemical potential decay. The approximations cease to be valid when higher order perturbations become nonnegligible. To see this in detail, we keep only the exponentially growing part of the solution and neglect the constant term, The relative weight between the linear to second-order perturbation terms in the equation of motion 3.7, and linear to third order terms, are defined by the ratios We expect the linear approximation to breakdown whenever either of these ratios become of order one. We thus define the breakdown time t b to be In figure 8, we show how this prediction performs, comparing the linearisation (dotted lines) to the numerical solution to Eq. (3.6). We see that t b in (3.13) gives a reliable prediction of the range of validity of the linearisation regime. In conclusion, the duration of the initial plateau can be estimated well with our simple MHD-inspired ansatz equation (3.2).
In the presence of an external magnetic field, the system of equations (3.3) does not admit a simple enough analytic treatment. We can solve them nonetheless numerically. We show the chemical potential obtained from numerical integration of these equations in figure 10. As in the previous case, we see that it captures reasonably well the dynamics of small chemical potentials. For weak magnetic fields, it is dominated by the chiral instability of the chemical potential. The system evolves as in the absence of magnetic fields, reaches a plateau and only at later times it decays to zero because of the presence of the magnetic background field. For a sufficiently large background magnetic field, the decay is driven by the presence of such field, driving quickly the system into a µ = 0 state. For intermediate external magnetic fields, we have a superposition of the two effects.

Self-similarity and inverse cascade
The power law decay ∝ t −1/2 of the chemical potential reported in section 2.2 was also predicted in the context of MHD in several works [25,26,[45][46][47]. More specifically, various analysis of the time evolution of the magnetic spectrum, have been carried out in refs. [25,46,47]. We will now compare these modelings with the outcome from our numerical simulations. First of all, we repeat the analysis of refs. [46,47]. The starting point in MHD are equations (3.1). Taking the curl of the first equation and neglecting the second time derivative of ∂ 2 B ∂t 2 ≈ 0, the equation for the magnetic field can be recast as (3.14) We can now decompose the magnetic field in an orthonormal helicity basis Q ± ( k, x), which corresponds to a basis 4 of eigenvectors of the curl operator, In that basis Eq. (3.14) reduces to which admits as a solution Furthermore, as soon as t t 0 dt µ(t ) becomes sizable, α − (t, k) becomes sub-dominant. Neglecting the latter, we can then write Equation (3.18) represents a prediction based on MHD, which we can use to compare against our data. A further simplification we make, following [46], is to neglect the k dependence 5 of the initial amplitude α + ( k) ≈ α + 0 . This leads us to write down the following function which will be used as a fit to |B(k, t)| 2 obtained from our simulated data at different times. This is precisely what we show in figure 11. In the left panel we plot |B(k, t)| 2 for three different times together with a best fit of the form (3.19). Despite the limited number of modes available, we see that the data is well fit by a functional form as in Eq. (3.19). Let us look now at the time dependence of the fit coefficients c 1 , c 2 and c 3 . From comparing (3.19) and (3.18), we expect c 1 to be a linear function of t. This is exactly what we observe in figure 11, where the top right panel displays a linear fit from our data as 4 The precise form of the basis depends on the base space the analysis is performed, e.g. Ref. [46] considers R 3 whilst Ref. [47] used a sphere. In the analytical description above we will simply work in a Euclidean base space R 3 , as we are describing the continuum theory. Were we to write down a basis for the torus, we could construct it nonetheless from a superposition of the basis of [46]. 5 Note that we could also get rid of it by considering ratios of |B(t, k)| 2 at different times. An analysis using this method turned out to be less accurate, we believe because of the small number of momenta we can use for the fits.  . The crossing of the fit with the y-axis encodes actually various effects. First of all, we learn that t 0 is not zero, and rather corresponds to an initial time at which |B(t, k)| 2 is equal to |α + 0 | 2 (we will estimate it shortly, in a self-consistent manner, from the fit to the second coefficient c 2 ). Secondly, a non-zero crossing may also reintroduce some k 2 dependence in the initial condition, which will partially compensate for the approximation α + ( k) ≈ α + 0 . From the fit to c 1 we can extract an effective conductivity parameter, but we postpone a discussion on this for section 4.
Let us look now into the other two coefficients of the fit Eq. (3.19). The coefficient c 2 is expected to go as M (t 0 , t) ≡ 1 t e 2 4π 2 t t 0 dt µ(t ). Its dependence on t 0 allows us to estimate this initial time, by comparing it to M (t 0 , t) for different values of t 0 , and looking for the best match. This is shown in the middle right panel of figure 11. This procedure leads to an estimated value of as t 0 ≈ 9.2/T . A last check is provided by the lower right panel of figure  11, which shows how close c 3 is to a constant. In summary, our data are well-modeled by (3.19), based on the MHD solution (3.18). As we will discuss, the inferred effective conductivity parameter from this functional fit leads however to a very discrepant value with respect to standard MHD calculations. We postpone this discussion for section 4.

MHD chiral magnetic rate
Yet another phenomenon which can be described by MHD is the external magnetic fieldinduced exponential decay of the chemical potential. Neglecting the electric field time derivative from 3.1, the evolution equation for µ can be rewritten as Taking B to be constant, we then obtain The conductivity also depends on the electric charge and can be estimated 6 [48,49] σ MHD 12 4 ζ(3) 2 π 3 (3π 2 + 32) T e 2 ln(4.2/e) . (3.22) We thus see that the leading behaviour of the rate is expected to go as In section 2.4, we observed an almost perfect B 2 dependence of Γ 5 , but a certain deviation from an exact ∝ e 6 scaling. In particular, fitting to a simple power law ∝ e p , we obtain a best fit with p = 11/2, whereas fitting to a form e 6 ln q e we obtain q = 7.9, see Table 1 7 . Both of these fits describe well our data although the logarithmic fit leads to a slightly worse χ 2 /dof . Given our limited range of values of e 2 , we cannot clearly discriminate between these two options.

Conductivity
As already mentioned, MHD is specified by only one matching coefficient, the electric conductivity. A way to assess to what extent our theory can be described by MHD is to try to perform the matching. In other words, we need to find different observables which are predicted in MHD, measure them on the lattice and then extract an effective conductivity by matching the result to the prediction. Were we in the validity range of MHD, all different effective conductivities would have to be actually the same and agree with the value used in MHD computations. As we are not, we will see that it is not the case here.
In principle, the most direct way to measure an electric conductivity from our simulations is to use linear response theory and the Kubo formula [50] for conductivity. For classical field theories on a lattice with finite volume V , for a homogeneous and isotropic plasma, can be written as [35] where j i are the spatial components of the electric current, and ... is an ensemble average over different realisations of our thermal initial conditions. This formula is readily to be 6 Ref. [48] computes the leading log correction as σMHD T e 2 ln(1/e) . The numerator in the log can be refined using the full leading log numerical result evaluated at e 2 = 4π 137.04 , see table 2 in Ref. [49]. Matching it to the expression σMHD implemented as an observable in our lattice simulations. The expression is however UV sensitive (in the free theory is even linearly divergent with momentum), so the results are expected to be sensitive to the natural UV cutoff imposed by the lattice k max ∼ 1/dx. The correlator in (4.1) is expected to decay exponentially in time, Σ(t) ∼ exp (−γt), and to oscillate with the plasma frequency. When obtaining numerically Σ(t) from our simulations, it exhibits an oscillatory pattern with frequency considerably smaller that 1/T , which is presumably attributed to the plasma frequency (plots of this can be found in figure 7 (right panel) of Ref. [35], so we do not reproduce them here again). After some time the oscillations occur around zero, indicating the dumping. However, the presence of short time oscillations, associated with the lattice UV cutoff, is also very noticeable. These contribute to make the behaviour of Σ(t) very 'noisy', making difficult to obtain a trustworthy attempt to extract σ by this procedure.
Fortunately, thanks to the qualitative agreement between MHD and our results, we have other observables at hand from which we can infer an effective value of σ. In particular, we saw in section 2. Comparing these rates against our numerical fits to Γ 5 (c.f. Table 1) or against the fits to Γ diff from Ref. [35] (c.f. Table 2), we observe that the numerically extracted rates are a factor ∼ 10 larger than the MHD counterparts. In particular for e 2 = 1 and imposing an exact scaling ∝ B 2 over our data (something that it is very well verified, recall the discussion in section 2.4), we obtain 8 The two ratios in Eq. (4.5) are, as expected, consistent with each other, even if these are numbers independently obtained. This is simply due to the fact that the fluctuationdissipation argument relating Γ 5 and Γ diff , is actually well verified in our data (within the errors of the diffusion rate), recall section 2.4. If we confront our multi-parametric fits of Γ 5 [Γ diff ] against Eq. (4.3) [Eq. (4.4)] and interpret them as in the MHD expression(s) scaling with the conductivity as Γ ∝ 1/σ, this leads us to conclude that the effective conductivity of our microscopic system is a factor ∼ 10 times smaller than the MHD value, c.f. Eq. (4.2). In summary, both the chiral decay rate or the Chern-Simons diffusion rate, obtained by completely independent simulations, lead to rates much larger than predicted in the MHD picture. Then, interpreting this in light of the MHD relation between rates and conductivity, assuming the MHD suggested electric charge scaling of the conductivity, it follows that the effective conductivity obtained in this way is a factor ∼ 0.04 times smaller than in the MHD regime. An alternative way to extract an effective electric conductivity of the system is based on the fit to the coefficient c 1 in the modeling Eq. (3.19) for the chiral magnetic field inverse cascade, see section 2.2. From the linear fit in time to this coefficient we can extract, from the slope, what would be interpreted as the MHD conductivity. From the fit to c 1 we actually obtain an effective conductivity σ ef f = (0.11 ± 0.01)T , where the uncertainties are assessed by the upper and lower estimates of our fit to c 1 . This measurement is therefore very precise, but it represents however a value ∼ 0.004 smaller than the MHD conductivity [c.f. Eq. (4.2)], i.e. an order of magnitude smaller than the value inferred from the rates. This shows that our data cannot be described by a single effective coefficient and that we are indeed in a regime different that MHD.
In summary, both determinations of the effective conductivity parameter, based on measuring the decay/diffusion rate(s) and fitting the gauge field inverse cascade dynamics, are different by an order of magnitude. Both procedures indicate that the basic phenomenological picture derived from MHD is preserved in our microscopic description at the qualitative level, but confronting numerically our results against MHD-based formulae, leads systematically to infer a O(10 −1 ) − O(10 −2 ) smaller effective conductivity than expected from pure MHD considerations.
A possible explanation as to why the extracted conductivities always seems to be at least an order of magnitude smaller than the MHD predictions is that rather than using fermions we used scalar fields as our charged particles. This however should be a realistic description of the hypercharge U(1) field of the electroweak theory above the electroweak crossover, as the theory does contain relatively light scalar fields -Higgses, which are charged under this U(1). We also expect the enhancement of the chiral/diffusion rate in the electromagnetic sector right below the electroweak cross-over, due to the presence of electrically charged W-bosons, introducing similar non-linearities as in the case of a scalar field.

Conclusions
In this work, we studied the evolution of the fermionic charge in an Abelian gauge theory at finite temperature, which can be a proxy either for fermion non-conservation or for chirality breaking. In section 2.1, we studied the dependence of the evolution on the initial value of the chemical potential. We observed that for small chemical potential, it develops long initial plateaus, before the decay is triggered. We saw in section 2.2 that for larger initial chemical potential, the decay happens through a self-similar process. This leads to a phenomenon of inverse cascade in the gauge field sector, with power transferred from UV to IR scales. We observed and quantified both of these phenomena, measuring the critical exponent of the chemical potential self-similar decay to be − 1 2 . In section 2.3, we moved on to the study of the effect of an external magnetic field on the chemical potential evolution. Both for large and small initial chemical potential, three different situations can happen. First, when the external magnetic field is small, the chemical potential relaxes to its critical value through in the same way that in the absence of chemical potential. Then, from the critical value, it eventually decays to zero. For large external magnetic field, the dynamic is fully determined by the magnetic field, and the decay happens through damped oscillations. For moderate magnetic fields, both effects can be observed at the same time.
A way to isolate the effect of the external magnetic field influence is to study configurations with initial chemical potential equal to the critical one, as presented in 2.4. Then, only the exponential decay is visible. Following this procedure, we could extract a chiral dissipation rate and study its dependence on the parameters of the theory, finding Γ 5 = 10 −2.10±0.02 B 2 e 11/2 /T 3 . We also compared this rate to the fluctuation rate of the topological charge in the absence of chemical potential, which was measured in Ref. [35]. These rates are related through a fluctuation-dissipation theorem. We found a good agreement between them, providing a solid self-consistency check for our framework.
In section 3 we analysed our results in light of MHD, which is the long-wavelength effective theory of our system. Despite being out of its range of validity, qualitative features of the chemical potential can be well described. For example, in section 3.1, we provide a simplified MHD modeling which reproduces the time-scale of the observed plateaus in the chemical potential decay. In section 3.2 we also showed that the spectrum evolution of the inverse cascade is well-modeled by a MHD-inspired description, as is the chiral rate parametric dependence, which we comment upon in 3.3. In section 4, we focused on the electric conductivity and showed that our system is indeed out of the MHD range of validity: our results cannot be described by a single linear-response coefficient, as the extracted values are different among them, and systematically smaller than the MHD prediction.
This work may be continued in different directions. First, the regime of extremely small chemical potential remains to be explored. This is however challenging from a technical point of view, as the smaller the initial values of the chemical potential, the larger the volume needs to be, given the existence of a critical value µ c . Moreover, one also need longer simulations, as the initial plateau gets longer. Our computer resources thus limited us to explore systematically this regime.
An important next step will be to understand if the results presented here are robust under quantum corrections. This can be done by taking into account the effect of hard thermal loops, see e.g. [15,51,52]. Another step could be to carry on a more detailed theoretical analysis of the conductivity, as of today, the difference between the conductivity in the presence of a scalar field and of a fermion system is not known in detail.

A Lattice set-up
The results presented in this work were obtained using the lattice discretisation presented in Ref. [38]. Initial conditions for the gauge fields and the scalar fields are drawn from a thermal ensemble, generated by a simple Metropolis algorithm. Gauss law is then enforced on the thermalised configurations. After that, the system is evolved along the classical trajectory specified by the set of discretised equations of motion i,+i , (A. 6) and ∆ ± µ f = ± 1 dx (f ±µ −f ), D ± µ f = ± 1 dx (e ∓iedx µ Aµ(n± 1 2 ) f ±µ −f ) the forward/backward finite difference operator and covariant derivatives. The notation f a,µ means that the component a of the vector field f is to be evaluated at the point n +μ, withμ a unit displacement in the direction µ, f a,µ = f a (n +μ). Notice that the equations are built out of composite field so that all the fields can be expanded consistently about the same point to reproduce the continuum equations to order O(dx 2 ). The scalar fields live on lattice edges while the gauge fields are link variables; they live on half-integer sites. The time differential operators evolve by half a step in times. More details are to be found in [38].
The constant background magnetic field is introduced through twisted boundary conditions [43], which imposes a constant flux. To specify, we modify the periodic boundary conditions of the first component of our gauge field as follow This corresponds to introducing a constant background magnetic field in the z-direction of magnitude B = 2πnmag (dxN ) 2 . For large volumes, it is essential to initiate the Monte-Carlo with configurations which already satisfy the twisted-boundary conditions (A.7). This can be achieved by taking as seed A init 3 (i, j, k) = 0, (A.10) where the non-trivial dependence of A 2 ensures that A is consistent. In more details, a gauge field which lives on a periodic manifold must transform as with potentially three different gauge transformations α i , one by direction. These gauge transformations cannot be completely arbitrary, they satisfy so-called compatibility conditions. These are the requirement that the relation of the gauge field value at a given point x + N ·î + N ·ĵ to its value at x better not depend on whether one first goes through the i or through the j boundary. In other words, we have The twisted boundary conditions (A.7) tells us that ∂ 1 α 2 (1, j, k) = −

B Fluctuation-dissipation theorem
The dissipative chiral decay rate Γ 5 is related to the diffusion rate of the topological charge Q 2 (t) = ΓV t. To derive this relation, we will use Zubarev's formalism [54], in the spirit of Ref. [55]. We treat µ as a dynamical variable and want to understand its out of equilibrium properties. To do so, we introduce a "local equilibrium" partition function ρ LE = exp (−βH + V T χ(t)µ(t)) (B.1) with χ(t) some Lagrange parameter which drives µ(t) locally out of equilibrium and the factor V T 2 was introduced to make χ(t) of the same dimensions than µ. The function H is the Hamiltonian associated with our system. The relevant µ-dependent part is This function is used to compute µ at any time [54] and thus determines χ(t) as a function of µ . with · · · the usual Poisson brackets. A way to fix this is to introduce a second density matrix with N such that ρ integrates to one. Integrating by part and neglectingχ(t), i.e. considering slow processes, this can be recast as Working at linear order in h, it becomes where ρ 0 = N 0 e − H T and h 0 is the average with respect to ρ 0 .
To obtain a relation between the chiral decay rate and the diffusion rate of the topological charge, we compute μ using eq. (B.9). Using the fact that in equilibrium we have μ 0 = 0 and μµ 0 = 0, we obtain where again we assumed that χ(t) as a weak dependence on time. Now we can make use of the anomaly equation (1.5) This last quantity also appears in the diffusion rate of the topological charge. In the absence of chemical potential, Q(t) follows a random walk and we can define its diffusion rate by, following the convention of [35] Q 2 (t − t 0 ) = Γ diff V (t − t 0 ), (B.15) for times t much larger than a reference time t 0 9 . We can differentiate this expression by t 0 to get Finally, plugging this in (B.18), we get Leading to the relation