Anomalous non-conservation of fermion/chiral number in Abelian gauge theories at finite temperature

We discuss the non-conservation of fermion number (or chirality breaking, depending on the fermionic charge assignment) in Abelian gauge theories at finite temperature. We study different mechanisms of fermionic charge disappearance in the high temperature plasma, with the use of both analytical estimates and real-time classical numerical simulations. We investigate the random walk of the Chern-Simons number $N_{\rm CS} \propto \int d^4x F_{\mu\nu}{\tilde F}^{\mu\nu}$, and show that it has a diffusive behaviour in the presence of an external magnetic field $B$. This indicates that the mechanism for fermionic number non-conservation for $B \neq 0$, is due to fluctuations of the gauge fields, similarly as in the case of non-Abelian gauge theories. We determine numerically the rate of chirality non-conservation associated with this diffusion, finding it larger by a factor $\sim 60$ compared to previous theoretical estimates. We also perform numerical simulations for the system which contains a chemical potential $\mu$ representing a fermionic charge density, again both with and without an external magnetic field. When $B=0$, we observe clearly the expected instability of the system for $\mu \neq 0$, as long as the chemical potential exceeds a critical value $\mu>\mu_c(L)$, which depends on the size $L$ of the system. When $B \neq 0$, the fluctuations of bosonic fields lead to the transfer of chemical potential into Chern-Simons number for arbitrary $\mu$.


Introduction
The triangular anomaly [1,2] in fermionic current leads to many important consequences in particle physics. For example, in Abelian gauge theories, such as quantum electrodynamics (QED), it describes the decay π 0 → 2γ. In non-Abelian theories like quantum chromodynamics (QCD), it plays a decisive role in the resolution of the U A (1) problem [3,4], whereas in the electroweak theory it leads to baryon and lepton number non-conservation [5,6].
In hot and dense matter in the early universe the fluctuations of gauge and scalar fields -sphalerons [7] -lead to rapid fermion number non-conservation in the Standard Model (SM) [8], and to chirality non-conservation in QCD [9]. The existence of these transitions is associated with the non-trivial vacuum structure of non-Abelian gauge theories [10,11]. The vacuum field configurations connected by large gauge transformations have the same energy but different Chern-Simons (CS) number, allowing the disappearance of fermion/chiral number.
The situation in Abelian gauge theories is visibly different. In the electroweak theory the anomaly in the fermionic and/or chiral current contains in fact a U(1) contribution.
In the symmetric phase of the SM model it is associated with the hypercharge field, and in the Higgs phase with the electomagnetic field of QED. However, in an Abelian gauge theory, there are no large gauge transformations, nor vacuum configurations with different Chern-Simons numbers. As a result, there is no irreversible fermion (or chiral) number non-conservation associated to an Abelian gauge sector, as in the case of non-Abelian theories. This does not prevent the fermion/chiral number to be transferred into gauge configurations carrying Chern-Simons number, and to reappear back again due to the changes in the gauge field background. These processes may have an important impact on the problems of baryogenesis [12][13][14], magnetic field generation in the early Universe [15], and chiral asymmetry evolution at temperatures in the MeV range [16]; they may also be visible in heavy ion collisions [17].
In spite of the fact that a lot of work on the dynamics of Abelian gauge theories with chiral fermions has been already done, a number of questions still remain unanswered. We briefly review these questions in what follows. To set up the scene and fix notation, let us consider scalar electrodynamics with a massless vector-like fermion field Ψ, described by Lagrangian 1 where D µ = ∂ µ − ieA µ , V (φ) = m 2 |φ| 2 + λ|φ| 4 , A µ is the Abelian gauge field, F µν = ∂ µ A ν − ∂ ν A µ its field strength, and we use metric signature (−, +, +, +). The chiral fermionic current, defined as is conserved at the classical level, i.e. ∂ µ J µ 5 = 0. When quantum effects are taken into account, it has however an anomaly [1,2], whereF µν = 1 2 µνρσ F ρσ is the dual of the field strength F µν , with µνρσ the Levi-Civita antisymmetric tensor in 4 dimensions, where 0123 = +1. The chiral analogue of Lagrangian Eq. (1.1) can describe the hypercharge sector of the SM at temperatures above the electroweak cross-over; if the mass of the scalar field m is sufficiently large, it can be integrated out so that we are left with a theory including only fermions and a U(1) gauge field, hence representing ordinary quantum electrodynamics. A number of results describing the anomalous dynamics of fermions and Abelian gauge fields are available for this type of theory. These can be split between zero and non-zero temperature results: i) Zero temperature i.1) Symmetric phase. Let us take first the symmetric phase of the model m 2 > 0, when the temperature is considered to be zero. We consider an initial state containing no gauge field, but a non-zero fermionic charge density. We characterize the latter by a chemical potential µ, adding to the statistical sum the term µJ 5 0 /2, with initial value µ = µ o = 0. The left-and right-handed fermionic chemical potentials are µ L = −µ R = µ/2. A state with non-zero chemical potential is unstable against creation of gauge fields with non-zero Chern-Simons number [18]. This can be seen as follows. The (free) energy for the gauge fields with fermions integrated out, contains the Chern-Simons (CS) term [19,20] where B ≡ ∇ × A is the magnetic field. The relation between F µνF µν and the CS number are given by e 2 16π 2 F µνF µν = ∂ µ K µ , K 0 = n CS .
(1. 5) In Fourier space the term (1.4) is linear in the momentum k of the gauge field and is not positive definite, contrary to the energy density ∝ B 2 , which is quadratic in k. For sufficiently infrared (IR) modes k < k inst ≡ α π µ , (1.6) the CS term H CS dominates the free energy, leading to an instability of the gauge field, which grows exponentially as [18] A k ∼ e ω k t , ω k = k(k inst − k) . (1.7) Eventually, the fermionic number disappears as it is converted into long-ranged configurations of the gauge field, with a density of Chern-Simons number equal to the initial chiral density n F , 8) and the typical scale of the order of the system size L, with typical magnetic field B ∼ 4π e n F L and gauge field A ∼ 4π e √ n F L [18]. While the initial energy density of the system was ρ F ∼ µ 4 o , the energy density of the gauge field carrying the same fermionic number is now smaller, as it is suppressed by the size of the system, ρ A ∼ µ 3 /Le 2 .
i.2) Higgs phase. In the Higgs phase, when m 2 < 0, an analogous instability can only be developed if [18] with M A = ev the mass of the vector field in the broken phase, where the scalar Higgs-like field takes a vacuum expectation value v 2 = − m 2 λ . The growth of the instability occurs in this occasion as for IR modes obeying the relation (1.11) For small M A k inst , we have M 2 A /k inst < k < k inst − M 2 A /k inst .
ii) Finite temperature ii.1) Symmetric phase. Let us take again the same initial state in the symmetric phase, but now at finite temperature. A linear analysis similar to that carried out at zero temperature, reveals the instability of the system for long-range gauge field modes with momenta k < k inst . Since Abelian magnetic fields in a high temperature plasma are not screened (this can be proven in all orders of perturbation theory [21,22], with non-perturbative lattice simulations confirming the same behavior [23]), the instability takes place for any IR mode k < k inst . The rate of the instability growth is, however, different from that at zero temperature. One can distinguish three qualitatively different regimes: a) In the first one, the typical instability length scale l inst ∼ 1/k inst is larger than the mean free path λ ∼ 1/(α 2 T ) of the fermions in the plasma. This happens if the chemical potential is sufficiently small, µ < e 2 T . For the analysis of instabilities one can use the effective description of long range modes in a plasma, namely magneto-hydrodynamics (MHD). The rate of the instability growth is suppressed in comparison with the zero temperature case by the electric conductivity σ of the plasma [12,15], like where k inst is still given by Eq. (1.6). For QED with N l fermionic flavours [24] σ 3ζ(3) log 2 , (1.13) where ζ is the Riemann ζ-function, so 3ζ(3)/ log 2 5.2. The result (1.12) can be derived easily from Eqns (2.10), which we introduce later on when discussing the magnetohydrodynamical regime of the system.
b) The second regime corresponds to the instability length shorter that the fermions mean-free path but larger than the Debye screening scale λ D ∼ 1/(eT ), λ D < l inst < λ. This corresponds to chemical potentials e 2 T < µ < T /e. The instability growth rate can be found using the thermal photon propagator, and is equal to [25] A k ∼ e ω k t , ω k = 4k 2 πm 2 where m D is the Debye mass. For the theory with Lagrangian (1.1) and in the one loop approximation, one finds m 2 D = 2e 2 T 2 /3. c) For even larger chemical potentials µ > T /e the instability length is shorter than the Debye radius, and the zero-temperature analysis is applicable. The development of an instability of IR modes of the gauge field is described, in this case, again by Eq. (1.7).
The transition between different regimes a), b) and c) is smooth, and the expressions for ω k are parametrically the same at the matching points µ ∼ e 2 T and µ ∼ T /e. In all cases the density of real fermions eventually disappears and is replaced by a Chern-Simons condensate of the gauge field.
ii.2) Higgs phase. To the best of our knowledge, a detailed analysis of instabilities, in what concerns the Higgs phase, has not been carried out. We expect that the condition for the instability to develop given in Eq. (1.9) still remains in force, but the rates Eqs. (1.12,1.14) are to be modified.
As we have already mentioned, these behaviors are different from what happens in non-Abelian theories in similar situations. For instance, non-Abelian gauge theories in the symmetric phase are confining at zero temperatures. At non-zero temperatures, still in the symmetric phase, the non-Abelian magnetic fields acquire a "magnetic" mass m mag ∝ g 2 T [26], where g is the non-Abelian gauge coupling 2 . We expect this to lead to a threshold similar to Eq. (1.9) with replacement of M A by the strong coupling scale Λ similar to Λ QCD at zero T , or by m mag if T >> Λ. On top of that, in the non-Abelian case the final state of the system after the development of the instability does not contain long-range non-Abelian fields, as they are screened by m mag . The system contains instead vacuum configurations with non-zero topological charge.
Moreover, at non-zero temperature, in addition to classical instabilities, there are thermal fluctuations -sphalerons -which also drive the system to a state with zero net fermionic charge. The rate of these fluctuations Γ sph does not depend on the chemical potential of fermions (at least, for small enough µ), and the chemical potential behaves as µ ∝ exp(−κΓ sph t), with κ ∼ 1 a known number which depends on the matter sector of the corresponding gauge theory [28].
A lot of work has been done for the study of fermion number non-conservation and sphaleron transitions in non-Abelian theories in the past. We mention just a few. The prefactor for the sphaleron rate in the Higgs phase of the SM was found in [29,30]. In [31] it was demonstrated that the sphaleron rate in the symmetric phase scales with the SU(2) gauge coupling like α 5 W , contrary to α 4 W expected from naive scaling [28]. In [32] Bödeker argued that the rate has an extra log(1/α W ) enhancement and suggested an effective field theory description accounting for this effect, which has been worked out further in [33,34]. The numerical simulations of sphaleron transitions in 1+1 dimensions were initiated in [35] and carried out in [36]. In 3+1 dimensions the first lattice simulations were done in [37], with many improvements accounting for hard thermal loops and Bödeker effective theory appearing later in [38][39][40]. The ultimate results for the sphaleron rate were reported in [41]. The combined dynamics of the fermionic chemical potential was addressed in lattice simulations in [37] and refined considerably in [42].
Much less efforts were invested to the study of the Abelian case. Besides the works we have already mentioned, Ref. [43] underlined the potential importance of fluctuations of the electromagnetic field for the problem of the chiral charge erasure. This question has been also studied recently in [44]; our discussion of the same phenomena in the present paper is considerably different.
To summarise, the aim of the present work is to elucidate the difference between the Abelian and non-Abelian theories in what concerns the behaviour of the fermionic number at high temperatures. This happens to be not as trivial at it may seem. As we will discuss in Section 2, the ground state of an Abelian theory at non-zero temperatures may have more in common with non-Abelian theories than normally considered, and thus potentially could lead to other possible mechanisms for anomalous fermion number nonconservation, in supplement to the instabilities discussed above. This certainly happens in the presence of non-zero magnetic fields, leading to the ground state degenerate with respect to Chern-Simons number. To account for fluctuations of Abelian gauge theories at non-zero temperatures, we propose to study numerically a classical Abelian Higgs model with a Chern-Simons term, replacing the theory described by Eq. (1.1) with fermions, by a theory in which fermions with non-zero density are integrated out. This theory and its lattice implementation are discussed in Section 3. In Section 4 we study the diffusion of the Chern-Simons number in the absence of the chemical potential, both without and with an external magnetic field. The rate of CS diffusion is known to give the rate of fermion number dilution in non-Abelian gauge theories [28], and a similar relation is expected to hold in the Abelian case as well. We find the rate of CS diffusion as proportional to the square of the magnetic field (an expected result [12,15] but quantified with novel inputs for the pre-factor). In Section 5 we carry out an analysis of the dynamics of the chemical potential when we assume its initial value is µ o = 0, again both in the absence and presence of an external magnetic field. In Section 6 we summarize our results, and discuss their implications and limitations.
2 Dynamics of the chiral charge and magnetic fields: analytical estimates

Abelian configurations with CS number
The Abelian Chern-Simons number density (1.4) vanishes for zero energy configurations. This means that the bosonic ground state of the system at zero temperature is unique (in the absence of an external magnetic field). This is opposed to the non-Abelian case, where pure gauge configurations with non-trivial topology carry an integer CS number. The situation changes if we have a non-zero magnetic field B. A non-zero uniform magnetic field can in fact play the role of the ground state of an Abelian theory, as this configuration satisfies the equations of motion and is stable due to magnetic flux conservation.
For definiteness, let us consider an external magnetic field in the direction of the third axis z, B = (0, 0, B 3 ). Now, if B 3 = 0, the ground state of the system acquires a degeneracy with respect to the CS number: a non-zero constant gauge field A 3 does not cost any energy but leads to a configuration with non-zero n CS ∝ A · B. Contrary to the case of a non-Abelian theory, N CS = d 3 xn CS is not quantized and can have any value. This type of configuration can serve as an infinite reservoir to absorb the fermionic charge, exactly as for the case of non-Abelian theories. One should expect, therefore, that the dynamics of the fermion charge and of the Chern-Simons number, is qualitatively the same for Abelian and non-Abelian theories at high temperatures, in the presence of a magnetic field in the symmetric phase.
The case when an external magnetic field is absent is more subtle. Consider the gauge field with the typical amplitude A and variation scale l. The energy density of this field is ρ A ∼ A 2 /l 2 , and it can carry a CS number density The dependence of the energy density on the CS number is therefore infrared-sensitive as In other words, a configuration with a fixed CS number density n CS = 0 may have arbitrarily small energy density for a sufficiently long-wave vector field. Since at non-zero temperatures the energy density of the system is different from zero, the thermal ground state may contain configurations with arbitrarily large CS numbers in the limit of infinite volume l → ∞. This looks pretty much similar to the non-Abelian case, where the same is true, though not only the energy density, but also the total energy of a configuration carrying a Chern-Simons number may be vanishing. One may wonder therefore whether in the absence of background magnetic field there might be yet another mechanism for fermion number transfer to Abelian fields, not related to the instabilities discussed in the Introduction. In particular, do the fluctuations of the Abelian field, similarly as the sphalerons of non-Abelian fields, play any role? We discuss this question next.

Instabilities or fluctuations?
The rate ω k at which the fermion number is 'eaten up' by the gauge field due to the instabilities described in Section 1, is proportional to µ n , with n changing from 1 to 3, depending on the situation. As it is well known, and it has been already mentioned in the Introduction, in non-Abelian theories there is yet another mechanism for fermion number dilution, related to thermal fluctuations -sphalerons -with a µ-independent rate, at least for small µ. This can be seen at a qualitative level, from the following considerations [31,33]. Let us take the fluctuation of a gauge field with a typical size l and amplitude A. It has a CS number N CS ∼ n CS l 3 and energy E ∼ ρ A l 3 , see Eqs. (2.1), (2.2). Fluctuation that can change the fermion number should have N CS ∼ 1, leading to a relation A ∼ π le . Their energy is, therefore, of the order E ∼ π αl . The probability to have such a fluctuation at temperature T is given by the Boltzmann exponent exp(− E T ), so the typical size of the required fluctuation should be at least of the order of l ∼ 1/(αT ), with a number density n sph ∼ (αT ) 3 . To get the rate of the diffusion of the topological number, what is left to do is to divide n sph by the typical time t sph ∼ 1/(α 2 T ) of the ∆N CS = 1 changing transition (one power of α comes from the dimensional analysis of effective 3-dimensional theory [28], while the extra power accounts for slowing down of the process by the medium [29,31]). As a result, the rate (per unit volume) of fluctuations changing the CS number by one unit, is parametrically given by Γ sph ∝ α 5 T 4 , up to a logarithmic factor that comes from a more refined treatment [32]. In a plasma carrying a non-zero fermionic charge, the rates leading to its increment or decrement differ by µ/T , as follows from Eq. (1.4). This leads to the conclusion that µ(t) ∝ exp(−Γ sph t/T 3 ), with a µ-independent exponential.
In the previous consideration the non-Abelian character of the gauge theory was not used at all. It seems that it may be equally applied to the Abelian theory as well. Qualitatively, the difference appears in the subsequent development of the fluctuation. In the non-Abelian case, after crossing the "sphaleron barrier" it may evolve to a new vacuum state with the different CS number but zero energy, thermalising its energy. In the Abelian case, such a vacuum state does not exist. The discussion in Section 2.1 indicates however that the effective degeneracy with respect to the CS number, appears in the limit of long wavelengths of the gauge field. It is an open question whether this effect may lead to dissipation of the fermionic number in the Abelian theory, with a rate that does not depend on µ. If the non-Abelian estimate were applicable to the Abelian case as well, the Abelian "sphaleron" rate Γ sph ∝ α 5 T 4 would supersede the rate associated with instability (1.12) for a chemical potential µ < αT . Whether this is true or not, this question can be studied, in principle, in real time lattice simulations. Unfortunately, due to limited computer resources, we could not get an answer to this, because of the reasons we explain in Section 3.1. We thus leave open the investigation of this matter for the future. We have got, however, indirect evidence (Section 3.1) that the Abelian rate may scale like Γ sph ∝ α 6 T 4 , i.e. with one extra power of α compared to the non-Abelian case. If true, this will exceed the rate associated to the instability for small µ < eαT .

Diffusion of CS number and chirality non-conservation
For physics applications, the main quantity of interest is the time evolution of the fermion number and of the magnetic fields. This is a complicated problem involving many different time and length scales operating at different stages of the equilibration process. However, in the limit of small chemical potential µ, the rate of fermion number non-conservation can be found within the pure bosonic theory, by considering the diffusion of the CS number and the use of the fluctuation-dissipation theorem [28,42].
To make a proper comparison between Abelian and non-Abelian cases, consider bosonic theories, without fermions, based on either a SU(2) gauge theory with Higgs doublet, or on a U(1) theory. The main quantity which determines the dynamics of the topological transitions is the CS diffusion rate Γ [28], which characterizes the expectation value is a gauge-invariant quantity related to the Chern-Simons number, with the index a omitted in the U(1) case, and running like a = 1, 2, 3 in the SU(2) case. Homogeneity in time and space leads to Q cont (t) 2 = ΓV t , (2.5) where for t → ∞, As we briefly reviewed in Section 2.2, it is well known that in the symmetric phase of a non-Abelian theory, Γ depends on the coupling constants as Γ ∝ α 5 log(1/e)T 4 . (2.7) In the Higgs phase it is suppressed by the Boltzmann exponent exp (−M sph /T ) [8], where M sph ∝ M A /e 2 is the sphaleron mass, and M A is the temperature dependent mass of the vector boson.
In Abelian theories, and in the absence of a magnetic field, we expect Q cont (t) 2 to become constant at large times t, because contrary to the non-Abelian case, to have a nonzero CS number now costs energy. This constant is estimated in Section 3.1. We expect however a diffusive behavior if an external magnetic field is present, since the CS number can grow without limit with the energy fixed.
We show now that the perturbative contributions to the diffusion rate at B = 0 vanish. For this end let us consider the expansion of the correlator Eq. (2.6) with respect to an external magnetic field in theẑ-direction, writing B z =B + δB. The lowest order term in B is zero due to the Abelian character of the theory, while the first order term vanishes as well, due to symmetry considerations (isotropy). We are thus left with the second order contribution inB, In perturbation theory this correlator is zero. This is most easily seen in the gauge A o = 0, where the electric field can be written as E 3 (x, t) = ∂ o A 3 . The correlator in Eq. (2.8) becomes where the first term is zero due to the cluster property, while the second can be written as 1 2 ∂ 0 A 3 (0) 2 , and hence it also vanishes.
To get a non-perturbative contribution, we start from an estimate of the diffusion rate which can be obtained from previous results [12,15] based on the equations of magnetohydrodynamics. In the presence of a homogeneous chemical potential for the chiral charge, the effective action for the electromagnetic fields, leads to a modification of Maxwell equations for wavelengths larger than the fermions mean free path λ ∼ 1/(α 2 T ) in the plasma. Defining the electric and magnetic fields as the modified equations read [12,15] [we use a metric signature (-,+,+,+)] where σ is the electric conductivity of the plasma, and we have assumed that the density of electric charge is zero, and the plasma has zero velocity. This system of equations is complemented by an anomaly equation like Eq. (1.3), which rewritten in terms of the chiral chemical potential reads Neglecting the time derivative of the electric field in Eq. (2.10), one gets an equation for µ as This shows that, in the presence of an external homogeneous magnetic field, the contribution of long-range fluctuations of gauge fields to the rate of chirality non-conservation, is (2.13) The second term in (2.12) can be considered as a source leading to non-zero chemical potential in the presence of helical magnetic fields, with interesting dynamics of exchange between the chiral charge and magnetic fields [12][13][14][15][16].
The fluctuation-dissipation theorem allows to relate the diffusion rate Γ (per unit time and volume) of the CS number with the dilution rate Γ 5 (per unit time) of the chemical potential. Standard considerations following [28,42], give for the theory (1.1) the following relation (2.14) Using (2.14) and Eq. (2.13), we finally arrive at The rate scales as α 3 and is proportional to B 2 . We present below yet another estimate inspired by [31], making use of the fluctuations discussed in Section 2.2 with typical size parametrically smaller than the mean free path, l sph ∼ 1/(αT ) λ. For this end we take as in Section 2.2 A ∼ π e l sph and ∂ o ∼ 1/t sph ∼ α 2 T , leading to the rate (2. 16) This result has the same parametric dependence on the magnetic field as eq. (2.15). However, it contains the square of the fine structure constant contrary to α 3 in (2.15). In principle, there is no contradiction, as (2.15) only accounts for sufficiently long range fluctuations (for which the Eqs. (2.10) are valid) and when ∂ E ∂t = 0. Still, this difference is perturbing and calls for reconsideration of the typical time scale in Abelian gauge theories, which may be different from the non-Abelian case because of the different vacuum structure in the absence of a magnetic field. Assuming that Eqs. (2.10) with µ = 0 and e j = −σ E were valid (at least parametrically) for short ranged fields as well, we get that the typical electric field E associated with magnetic field of strength B ∝ A/l ∝ e 3 T 2 and size l sph , is of the order E ∼ B/(σl) ∝ A/(σl 2 ) ∝ ωA with ω ∝ α 3 T , i.e. α times slower than the non-Abelian one. If correct, the rate from fluctuations is proportional to an extra power in α, and has the same parametric dependence as (2.15). If we replace the external magnetic field by a typical magnetic field B ∼ A/l ∝ e 3 T 2 of the fluctuation which carries N CS ∼ 1, we will get an "Abelian sphaleron rate" as which is suppressed in comparison with the non-Abelian one by an extra power of α. In the estimates above we assumed that the electric conductivity scale like σ ∝ 1/α even at small distances. The real time lattice simulations of the CS diffusion in the presence of magnetic field are considered in Sections 4 and 5. We will see there that the parametric dependence on B and α coincides with (2.13, 2.17), but with a numerical pre-factor exceeding substantially that in (2.15), indicating that the small scale fluctuations are more important.

A model for real-time classical simulations
It is notoriously difficult to make lattice simulations with dynamical massless fermions. We hence use an effective field theory approach, where fermions are integrated out, but their presence is accounted for by a homogeneous chemical potential µ for their chiral charge. As mentioned before, in this approximation the energy of the system acquires the term Eq. (1.4), and the theory becomes purely bosonic, so that it becomes suitable for classical real time simulations. We can then put the resulting bosonic theory in a finite volume V , with periodic boundary conditions for gauge-independent quantities in all spatial directions.
An economic way to derive, in an explicit gauge-invariant way, the equations of motion for the different fields involved in the bosonic approach, is to use an (auxiliary) axion field a(x) with action where c 2 s and M are respectively dimensionless and dimensionfull parameters, to be fixed later on. Due to the fact that F µνF µν is a 4-divergence, this action is invariant under the shift symmetry a → a + const. The equations of motion following from action (3.1) are where the (unit-charge) current is defined as j µ = 2Im{ϕ * D µ ϕ}. As expected, these equations only contain derivatives of the field a(x). Now, the anomaly equation Eq. (1.3) for the chiral current J 5 µ can be compared with Eq. (3.4), leading to an identification We consider only the homogeneous 3 case ∇a = 0, so that J 5 = 0. As a result, the equation of motion for the axion field Eq. (3.4), turns into At the same time, the equation of motion of the gauge field, in the presence of a homogeneous axion (chemical potential), changes to (3.7) The relation between the chemical potential and the chiral charge at T = 0, for small µ, reads where we have used F µνF µν = 4 E B andF i0 = −B i . Besides, the system is characterized by an integral of motion -the total energy -given by This quantity is important for monitoring the numerics as any time discretization induces always some degree of violation of energy conservation. This system of equations ensures also the magnetic flux conservation, for example ∂ ∂t dxdyB 3 = 0 . (3.14) Thus, the initial conditions for the time evolution can be characterised by giving the energy of the system and of the magnetic flux through one of the planes, which we can always choose to be that of xy. The equations of motion Eqs. (3.9)-(3.12) and the Hamiltonian Eq. (3.13), serve as the theoretical basis for our numerical simulations. The way how these equations can be put on the lattice while keeping the gauge invariance, the topological character of FF , and a non-zero magnetic flux, is discussed in detail in [45]. In Sect. 3.1 we discuss some of the limitations naturally arising in numerical simulations in a lattice with a finite volume. In Sect. 3.2 we summarize the lattice formulation of the theory based on [45], whereas in Sect. 3.3 we discuss the choice of parameters to capture appropriately the correct physical regimes in the numerical simulations. In Sect. 3.4 we discuss how we set up the initial conditions of the system in the lattice.
In Section 4 we present our numerical results, studying first the diffusion of the Chern-Simons number in the absence of fermions (i.e. with µ = 0). With standard Monte-Carlo techniques we create a set of configurations with probability exp(−E/T ), which are then evolved in time with the use of the lattice equations of motion. We consider both cases when an external magnetic field is absent (Sect. 4.1) and present (Sect. 4.2). We determine the correlator Eq. (2.5) from which we extract the diffusion rate. In Section 5 we study the dynamics of the chiral charge µ when an initial value µ o = 0 is provided. To follow the non-equilibrium evolution of µ we first prepare an equilibrium configuration with (a lattice version of the) Hamiltonian given by Eq. (3.13) with µ fixed to zero. Following, we introduce µ = µ o = 0 as an initial condition, and evolve the fields and µ with the use of the discretized version of Eqs. (3.9)-(3.12). We analyze the evolution of the chemical potential, and of the scalar and gauge fields, considering again both the absence (Sect. 5.1) and presence (Sect. 5.2) of an external magnetic field.

Finite volume effects
As the numerical lattice simulations can only be done in a finite volume, there are peculiarities of the system which are modified with respect to the continuum. A main feature of this kind is the following fact: in the absence of an external magnetic field, the transfer of the chiral charge (chemical potential µ) into the gauge field by any type of mechanisms discussed (instabilities or fluctuations), has a volume-dependent threshold effect. Namely, considering a simulation box of length L and volume V = L 3 , there will be no dilution of the chemical potential for values below a critical amplitude µ < µ c (L). To find µ c (L), let us consider a system with initial µ o = 0 and zero CS number density. The change of the energy of the system due to the small change of the chemical potential δµ is see Eq. (3.13), and the change of the CS number of the system is see Eq. (3.12). Now, the transfer of the fermionic charge into bosonic fields is energetically favourable only if δE is greater than the change of the magnetic energy associated with CS number (3.16).
With the use of different functional inequalities one can derive the following bound for the magnetic energy of configuration carrying non-zero CS number (see [46] and references therein), 18) and the minimal value of the numerical constant C that can be found in the literature [47] is C = 1/(256π 4 ). This inequality is true for an infinite space where the gauge potentials are continuously differentiable functions with compact support; it is also valid for strictly periodic boundary conditions for gauge potentials in a finite box. With the use of inequalitȳ Eq. (3.17) can be put in the form The requirement δE > E M for CS number given by (3.16) leads to the determination of the critical chemical potential, giving if we take the value of C quoted above. Interestingly, pretty much the same value of µ c comes from the trivial analysis of instabilities in the free system. The main difference with the infinite volume case is that the spatial momentum is now discrete, where n is an integer number. The instability discussed in Section 1 can only develop if k min = 2π/L is smaller than k inst , leading to conclusion that only the states with sufficiently large chemical potentials are unstable, i.e.
This is larger than (3.22) just by 14%. We suspect that the value of C for the specific geometry of the finite cubic box can be lowered down to C → C 2/π, as our numerical simulations described in Section 5 are consistent well with the bound (3.24). If the initial chemical potential exceeds µ c , then in the course of the system evolution µ(t) will be bounded from below by the same value µ c (in reality, it will typically oscillate around µ c , until it finally relaxes to it). Of course, in the infinite size limit L → ∞ the threshold value vanishes.
The inequalities (2.15,3.24) make very demanding the non-perturbative lattice analysis of the system with physically interesting small chemical potentials µ T . Indeed, the lattice size must verify Now, let us discuss whether we can disentangle the change of chemical potential due to instabilities from that due to fluctuations. Parametrically, the rate of change of the chemical potential from fluctuations is expected to be given by Eq. (2.18), whereas the rate due to instabilities for small enough chemical potentials is coming from Eq. (1.12) is Γ inst ∼ k 2 inst /(4σ). The rate associated with fluctuations dominates for µ T α Unfortunately, due to lack of computing power, we could not reach this condition in our numerical simulations. Note that the analysis presented above does not apply if our finite system contains a constant background magnetic field. In this case the inequality (3.17) is not valid, and the CS number of the system is not bounded from above for a fixed value of the system energy. We will explicitly show this in Sect. 5.2, where µ will decay completely down to zero due to the presence of an external magnetic field.
We finish this subsection with the estimates of the fluctuations of the CS number in an Abelian free field theory without fermions, to test the correctness of the lattice simulations in a limiting case. In the absence of background magnetic field the computation has been done in [37] and reads The fluctuations do not grow in time contrary to non-Abelian theories. If the magnetic field B is present one easily finds where the growth as time square comes from the homogeneous electric fields, always present as fluctuations in the system.

Lattice formulation
In the lattice we do not consider summation over repeated indexes, as this can lead to confusion. We work in a cubic lattice of length L = N dx, with dx the lattice spacing and N the number of points per dimension. We consider periodic boundary conditions. A lattice point n = (n o , n) = (n o , n 1 , n 2 , n 3 ) displaced in the µ−direction by one unit lattice spacing, n +μ, will be referred simply as n + µ or by +µ. For example, ϕ +µ ≡ ϕ(n +μ).
Components of gauge fields live in between lattice sites in the direction of the component, We define lattice ordinary and covariant derivatives, forward (+) and backward (-), as where a link is defined as ) . An Abelian U (1) gauge transformation in the lattice corresponds to with β(n) an arbitrary function of the space-time site, so that the links and covariant derivatives transform as Using these transformation rules we can build a gauge invariant lattice action mimicking the continuum action Eq. (3.1). We first build gauge invariant electric and magnetic fields as For convenience, we also define the following field combinations Considering the auxiliary axion as a homogeneous field so that its time derivative represents the chemical potential of the system, We derived in Ref. [45] a lattice action built out of two pieces, an Abelian-Higgs part and a Chemical Potential part, given by The lattice action Eq. (3.40) exhibits exact gauge and shift symmetries in the lattice, and reproduces the continuum action Eq. (3.1) to order O(dx 2 µ ). Besides it gives rise to a set of lattice equations of motion compatible with i) the Bianchi identities 4 , and ii) an explicit iterative scheme to solve them. We note that the choice of the lattice operator , one obtains the set of lattice equations of motion mimicking a continuum system with chemical potential at finite temperature T . In the Coulomb gauge A o = 0 (so that U o = 1), these equations read [45] Equation Natural Site We have indicated in the rhs the common natural space-time site to all terms belonging to a given equation of motion, around which we can expand each equation and reproduce the continuum analogue Eqs. (3.9)-(3.12) to order O(dx 2 µ ). A lattice definition of the Chern-Simons number follows naturally from Eq. (3.42) as In Ref. [45] we showed how other lattice representations of FF often fail to fulfill the Bianchi identities, introducing extra terms in the lattice equations of motion which do not represent well the continuum equations. 5 For this it is crucial that we use a non-compact formulation for the gauge sector.
is an initial constant. This lattice representation of the Chern-Simons number reproduces to order O(dx 2 µ ) the continuum expression More importantly, our lattice definition admits a total (lattice) derivative representation as which is a local identity at every lattice site, with K µ defined by components like Given that we consider periodic boundary conditions, n i ∆ + i K i = 0 (no matter the expression for K i ). We then arrive at a expression for the Chern-Simons number (say after p time iterations) as depending only (as it should) on the difference between the final and initial values of K 0 , mimicking the continuum result 4π 2 N cs = 1 Let us also note that Eq. (3.45) represents an exact solution to the lattice equation of motion for the chemical potential (say after p time steps) as which mimics the continuum relation given by µ(t) = µ(0) + 12 Therefore, in our lattice simulations, whenever we want to obtain the Chern-Simons number, we can either call Eq. Finally, let us mention that if we require an external magnetic field in the system, we can introduce it in the lattice, following [48], by using twisted boundary conditions for the (relevant) component of the gauge field. For instance, without loss of generality, we can consider an external magnetic field in theẑ direction, B ext = (0, 0, B ext ). In this case we only need twisted boundary conditions for A 1 : whenever we want to calculate thê z-component of a magnetic field B 3 at a location (n 1 , n 2 , n 3 ) = (1, N − 1, n 3 ), instead of simply computing (∆ + 1 A 2 − ∆ + 2 A 1 ), we should really make the calculation where n mag is a positive integer number. Condition Eq. (3.50) is, in fact, equivalent to assume that the first component of the gauge field makes a 'discrete jump' dx , which represents precisely the required twisted boundary condition for the gauge field, in order to introduce an external magnetic field in theẑ direction with magnitude B ext ≡ 2πnmag (dxN ) 2 . See [48] for details. At the same time, this requires a magnetic correction in the expression of the CS number like where N L cs in the rhs is the CS number given by Eq. (3.44) in the absence of an external magnetic field, and M 0 , M p are initial and final (after p time steps) magnetic field corrections, given by Let us notice that in the lattice, our choice of covariant derivative Eq. (3.30), mimics the continuum derivative D µ = ∂ µ − iA µ , without the gauge coupling e multiplying the gauge field A µ , contrary to our continuum definition introduced in Sects. 1, 2. Correspondingly, the lattice kinetic term of the gauge fields is divided by e 2 , see the lattice action (3.41). This implies that the electric and magnetic fields in the lattice should be converted to their continuum counterparts as E → e E cont , B → e B cont . This applies also the external magnetic field we introduced before, ext . This explains as well why we did not multiply by e 2 our lattice expression(s) for the CS number in Eqs. (3.44), (3.51).
For further details and clarifications about of the lattice formulation presented above, both in the absence and presence of a magnetic field, we refer the reader to Ref. [45].

Choice of parameters for lattice simulations
For the discretised Eqs. (3.43) to capture well the physics described by the continuum set of Eqs. (3.9)-(3.12), a number of conditions must be met in the lattice. In particular, to reduce the effects of discretisation, all relevant length scales of the continuum system λ i , must be larger than the lattice spacing dx. At the same time, to reduce finite volume effects, the same scales must also be smaller than the lattice size L = N dx, where N is the number of lattice points per dimension. In other words, relevant scales λ i must fulfill the condition In this section we discuss the choice of parameter values so that these two simple conditions are satisfied as best as possible.
Let us recall that we consider an effective potential for the Higgs-like field as with λ eff , m 2 eff the physical self-coupling and the physical mass of the scalar excitation. In the the simplest case of the Higgs phase of the theory, when the mass parameter is m 2 eff < 0 (at tree level m 2 eff ≡ −λv 2 < 0), we have two particle excitations: one corresponding to the Higgs boson with the mass m 2 H = 2λv 2 , and another to the vector boson with the mass M 2 A = e 2 v 2 , both given at tree level approximation. To minimize the number of scales, it is then natural to make the choice m 2 H = M 2 A , leading to 2λ = e 2 . We will keep this relation all through our investigation, even when considering the symmetric phase as well.
We are however mostly interested in the symmetric phase, when the effective mass parameter of the Higgs field is m 2 eff > 0. At tree level, the 3d Higgs model in this case contains a scalar excitation with mass m 2 eff and a massless photon. This requires In addition, we expect to have a length scale l π 2 e 2 1 T associated with the fluctuations carrying a CS-number N cs ∼ 1, see discussion in Section 1. This requires For instance, for e 2 1, dxT = 1 satisfies that T dx = 1 π 2 e 2 10 e 2 , and at the same time that LT = N π 2 e 2 10 e 2 (for as long as N > 10). The smaller we make dxT the larger we need N in the right hand side of the inequality (3.56). However, due to the need to simulate some cases up to very large physical times, very large values for N (typically expressed in powers of 2 n ) would require (unfeasible) long simulation times. Hence, as long as we take a value of e 2 neither much smaller nor much larger than ∼ 1, a good compromise is to fix the lattice spacing to the value dxT = 1. If instead of we choose typical values of the gauge coupling say as ∼ 10 −n with |n| 1 (with n either positive or negative), then we simply need to re-scale the lattice spacing to dxT = 10 n . From now on we stick to e 2 = 1 as a canonical value, and hence consider values for the gauge coupling somewhat smaller or larger than unity, but not orders of magnitude different. We then fix dxT = 1 for all our simulations. In that way we verify the left hand side inequality (3.56) with a factor of margin 1 10/e 2 , while the right hand side with a factor 1 10/(N e 2 ) 0.62/e 2 , 0.31/e 2 , 0.16/e 2 , ... for N = 16, 32, 64, ... respectively. If we literally take e 2 = 1, then N = 16 is really border line for verifying the rhs of inequality Eq. (3.56), as 1 is only marginally bigger than 0.62. Therefore, for N = 16 we will not simulate smaller gauge coupling values than unity, and hence e 2 ≥ 1. Similarly, if we double the points per dimension to N = 32, we must consider gauge coupling values e 2 ≥ 0.5, for N = 64 then e 2 ≥ 0.25, and so on. We have thus considered a sample of gauge coupling values taken by powers of 2, i.e. e 2 = 2 p , p = ..., −2, −1, 0, 1, .... In practice, N = 64 is the maximum number of points per dimension we take (otherwise the simulation time becomes too large), so the smallest coupling value we have really considered is e 2 = 2 −2 = 0.25.
At the same time, and independently of N , the gauge coupling cannot be very large (given that we have already fixed dxT = 1), because otherwise the lhs of inequality (3.56) becomes only roughly satisfied. Since we are taking the gauge coupling values in powers of 2, e 2 = 2 2 = 4 makes already too rough the verification of the lhs of the inequality. Hence, in practice, the maximum value we have considered is e 2 = 2 1 = 2. In summary, for the lattice spacing dxT = 1 we chose, we can capture well the scales of the fluctuations carrying out a CS-number of the order ∼ 1, for the parameters N = 16 : e 2 = 1, 2 (3.57) We need of course that such choice of parameters verifies as well the inequality (3.55). Some care is however needed when identifying in the lattice the Higgs mass parameter in the symmetric phase. The 3d Higgs model is a super-renormalisable theory, with the only parameter requiring renormalisation being the Higgs mass. As is immediately seen by the use of power counting, the Higgs mass is linearly divergent at the one-loop level and log-divergent at the two loop level. The relationship between the parameters in the M S renormalisation scheme m 2 (µ * ) and the lattice m 2 lat has been worked out in [49], and reads where µ * is the scale parameter of dimensional regularisation, the linear term ∝ e 2 , λ on the rhs is coming from one loop diagrams, and the terms proportional to ∝ e 4 , λe 2 , λ 2 from the two-loop diagrams. The latter can be neglected if m 2 (µ * ) e 2 4π 2 T 2 . Even though this condition is not really necessary, we will adopt it to make our analysis more transparent. At the same time, the relation between the physical and the lattice expectation value of the Higgs squared amplitude (which represents the 'order parameter' in the phase transition), was also computed in [49] as The non-perturbative lattice simulations of [50][51][52] have demonstrated that the physical correlation length in the symmetric phase is close to the M S mass m 2 (µ * ). So, in order to be in the symmetric phase but still sufficiently close to the phase transition, the lattice mass m 2 lat should be tuned to the linear counter-term ∝ 1 dx , to get the physical mass m 2 (µ * ) positive but sufficiently small. At the same time, we must satisfy the condition (3.55). In practice, we parametrize the physical mass as m 2 (µ) = y * e 4 T 2 , and choose the value of the dimensionless parameter y * , according to the phase diagram of the 3d Higgs + U(1) model [50][51][52]. We decide in this way whether we are in the symmetric or in the broken phase, depending on the value of |φ| 2 . Given our chosen values for e 2 , and based on the numerical outcome we obtain for Eq. (3.61), the typical values required to be in the symmetric phase, but close enough to the phase transition, are y * ∼ 0.01 − 0.05. This implies that e 2 T m −1 (µ) ∼ 1/ √ y * ∼ 5 − 10, and hence the condition (3.55) becomes, modulo a factor ∼ 0.5−1, identical to condition (3.56). Therefore, our set of {e 2 , N } values listed above in Eqs. (3.57)-(3.59), essentially satisfy both inequalities (3.55) and (3.56) in a reasonable way, given our choice dxT = 1. We will therefore stick to such parameter values in the following sections.
As a final remark about parameters let us note that we have worked in all simulations a time step dt sufficiently small, so that when decreasing it further, the results would be the same within the statistical error. Whereas dtT = 0.1 was borderline for some cases, dtT = 0.01 proved to be a good choice in all the cases. We thus adopted dtT = 0.01 as a fiducial value, and only increased it (hence speeding our runs) when we saw explicitly (by trial and error) that the observables were not sensitive to it.

Initial conditions
In Section 4 we start our numerical study by considering the diffusion of the Chern-Simons number in the absence of fermions, i.e. fixing µ = 0. In order to set up the initial thermal configuration, we use a standard Monte-Carlo Metropolis-Hasting algorithm to create a set of configurations according to the Gibbs distribution exp(−H L /T ), where H L is the lattice Hamiltonian The updates are such that we first update, in a given site, a given field amplitude, and then we determine whether we accept the change or not (if not, we undo the update). Following, we update in the same site, another field amplitude, and determine whether we accept it or not. We repeat this procedure in the same lattice site until we have considered all field amplitudes {φ, A}. Afterwards, we do the analogous updates, still in the same site, of the conjugate momenta {π, E}. After we have updated, in a give site, all field and conjugate momenta variables {φ, A, π, E}, then we move to the next lattice site, and repeat the procedure. We perform the updates at sites selected sequentially, i.e. at each sweep we survey the N 3 lattice sites always in the same order. By tuning (by trial and error) the typical size of the field and conjugate momenta updates, we achieve a ∼ 50% acceptance rate, and typically within O(10 2 ) − O(10 3 ) sweeps, we obtain the desired configuration. We consider both the case when an external magnetic field is absent (Sect. 4.1) and present (Sect. 4.2). In the case when there is an external magnetic field, say in theẑ direction B ext = Bê z , we simply repeat the procedure just described, but imposing the condition (3.50), whenever evaluating the magnetic field energy density in the lattice Hamiltonian (3.62).
The configurations created with the Monte-Carlo method do not preserve, however, the Gauss law G ≡ i (E i −E i,−i )−2e 2 dxIm{ϕ * π} = 0. As this is a constraint between the field amplitudes and conjugate momenta, the Metropolis-Hasting updates does not know about this apriori. One way to enforce the Gauss law during the Monte-Carlo process, is to perform the updates with respect to a new lattice Hamiltonian density H L → H L + ξG 2 , with ξ some dimensionless parameter. By tuning ξ to increasingly large values, one can obtain configurations that respect the Gauss law to an increasingly better and better degree. However, the larger the value of ξ, the larger the Monte-Carlo the thermalization time is required, making the method impractical if one desires to obey the Gauss law to a high degree of accuracy. In order to achieve a better accuracy we have used instead a cooling procedure, immediately after the the Monte-Carlo sweeps based on the Hamiltonian (3.62). In particular, we have solved the equations of motion derived from considering an effective action S cooling = G 2 . Such a cooling evolution hardly changes the energy of the gauge field transverse modes if the configuration was generated by the Metropolis-Hasting algorithm with Hamiltonian (3.62), but it kills the longitudinal components responsible for the violation of Gauss law [37]. In this way we can get the Gauss law satisfied to any desired accuracy for the initial configuration (most of our simulations reached machine precision). After the cooling, we finally start the real time evolution of the fields according to the set of discrete lattice equations of motion Eqs. (3.43) [with the chemical potential fixed to µ = 0].
In Section 5 we finally study the dynamics of the chiral charge µ when an initial value µ o = 0 is provided. As the classical equilibrium state of the system in the limit of large volume V → ∞ contains no chemical potential, we proceed as follows in order to study the non-equilibrium evolution of µ. We first prepare an equilibrium configuration with identical procedure as described above in the absence of µ, with Hamiltonian given by (a discrete version of) Eq. (3.13). Then we apply the cooling procedure, and following, before we start iterating the discrete equations of motion, we introduce µ = µ o = 0 as an initial condition. We then evolve the scalar and gauge fields, together with µ, with the use of Eqs. (3.43). We analyze the evolution of the chemical potential in Sect. 5, considering again both the absence (Sect. 5.1) and presence (Sect. 5.2) of an external magnetic field.

Chern-Simons number diffusion
In this section we consider the diffusion of the CS number in the absence of fermions, hence fixing µ = 0 at all times. We first consider the case with no external magnetic field in Sect. 4.1, where we expect the correlator Q 2 to become constant. In Sect. 4.2 we then 'switch on' an external magnetic field and study the diffusion of Q 2 , which now grows continuously as time goes by. In both cases, without and with external magnetic field, based on the numerical outcome from the lattice simulations, we obtain fits to the appropriate coefficients characterizing the behavior of Q 2 in each regime. We then compare with the analytical estimates presented in the previous sections.
Let us note that all through the sections where we report our numerical results, i.e. Sects. 4-5, we will consider the quantity Q defined as which is a factor 1/4 smaller than the CS number introduced in the theoretical discussions in Sects. 1-2. As explained in Sect. 3.2, given our lattice definition of the covariant derivative, we do not include the factor e 2 multiplying in the expression (4.1), as our lattice gauge fields (and hence the electric and magnetic fields) are a factor e bigger than the continuum counter parts. See discussion at the end of Sect. 3.2.

Zero magnetic field
In this section we show the results obtained from our numerical simulations in the absence of fermions and when we consider no external magnetic field (i.e. we fix n mag = 0, see Sect. 3.2 for the definition of magnetic number n mag ). We have obtained Q 2 for a fixed lattice spacing T dx = 1, varying the gauge coupling e 2 and the number of points per dimension N [hence varying the lattice volume 6 V = N 3 dx 3 ]. For each pair of parameters (e 2 , N ) we have obtained N R = 20 independent realizations of Q 2 as a function of time.
We have built the desired correlator Q 2 simply by statistical averaging, where the index i labels the different realizations. The statistical error in the measurement of Q 2 is determined by where Σ 2 Q 2 represents the statistical variance of the ensemble of N R realizations. In most plots we represent with a central solid line the mean value Q 2 , whereas a shaded area between some upper Q 2 +∆ Q and and lower Q 2 −∆ Q curves around a given central line, represents the error in the assessment of Q 2 . See for instance Fig. 1. We will maintain this plotting style all through the paper.
In the left panel of Fig. 1 we plot Q 2 for the gauge coupling e 2 = 2 and varying the volume V ∝ N 3 with N = 16, 32, 64 (T dx = 1 fixed). As expected in the absence of external magnetic field, the correlator saturates rapidly to a constant asymptotic amplitude. The saturation regime is reached rather quickly, in a time t sat 5T −1 , independently of N . We observe a monotonic growing dependence of the asymptotic amplitude of Q 2 with N . Based on the theoretical calculation Eq. (3.28), we expect the amplitude after saturation to scale linearly with the volume, Q 2 ∝ V ∝ N 3 . Hence, if we plot Q 2 V from the same data used in the left panel of Fig. 1, we expect the horizontal asymptotes of this quantity to sit on top of each other. This is precisely what we observe in the right panel of Fig. 1. There we clearly appreciate that the amplitude of the (volume-normalized) correlators nicely overlap with each other, within the statistical error.
In the left panel of Fig. 2 we plot Q 2 V for a fixed volume (T dx = 1, N = 64), and various gauge couplings e 2 = 0.25, 0.5, 1 and 2. The asymptotic amplitude of the correlator in the saturation regime, exhibits as expected, a monotonic growth with e 2 . Simple scaling arguments [37] led to expect the correlator to scale quadratically in α ∝ e 2 , Q 2 ∝ α 2 ∝ e 4 , V e 4 , based on the same data used in the left panel. We observe again that within the statistical error, the expected scaling Q 2 ∝ e 4 is nicely verified.
In Fig. 3 we show the asymptotic amplitude of Q 2 V as obtained for all the gauge couplings and volumes considered. Looking at Fig. 1 and Fig. 2, it is noticeable that once the growth of Q 2 saturates, its amplitude remains constant but oscillatory around some value. In order to assign a specific value to these wiggly asymptotic amplitudes, we average the signal over arbitrary long times in the saturation regime, like (4.4) We can assign two uncertainties to this quantity. First, a statistical error obtained by identical time averaging procedure over the upper and lower bound signals, and secondly a dispersion due to the oscillatory behavior of the signal (4.6) We will assign as the error of the asymptotic amplitude in the saturation regime, the maximum of the previous two errors, The data points and error bars in Fig. 3 represent the amplitudes Q 2 sat ± ∆ sat for each pair of values (e 2 , N ). Let us notice that even though the value Q 2 sat depends on the time t, see Eq. (4.4), the changes are always much smaller than the typical error ∆ sat , as long as t t sat . We also note that for the majority of the cases the statistical error Eq. (4.5) is larger than the oscillatory error Eq. (4.6), but not always.
In the left panel of Fig. 3 we show Q 2 sat V versus N . As expected, for each gauge coupling value (e 2 = 2, 1, 0.5 and 0.25), the amplitude is the same, within the errors assigned. Solid horizontal lines indicate the best fit to Q 2 sat V for a fixed value of e 2 , whereas dashed lines indicate the expected amplitude based on the theoretical scaling Q 2 sat V ∝ e 4 (taking the e 2 = 1 amplitude as a reference). In the right panel of Fig. 3 we show Q 2 sat V as a function of e 2 . Consistently with the left panel, the amplitude is the same (within the error) for a given gauge coupling e 2 , independently of the volume. We observe that the theoretical expectation Q 2 V ∝ (e 2 ) 2 is verified well within the errors, though the best fit of the data actually prefers a scaling as Q 2 V ∝ (e 2 ) 1.93 . The solid line in the figure represents an average of the best fits to each choice of N , whereas the shade area represents the dispersion among those best fits. In particular, for each case of N we obtain a fit to the data in the form 6.06 · 10 −6 · (e 2 ) 2.02 , N = 16 and We expect our results to agree well with the theoretical expectation quoted in Eq. (3.28), which can be re-written as Q 2 e 4 V = 1 16 (ζ(2)−ζ(3)) 2 3 π 6 3.6·10 −6 for dxT (30/π 2 ) 1/3 1.45 [37]. Our fit Eq. (4.9) was however obtained for dxT = 1, so the comparison cannot be made exact, and yet it only differs by a factor ∼ 2 with respect to the theoretical expectation.

Non-zero magnetic field
Now we turn our attention into the case where we consider a constant external magnetic field. As we explained in previous sections, without loss of generality we can consider a magnetic field in theẑ direction, B = (0, 0, B), with a flux at a given orthogonal lattice (4.14) We recall once again that, due to our choice of lattice gauge derivative [see Eq. (3.30)], this lattice magnetic field is related to the continuum one as In all the figures of this subsection, Figs. 4-10, we plot Q 2 and its statistical error ∆ Q 2 , as obtained from our simulations with an external magnetic field with magnetic seeds n mag = 4, 16 and 64. For each set of values (e 2 , N, n mag ) we obtain Q 2 ± ∆ Q 2 from N R = 20 independent realizations. Let us notice that due to the scaling of the magnetic field as B ∝ N −2 , the strength of the magnetic field weakens for larger the N (for a fixed magnetic seed n mag ). Hence, if we take e.g. n mag = 4, the magnetic field for the simulations with N = 32 will be a factor 1 4 weaker than the magnetic field for the simulations with N = 16 for the same magnetic seed. However, the strength of the magnetic field will be the same for N = 16 and N = 32 simulations if we take n mag = 4 for the former, but n mag = 16 for the latter.
In Fig. 4 we plot Q 2 as obtained from our simulations with n mag = 4. In the left panel we show Q 2 V for a fixed volume (T dx = 1, N = 16) and two gauge coupling values e 2 = 1, 2. As expected the correlator amplitude grows linearly in time as ∝ t, corresponding to the diffusive regime of CS number in the presence of a constant external magnetic field. The system actually attempts first to relax into a constant asymptotic amplitude (similarly as in the case of absent magnetic field), until the diffusion becomes noticeable. For such small volume simulations with N = 16, at a time t ∼ O(1)T −1 the system has already entered into some kind of diffusion. However it is not until a time t ∼ O(100)T −1 (gauge coupling dependent), that Q 2 enters into its final asymptotic regime, linearly growing in time. We have fitted this growth (averaging over the wiggly behavior) as Q 2 V = Γt, with Γ 1 7.11 · 10 −7 T 4 for e 2 = 1 and Γ 2 2.54 · 10 −6 T 4 . As it is appreciated in the plot, the linear behavior (once started) is well sustained during all our simulation time, and it seems a robust feature against the oscillatory behavior of Q 2 . The amplitude of the slopes is clearly bigger the larger the coupling e 2 , however the ratio Γ 2 /Γ 1 3.57 deviates from the  Figure 5. Growth of Q 2 in the presence of an external magnetic field with n mag = 4, for a variety of gauge couplings and volumes. Left: Volume and magnetic field normalized correlator, Q 2 V B 2 , for a fix coupling e 2 = 1, sampling the volume V = (N · dx) 3 for N = 16, 32 and 64. The straight line corresponds to a fit to the data, denoting a linear growth well sustained once diffusion is onset. Right: Volume, magnetic field and time normalized correlator, Q 2 V B 2 t , shown at times t t diff , once diffusion has been well established for all the cases. The horizontal lines correspond to fits to the data.
theoretical value 4, based on the expected behavior Γ ∝ (e 2 ) 2 (for fixed n mag ). We will quantify this deviation later when exploring other volumes and gauge coupling values, and taking into account statistical and averaging errors.
In the right panel of Fig. 4 we plot Q 2 V from our simulations with n mag = 4, in this occasion for a volume with N = 64 (T dx = 1), and the gauge coupling values e 2 = 0.25, 0.5, 1 and 2. The volume in the simulations used in the right panel is a factor (64/16) 3 = 2 6 = 64 larger than in the simulations in the left panel, so the physical magnetic field is, correspondingly, a factor 1 16 weaker (n mag = 4 is the same for both panels). This fact can be clearly appreciated by noting that the moment when the correlator amplitude starts growing is visibly delayed in the right panel, compared to the smaller volume case of the left panel. We quantify the onset of diffusion time scale later on. Let us now simply note that Q 2 also grows in time in the larger volume simulations, and its amplitude is larger the bigger the coupling e 2 . However, it is not clear at this point, whether the temporal growth is exactly linear, or what is exactly the scaling power with e 2 . As a guidance to the eye, we have added in the figure four straight lines corresponding to a linear growth with different slopes. These lines do not correspond to fits to the data, but rather represent a guidance to the eye for comparing the growth of Q 2 against a linear behavior.
In Fig. 5 we plot again the correlator Q 2 in the presence of a external magnetic field with n mag = 4, for a variety of gauge couplings and volumes. In particular, in the left panel we plot the volume and magnetic field normalized correlator amplitude, Q 2 V B 2 , for a fix coupling e 2 = 1, varying the volume V = (N · dx) 3 for N = 16, 32 and 64 (T dx = 1). Based on theoretical considerations, we expect the correlator to be directly proportional to external magnetic fields as Q 2 ∝ B 2 . We then expect the slopes of Q 2 V B 2 for different volumes N , to sit on top of each other. This is precisely what we observe in the left panel of Fig. 5, within the statistical error, after a time t O(10 3 )T −1 . As n mag = 4 is the same for all cases depicted, but N is different, the magnetic field B ∝ nmag N 2 is weaker the larger the volume simulated. Hence, the diffusion process starts later the weaker it is the external magnetic field B (i.e. the larger the N ). We find that diffusion is established at times T t diff 40, 400 and 4000 for N = 16, 32 and 64, corresponding respectively to a magnetic field B 16 = π 32 T 2 0.0982 T 2 , B 32 = 1 4 B 16 0.0245 T 2 and B 64 = 1 16 B 16 0.0061 T 2 . The dependence on the magnetic field of the time scale to onset diffusion, can be fitted like In summary, we can clearly appreciate two effects here, 1) the scaling Q 2 ∝ B 2 is pretty much verified by our simulations, and 2) the time scale to onset diffusion is larger the weaker the magnetic field (i.e. the larger the volume for a given magnetic seed n mag ). In the right panel of Fig. 5 we plot Q 2 V B 2 t for times t t diff , i.e. once the diffusion behavior is well established. If, as expected, the correlator grows linearly in time as Q 2 ∝ t, then Q 2 V B 2 t should turn into horizontal plots, lying on top of each other for a given gauge coupling strength e 2 , independently of the magnetic field B(n mag , N ) and the volume V (N, dx). This is precisely what we observe the right panel of Fig. 5, within the statistical error. To quantify deviations from the linearly growth behavior, we could introduce a new parameter Q 2 ∝ t 1+α . However, our attempts do so, show clearly that, within the statistical error, one can take safely take α = 0 for every set of values of (B, V, e 2 ) considered.
Let us also note that Q 2 grows with the gauge coupling. Based on the same scaling arguments as before, we expect Q 2 V B 2 t ∝ (e 2 ) 2 . In order to quantify possible deviations from this scaling, we have proceed similarly as in the absence of magnetic field, taking into account in this occasion that Q 2 grows in time. We have defined an amplitude for the correlator in diffusion as and assigned two uncertainties to this quantity, a statistical error as and an error due to the oscillatory behavior as We will assign as the error of Γ diff in the diffusion regime, the maximum of the previous two errors, ∆Γ diff ≡ max{∆Γ stat , ∆Γ osc } . Note that the factor 16 introduced in the definitions of Γ diff and ∆Γ diff is simply to compensate for the fact that we define Q as a factor 1/4 smaller than the standard N CS , see Eq. (4.1). This will make easier to compare later on the numerical rates with the theoretical predictions.
In the left panel of Fig. 6 we present Γ diff ± ∆Γ diff for all our data sample for n mag = 4 (i.e. considering all our volume and coupling constants pool). For each volume N , we have found a best fit to the data points In order to make connection with the continuum, given that B cont = B/e, we can also express the fits as Γ cont , as this will reveal the real dependence with respect to the gauge coupling. We obtain We have then averaged over the different fits, weighting each case by the number of gauge coupling values for each N [similarly as in Eq. (4.9) in the absence of a magnetic field]. We find an averaged fit of the form Γ ≡ Q 2 V t = (Ā ± ∆A) · (e 2 )m ±∆m , Γ diff (1.36 ± 0.07) · 10 −3 · (e 2 ) 2.92±0.14 · B 2 cont , (n mag = 4) , (4.22) where the errors in the amplitude and the exponent are obtained similarly as those in Eq. (4.9), in this case byĀ For a given gauge coupling value we find the amplitudes of this quantity to be similar to each other within the statistical error 7 , independently of the volume. Solid horizontal lines indicate the best fit to Γ diff B 2 for a fixed value of e 2 , whereas dashed lines indicate the expected amplitude based on the theoretical scaling Γ diff B 2 ∝ e 4 (taking the amplitude 7 The case (e 2 , N ) = (2, 32) shows however some stronger deviation.  Figure 7. Correlator Q 2 /V for a fixed gauge coupling e 2 = 1, sampling different volumes with N = 16 and 32, and magnetic fields with seed n mag = 4, 16 and 64. For comparison we also show Q 2 /V for e 1 and n mag = 0, for two volumes N = 16, 32. The straight lines are just meant to guide the eye (they are not fits), representing a linear growth ∝ t with different slopes for each case.
for e 2 = 1 as a reference). In the right panel of Fig. 6 we show Γ diff B 2 as a function of e 2 . Consistently with the left panel, the amplitude is of the same order for a given gauge coupling e 2 , independently of the volume. The solid line in the right panel represents the average fit Eq. (4.22) from the best fits Eqs. (4.21) to each choice of N . The shade area between the dashed lines represents the dispersion among such best fits. We conclude that the theoretical expectation Γ diff B 2 cont ∝ (e 2 ) 3 (equivalently Γ diff B 2 ∝ (e 2 ) 2 ) is verified well within the associated errors. The best fit of the data actually prefers a scaling as Γ diff  Let us note that we have introduced so far the B 2 factor in the fits Eqs. (4.21), (4.22) by hand, based on the theoretically expected scaling Q 2 ∝ B 2 . In other words, our fits above are obtained over the normalized correlators Later on, we will quantify numerically the deviations of our data from the theoretically expected behavior ∝ B 2 . Before, we shall present a similar analysis to the data obtained from our simulations with magnetic seeds n mag = 16 and 64 (i.e. with a magnetic field a factor 4 and 16 times larger than the case just analyzed so far).
In Fig. 7 we show for a fixed gauge coupling e 2 = 1, Q 2 /V for different volumes (N = 16, 32) and magnetic seeds (n mag = 4, 16, 64). If the scaling behavior Q 2 ∝ B 2 is correct, we expect some of the signals to overlap, as it is clearly appreciated in the figure. This is because the magnetic field scales as B ∝ n mag /N 2 , and hence, theoretically, the correlator scales as Q 2 /V ∝ B 2 ∝ n 2 mag /N 4 . Thus, the signal for a given pair of parameter values (n mag , N ) = (n 1 , N 1 ), should be, in principle, of the same size as the signal for other parameter values (n 2 , N 2 ), for as long as n 2 /n 1 = (N 2 /N 1 ) 2 . Our set of parameter values {e 2 , N } allow precisely for these type of combinations, for instance when comparing the signal obtained for (n mag , N ) = (16,16) [green signal in the figure] versus that for (n mag , N ) = (64, 32) [yellow signal], or the signal for (n mag , N ) = (4, 16) [dark blue signal] versus that for (n mag , N ) = (16,32) [light blue signal]. Within the statistical error these two pairs clearly overlap, once the diffusion regimen is well onset. As a guideline for the eye we also include straight lines depicting a linear growth, demonstrating that a linear growth in time remains a robust feature even in the presence of a strong magnetic field. For comparison, we also show the signal discussed in Section 4.1 in the absence of magnetic field (n mag = 0), for e 2 = 1 and for two volumes (N = 16, 32). In this case the growth of Q 2 /V ceases very soon, reaching a given amplitude fixed by the strength of e 2 = 1, independently of N . It is illustrative to note in Fig. 7 the growth of the correlator Q 2 /V for the parameters (n mag , N ) = (4, 32), which corresponds to the weakest magnetic field case of out sample [purple signal in the figure]. This correlator actually saturates initially to a similar amplitude to the signal with no external magnetic field, and only after quite some time ∆t ∼ 300 T −1 , it finally enters into a diffusion regime with its amplitude growing linearly in time, due to the presence of the weak magnetic field. In Fig. 8 we present a series of panels where we plot the correlator Q 2 for various magnetic seeds (n mag = 4, 16, 64) and gauge coupling values (e 2 = 0.5, 2), for a fixed lattice size N = 32 (T dx = 1). The color coding is actually shared by the four panels of the figure, as the plots of each panel are based in the same data. This way we can compare easily the effect of the different normalization of the correlators by simply comparing the signals of the same color from panel to panel. In the top left panel we show Q 2 V B 2 . If the scaling Q 2 ∝ B 2 is roughly correct, the six signals considered should be split in two groups, a higher amplitude corresponding to the highest gauge coupling e 2 = 2, and a lower amplitude corresponding to e 2 = 0.5. The top left panel clearly shows this two-group splitting pattern, indicating that a scaling ∝ B 2 is roughly correct [we will quantify this shortly]. In the top right panel of we plot Q 2 V e 4 . Based on the scaling Q 2 /B 2 ∝ (e 2 ) 2 , we expect the signals to split in three groups, depending on the strength of n mag (as N is common to all of them). For each magnetic seed n mag = 4, 16 and 64 we should see a pair of overlapping correlators. This three-group pattern is clearly observed in the top right panel. In the bottom left panel we plot Q 2 V B 2 e 4 . Based on the scaling Q 2 ∝ B 2 (e 2 ) 2 , we expect the whole set of signals to overlap into a single growing pattern. This is precisely what we observe in this panel, where all signals overlap within a scatter of ∼ 1 order of magnitude, compatible with the statistical errors. For further clarity, we plot again these set of overlapping signals in the bottom right panel, but this normalizing by a linear time growth as well, Q 2 V tB 2 e 4 . If the scaling Q 2 ∝ B 2 (e 2 ) 2 t is correct, then all the signals should not only overlap but also remain constant in time. Within the statistical error, this is precisely the pattern we observe in the bottom right panel, where all signals oscillate and scatter around a common value ∼ 10 −4 .
We have analyzed our n mag = 16 and n mag = 64 data similarly as we did for n mag = 4. For each pair of values (n mag , N ) we have determined Γ diff ± ∆Γ diff according to Eqs.
(4.20). Then, for each magnetic seed n mag , we have found a best fit to the data points (e 2 , Γ diff ) for each volume N , of the form Γ For each n mag value, we have averaged over the different fits, weighting each case by the number of gauge coupling values for each N (similarly as we did for the n mag = 4 data). We find an averaged fit of the form Γ diff ≡ Q 2 V t = (Ā ± ∆A) · (e 2 )m ±∆m as [for completeness, we reproduce again the fit to n mag = 4 Eq. In the left panel of Fig. 9 we present Γ diff ± ∆Γ diff as a function of e 2 , using all our data sample for n mag = 16 (i.e. considering all our volume and coupling constants pool). In the right panel of Fig. 9 we present the analogous plot but for all our data sample for n mag = 64. The solid line in the left panel represents the average fit Eq. (4.26) from the best fits Eqs. (4.23) to each choice of N . In the right panel of Fig. 9 we present the analogous plot but for all our data sample for n mag = 64. The solid line in the right panel represents the average fit Eq. (4.27) from the best fits Eqs. (4.24) to each choice of N . The shade area between the dashed lines in both left and right panels, represents the dispersion among the best fits for each n mag case. We conclude from these fits that the theoretical ∝ (e 2 ) 3 ) is still well verified within the associated errors for both n mag = 16 and n mag = 64. The best fit of the data for n mag = 16 prefers a scaling as Γ diff B 2 ∝ (e 2 ) 1.99 , whereas for n mag = 64 it prefers Γ diff B 2 ∝ (e 2 ) 1.76 , see Eqs. (4.26)-(4.27). The uncertainty however in the exponent encompasses well, in both cases, the theoretical scaling.
Finally, let us return to the fact that in all our previous fits as a function of e 2 , we have always assumed that the scaling Q 2 ∝ B 2 was exact. Even though our plots in Fig. 8 showed clearly that, if not exact, such scaling must be roughly true, we can of course do much better, test this simple scaling behavior against our data. In order to do so, we have now divided our data sample in such a way that, for a given gauge coupling value e 2 , we collect the amplitudes Γ diff as a function of volume N and magnetic seed n mag . Given a given pair of values (e 2 , N ), we can then fit Γ diff as a function of n mag , theoretically expecting that Q 2 ∝ n 2 mag . For each coupling e 2 we have found a best fits for each volume N considered Γ diff (4.50 ± 0.11) · 10 −4 · B 2.08±0.08 , e 2 = 0.5 (4.33) We note that there is clearly some deviation with respect an exact ∝ B 2 scaling. In the case of e 2 = 2.0 and e 2 = 0.5, the deviation is not significant, whereas for e 2 = 1.0 we observe a systematic deviation of f ew % from an exact quadratic scaling. Let us note that whereas the fits Eqs. (4.25)-(4.27) of Γ diff vs e 2 are obtained for a fixed magnetic seed n mag , the new fits Eqs. (4.31)-(4.33) of Γ diff vs B are obtained for a fixed gauge coupling e 2 . Another approach we may consider is therefore a multi-dimensional fit to Γ diff , as function of both e 2 and B. Theoretically Γ diff ∝ B 2 (e 2 ) 2 , so we can consider our data sample of Γ diff values obtained for each set of parameter values (e 2 , N, n mag ), and attempt a multi-dimensional fit of the whole array of data to a functional form Γ diff = A · (e 2 ) me (2πn mag ) m B N −2m B +3m N , where m e will characterize deviations from ∝ (e 2 ) 2 , m B will characterize deviations from ∝ B 2 , and m N will characterize deviations from Q 2 ∝ V . Using our full set of data obtained for N = 16, 32, 64, e 2 = 0.25, 0.50, 1.0, 2.0 and n mag = 4, 16, 64, and expressing the fit in terms of the continuum magnetic field, we obtain 8 with the errors representing 68% confidence level intervals. Our multi-dimensional fit indicates that there is no residual volume dependence, so that Q 2 ∝ V is pretty much exact, whereas the scalings Q 2 ∝ B 2 cont and Q 2 ∝ α 3 are pretty much compatible with our data set, within the statistical errors. As the theoretical estimation scales exactly as ∝ α 3 ∝ e 6 , and ∝ B 2 cont , in order to make a fair comparison with our numerical outcome, we present below the numerical diffusion rates obtained when assuming an exact scaling of the data as ∝ e 6 B 2 cont , Comparing the prefactor from Eq. (4.39) with that in Eq. (4.35), we see that we obtain numerically a rate Γ (58 ± 9) times larger than the theoretical expectation 9 . This is of course one of the most relevant results of this paper.

Dynamics of chiral charge
Let us now turn our attention into the case where a non-zero fermion charge density is initially present in the system, characterized by a chemical potential µ, with initial value µ o . As discussed in Sect. 3.1, in a cubic lattice with length L = N dx, we there is a gauge coupling and volume dependent critical value of the chemical potential, Eq. (3.24). In the absence of magnetic field given our choice dxT = 1, it reads This value signals a threshold below which the chemical potential is not expected to be transferred anymore into the gauge field sector for B = 0, see Sect. 3.1 for details. The value of the chemical potential to see the instabilities is rather large, and for the relatively small lattices we used it induces the large perturbation of the system. Moreover, for µ ∼ µ c the number of unstable modes is quite small. In other words, we are quite far from the continuum limit and far from the physically interesting domain µ/T 1. Therefore, it is difficult to compare the results in this section with the expectations from the continuum theory. So, our study here is qualitative rather than quantitative. In Sections 5.1 and Section 5.2 we consider the behavior of the system with an initial chemical potential both when µ o > µ c and µ o < µ c , in the absence and presence of an external magnetic field, respectively.

Zero magnetic field
In the absence of an external magnetic field, if the initial chemical potential verifies µ o ≤ µ c , the chemical potential is expected to remain equal to its initial value at all times. However, for µ o > µ c , we expect the chemical potential to relax its amplitude down to a final value µ → µ c . The decay of the chemical potential from an initially large value µ o down to µ c is realized through exponentially damped oscillations of µ. The damping envelope of the oscillatory pattern of µ characterizes a decay rate, corresponding to an energy transfer of the initial fermion charge into long wavelength gauge fields. In Fig. 11  Given that initial condition, and since the critical value scales as µ c ∝ 1 e 2 , the smaller the value of the gauge coupling the larger the initial chemical potential 9 Of course, if we consider the multidimensional-fit Eq. (4.34), the numerical rate can be considered to be closer to the theoretical, given the great dispersion obtained in the multidimensional fitted amplitude. This large dispersion is however only a consequence of having fitted the data against multiple parameters, and thus we consider more fair to make a comparison between theory and numerics by using Eq. (4.39), extracted from the data once we assume (for the shake of the comparison) a scaling Γ diff ∝ α 3 B 2 cont as exact.  Figure 11. Relaxation of the chemical potential µ when its initial amplitude µ o is larger than the critical value µ c , given by Eq. (5.1). The chemical potential relaxes through damped oscillations to µ → µ c , as expected. In all plots we take µ o = 2µ c . Horizontal dashed lines represent either initial or final chemical potential amplitudes. is in Fig. 11 (for a given volume N ). Hence the larger the final chemical potential value to which the system relaxes to, after several oscillations. This can be clearly appreciated in each of the panels of Fig. 11, when comparing the plotted lines within each panel. At the same time, as the critical value scales with the volume as µ c ∝ 1 N , the larger the volume is, the smaller the initial chemical potential is (for a fixed gauge coupling), and hence the smaller the value the chemical potential relaxes to. This can also be seen in the output from our numerical simulations in Fig. 11, comparing the initial and final values of µ for a common gauge coupling value (e 2 = 1 or 2 in the figure), between the left and the right panels.
Let us also note, looking at Fig. 11, that the chemical potential remains constant and equal to its initial value for brief period of time, before the oscillations start. The time scale t dec for the onset of the damped oscillations is mostly independent of the value of e 2 . This can be clearly seen by the eye in the panels of Fig. 11, where t dec is very similar among the different cases within each panel. However t dec depends clearly on the volume N , as t dec for N = 32 is roughly double than for N = 16, as appreciated in Fig. 11 when comparing the left to the right panel.
In the left panel of Fig. 12 we show 10 different realizations of the relaxation of the chemical potential µ when its initial amplitude is larger than the critical value. All realizations are obtained for (e 2 , N ) = (0.5, 32) and an initial chemical potential µ o = 2µ c . The behavior of µ in all realizations is rather similar, e.g. starting their decay at a similar time t dec , or performing a similar number of oscillations until the system relaxes down to µ c . The exact details of the damped oscillations differ however from realization to realization, as the oscillations are not necessarily in phase, neither they reach necessarily the same height at local maxima. Due to this, extracting a mean value µ of the chemical potential averaged over realizations [for a given set of fixed (e 2 , N ) values], may not be a good estimator of its damping: different realizations may interfere whenever they are  , when its initial amplitude µ o is larger than the critical value µ c given by Eq. (5.1). In all plots we take µ o = 2µ c , and the damped amplitude relaxes to µ → µ c , as expected. Horizontal dashed lines represent the initial (top) and final (bottom) chemical potential amplitudes. Right: Extraction of local maxima of a given realization of the evolution of chemical potential when the initial amplitude µ o is larger than the critical value µ c . Green triangules represent all local maxima, blue squares all local maxima and the chemical potential at the onset of oscillations, and red circles are the same as the blue squares but removing those points followed by a higher maximum, so that maxima in this subset follow always a monotonic decaying pattern. out-of-phase, hence resulting in a pattern not following a smooth and monotonic damping envelope. This can be appreciated by the thick red line in the left panel of Fig. 12. Even though the averaged value µ still follows the general trend of damped oscillations, the pattern of the averaged oscillations exhibits a lower amplitude in general, and even a noticeable suppression at some of its maxima. For a larger number of realizations these effects may become even more noticeable. Therefore, in order to characterize the relaxation decay time of the system in our numerical simulations, we will fit the oscillations envelope in a realization by realization basis, and then average over the fits, rather than fitting the averaged µ function.
In the right panel of Fig. 12, we plot one realization of the relaxation of µ from an initial value larger than the critical one. We expect the system to decay exponentially fast down to µ µ c . Initially the chemical potential is frozen. After some time t * [t * 12 for the particular example in this figure, corresponding to (e 2 , N ) = (1, 16)], the chemical potential starts oscillating. To be precise, we extract the time scale of the onset of the oscillations t * by determining the time when the amplitude of µ has fallen a give small * fraction with respect it initial amplitude µ o . We chose * = 0.03, but note that our following results are rather insensitive (within the statistical error) to any choice of this initial threshold between * ∼ 0.01 − 0.05. We take the amplitude µ * = µ(t * ) ≡ (1 − )µ o as the 'onset-of-oscillations' point from where to determine the decay envelope of the chemical potential. We then determine all local maxima of the oscillations of µ, but removing the points followed by another maximum of higher amplitude. In general we can think of the chemical potential as a product of two functions, µ(t) = Υ(t) · f (t), with f (t) representing the quasi-periodic oscillations, and Υ(t) their damped envelope amplitude (i.e.Υ < 0). In the right panel of Fig. 12, green triangles represent all local maxima of the shown realization, whereas blue squares are equal to the green triangles plus the 'onset-of-oscillations' point (t * , µ * ). Red circles represent the same as the blue squares but removing those points with a local maxima followed by a maximum at a higher amplitude. The set of red circles represents therefore a selection of all maxima [including the 'onset-of-oscillations' point (t * , µ * )] which form a monotonic decreasing pattern of the chemical potential's envelope amplitude.
behavior of the chemical potential. We can then build statistically averaged quantities For each set of parameters (e 2 , N ) we characterize therefore the chemical potential decay by an envelope function Υ = µ c [1 +Ā exp{−Γ(t − t * )}], where the fitted parameters 10 admit an uncertaintyΓ ± ∆ Γ andĀ ± ∆ A . For µ o = 2µ c , we have µ * = 2(1 − * )µ c , so we expect the initial amplitude A r to be of order of (but somewhat smaller than) unity. If we forced our fits to pass through (t * , µ * ), then there would be no need to fit A r , as this will be given automatically by A r = (1 − 2 * ) for every realization 11 . We prefer however to fit our data leaving A r free, so that the fitted envelope amplitude considers the point (t * , µ * ) for the fitting analysis, but does not force the fitted curved to pass necessarily through such starting point. Putting all our fits together, we obtain Table 1, where the decay widthsΓ ± ∆ Γ are listed for each pair of values (e 2 , N ). The time scale of the decay is characterized by t dec ≡ 1/Γ. We consider the system has relaxed into its final asymptotic chemical potential value µ → µ c at a time ∼ (4 − 5)t dec , where the envelope amplitude has been damped already by a factor Υ µc − 1 ∼ O(10 −2 ). From the theoretical considerations discussed in Sect. 1, given that our initial chemical potentials are larger than the Debye screening (µ > T /e), we expect the decay width to behave according to the theory at zero temperature, i.e. asΓ ∼ k inst ∝ αµ o ∝ α 0 [see Eq. (1.6) for k inst definition], where in the last proportionality we have used that µ o = 2µ c ∝ 1/e 2 . We obtain however a behavior asΓ ∼ e 2 , as shown in Fig. 14 (2,32), because in that case the decay is very strong, and it cannot be well modelled 12 with an exponential decaying envelope function, see e.g. left bottom panel of Fig. 13. In continuum theory, one expects the rate to behave asΓ ∝ e 2 only when the chemical potential is sufficiently small µ < e 2 T , as this guarantees that the instability scale 1/k inst is larger than the fermion's mean free path λ ∼ 1/(α 2 T ), see Sect. 1. This correspond to the magneto-hydrodynamical regime, but we are far from it. We make no attempt to understand these discrepancies and attribute them to the lattice artifacts and to essentially non-linear character of our system at large chemical potentials.

Non-zero magnetic field
Let us finally consider the most complex case when there is an initial chemical potential, in the presence of an external magnetic field. Let us consider first the case when the initial chemical potential µ o is larger than the critical value µ c , recall Eq. (5.1). As just described in Sect. 5.1, in the absence of magnetic field the chemical potential will relax down to µ c , via damped oscillations, with decay width Γ depending on the parameter values (e 2 , N ). However, we now consider the presence of a background magnetic field given by B = 2πn mag /(dxN ) 2 = eB cont , recall Eq. (4.14). Hence, we also expect that after a while, due to the external magnetic field, the chemical potential will tend to decay to zero 12 In such a case the chemical potential is very rapidly damped down to µc, and hence a fitting in the form of sustained damped oscillatory pattern is not a good description, as the chemical potential relaxes into the final value only after 1-2 oscillations. exponentially. Besides, from consistency, we expect the exponent of this magnetic−induced decay to be related to the corresponding diffusion rate measured in Section 4.2, for small enough chemical potentials. These types of relations have indeed been verified for non-Abelian theories by Moore [42]. In our study we did not attempt to do so because of the quite large chemical potentials we could only simulate numerically. In Fig. 15 we plot the evolution of the chemical potential from an initial super-critical value, showing results for all our volume and gauge coupling sampling. We consider as in the previous section µ o = 2µ c , without loss of generality. Let us begin by analyzing the weakest magnetic seed we considered, n mag = 4. In the top panels we show the chemical potential evolution for µ o = 2µ c , n mag = 4 for N = 16 (left panel) and N = 32 (right panel). Let us recall that the magnetic field scales as B ∝ n mag /N 2 . Hence, for the same n mag , the magnetic field is weaker the larger the volume N . This explains the main difference between the left and right top panels. In the top left panel, the magnetic field is B(n mag = 4, N = 16) = 0.0982 T 2 , which represents a value sufficiently strong to trigger the final 'magnetically induced' decay of the chemical potential down to zero. In this figure we distinguish two regimes: first when the chemical potential relaxes down to its critical value through exponentially damped oscillations (like like in Section 5.1, in the absence of a magnetic field), and secondly, when the chemical potential begins its exponential decay down to zero, a time afterwards it had reached the critical amplitude. This second magnetic-induced decay can be characterized by a fit with t B the starting point of this magnetically induced exponential damping. As the evolution of µ is shown in semi-log scale, the fits (when applicable) correspond to straight lines, see e.g. top left panel in Fig. 15. The decay rates Γ 5 measured in the top left panel in Fig. 15, and in general in all panels of this figure, are reported in Table 2.
In the right top panel of Fig. 15 we see the chemical potential relaxing to its critical value as in the absence of magnetic field, but then we observe that it remains constant at that value. To guide the eye, we added horizontal (dashed) lines, to emphasize that the chemical potential remains really constant after its initial relaxation down to µ c . The magnetic field in this case, B(n mag = 4, N = 32) = 0.0245 T 2 is too weak (4 times smaller than in the analogous left panel case), so the regime where the chemical potential is expected to decay from µ c down to zero, is never reached. Hence no fits are provided.
In the central panels of Fig. 15 we show the chemical potential evolution for initial super-critical value µ o = 2µ c , when n mag = 16, for N = 16 (left panel) and N = 32 (right panel). The magnetic field B in these central panels is now 4 times stronger than in the corresponding top panels. In the left central panel we observe a very strong decay of the chemical potential down to zero. Even though the chemical potential attempts to oscillate, presumably around µ c , the magnetic-induced decay is too strong, forcing the chemical potential to decay to zero almost immediately after the onset of the initial damping. In the right central panel, however, the magnetic field is 0.25 times weaker than in the left panel, so we observe a similar situation as in the top left panel. The chemical potential first relaxes down to the critical value as in the absence of magnetic field, and then it decays down to zero after a period of time (which depends, non-monotonically, on on the gauge coupling). The magnetically induced final decay rates of these central panels are reported in Table 2.
is the highest we have considered in our simulations, so strong that the chemical potential is immediately dumped to zero, skipping the relaxation phase down to the critical value. The inertia of the magnetic induced damping is so strong, that the chemical potential becomes negative and oscillates around zero. It is still possible, however, to fit the envelope amplitude of these oscillations with an exponential decay as in Eq. (5.5). We report the fits to these strong decay rates in Table 2. In the right bottom panel of Fig. 15, the magnetic field is like before, 0.25 weaker than in the analogous left panel. We observe nonetheless that the magnetic field is yet sufficiently strong to drive the chemical potential down to zero, again skipping an initial relaxation down to µ c . However, the effect is not as strong as in the bottom left panel, so that µ only becomes negative in the case of e 2 = 1.0, whereas it remains positive for e 2 = 0.5 and 2.0. It is therefore not monotonic in e 2 whether the system will be driven or not towards a chemical potential of negative value. Positive or negative, what is common to all gauge coupling values in this case, is that the system behaves as an over-damped oscillator, since after the initial descent of the chemical potential value [which can be fitted by an exponential decay with Eq. (5.5) by a brief period of time], the system does not oscillate afterwards, simply approaching zero from either positive or negative chemical potential values. The dumping of the chemical potential in all these cases can then be split in two decay rates, but only the first stage admit a fitting as Eq. (5.5), the decay rate of which is reported in Table. 2. Let us study now the case when the chemical potential is initially equal or smaller than the critical value. In Fig. 16 we show a series of panels plotting the evolution of the When the chemical potential is initially critical, µ o = µ c , we observe that for e 2 = 1, the chemical potential decays exponentially fast with a given decay rate Γ 5 during a finite period of time, and then changes its decay to a slower rate (yet exponential) Γ  Table 3. A similar 'double' decay rate behavior is observed for e 2 = 1 when the chemical potential is initially marginally sub-critical, µ o = 0.75µ c , although the corresponding decay rates Γ 5 are slightly different than those measured when µ o = µ c , see Table 3. The decay of the chemical potential for e 2 = 2 exhibits however a single exponential decay both when µ o = 1.0µ c and µ o = 0.75µ c , though with different damping rate depending on µ o . These decay rates are provided in Table 3. Finally, for smaller initial chemical potentials, µ o = 0.5µ c and µ o = 0.25µ c , the decay of µ is also characterized by a single exponential damping rate, independently of the gauge coupling strength. See Table 3 for the specific decay rate values.
We continue studying the case when the chemical potential is initially equal or smaller than the critical value, but this time increasing the strength of the external magnetic field. In Fig. 17 we show a series of panels plotting the evolution of the chemical potential in the  Table 3. Fits to magnetic field -driven decay rates of initially critical and sub-critical chemical potentials, in the presence of a physical magnetic field B = 0.0982 T 2 , corresponding to a magnetic seed n mag = 4 and volume N = 16. These fits correspond to the different cases plotted in Fig. 16.   Table 4. Fits to magnetic field -driven decay rates of initially critical and sub-critical chemical potentials, in the presence of a physical magnetic field B = 0.393 T 2 , corresponding to a magnetic seed n mag = 16 and volume N = 16. These fits correspond to the different cases plotted in Fig. 17. The upper bound in the time interval of application of these fits is only due to the fact that the data becomes too noisy, so it cannot be fitted by Eq. (5.5).
All four plots correspond again to a N = 16 lattice volume, so the physical magnetic field is B = 0.393 T 2 in all cases. When the chemical potential is initially critical, µ o = µ c , we observe that for both e 2 = 1 and e 2 = 2, the chemical potential decays extremely fast, turning negative in a brief interval of ∆t ∼ 70 T −1 . Once it becomes negative, it soon stops decaying, and starts approaching slowly and asymptotically to zero, while remaining negative. We cannot fit this decay as an exponential damping, not even while µ remains positive, as the chemical potential decreases simply faster than exponentially (until it bounces back once it has already become negative). A similar trend is actually followed for e 2 = 1 in all sub-critical cases, µ o /µ c = 0.75, 0.50, 0.25, with the chemical potential falling always down to negative values in a short interval of time ∆t ∼ 70 − 90 T −1 . The only difference among the sub-critical scenarios is that the smaller the initial chemical potential, the less negative the values reached by µ, and the longer it takes the slow asymptotic approach to zero from negative values. In the case of a larger gauge coupling e 2 = 2, the behavior of the chemical potential is however very different. For the sub-critical cases µ o /µ c = 0.75, 0.50, µ first starts a faster than exponential decay, but this is stopped after a short time (around the same moment when µ for e 2 = 1 turns negative) and from then on the chemical potential simply decays following an exponentially damped behavior like Eq. (5.5), the rates of which are reported in Table 4. For the most sub-critical case µ o /µ c = 0.25, µ does not exhibit any faster than exponential stage, and simply starts a exponential decay towards zero, the rate of which can also be fitted by Eq. (5.5), as listed in Table 4. We also study the case when the chemical potential is initially critical or sub-critical, in the presence of the strongest external magnetic field that we have simulated. In Fig. 18 we show again a series of panels plotting the evolution of the chemical potential, but this time in the presence of a magnetic field with n mag = 64, where as before we consider  This regime is sustained during several oscillations, before the chemical potential relaxes to zero. The envelope amplitude of damped oscillations can be fitted in each case by an exponential decay as in Eq. (5.5), the rates of which we report in Table 5 for the different initial chemical potentials. For the largest gauge coupling e 2 = 2, a similar pattern is developed. The amplitude of the chemical is however damped faster, exhibiting only only 2-3 semi-oscillations for µ o /µ c = 1.0, 0.75, 0.50, and even even less, 1-2 semi-oscillations, for the smallest initial chemical potential µ o /µ c = 0, 25. The envelope of these strongly damped oscillations for e 2 = 2 can also be fitted by an exponential decay as in Eq. (5.5), the rates of which are reported as well in Table 5 for all the initial chemical potential values. These last fits are however not as good as for e 2 = 1, and should be taken only as indicative of the time scales, particularly for the case µ o /µ c = 0.25, where the damping effects is so efficient that we barely barely capture 1 semi-oscillation of µ.
In Fig. 19, we show the evolution of the chemical potential when its initial value is critical or sub-critical, in the presence of magnetic field with seeds n mag = 4 or 16, but for a larger lattice volume N = 32 than before. The physical magnetic fields are therefore weaker by a factor 0.25 compared to those in the analogous plots in Figs. 16-17. In the upper left panel of Fig. 19 we show the results for n mag = 4, where the magnetic field B 0.0245 T 2 is relatively weak. The chemical potential in this case starts from the beginning a slow  Table 5. Fits to magnetic field -driven decay rates of initially critical and sub-critical chemical potentials, in the presence of a physical magnetic field B = 1.57 T 2 , corresponding to a magnetic seed n mag = 64 and volume N = 16. The rates tabulated correspond to the different cases plotted in Fig. 18, obtained via fitting the envelope function of the µ oscillations with Eq. (5.5).  Figure 19. Evolution of the chemical potential when the initial value is critical (µ o = µ c ) or subcritical (µ o < µ c ), in the presence of magnetic fields with seeds n mag = 4, 16. All plots correspond to a lattice volume N = 32, so the physical magnetic fields are weaker by a factor 0.25 compared to those in the analogous plots in Figs Table 6. Fits to magnetic field -driven decay rates of initially critical and sub-critical chemical potentials, in the presence of a physical magnetic field B = 0.0245 T 2 (n mag = 4) and B = 0.0982 T 2 (n mag = 16), all corresponding to a lattice volume N = 32. The rates tabulated correspond to the different cases plotted in Fig. 19, obtained via fitting of the envelope function of the µ oscillations with Eq. (5.5).
decay, yet exponential, independently of the initial condition for the chemical potential (as long as µ ≤ µ c ). The decay rate is faster the larger the gauge coupling, as reported in Table 6. A similar pattern of the decay of the chemical potential is repeated in the case n mag = 16 for the smallest initial chemical potential amplitudes µ o /mu c = 0.50, 0.25, as can be seen in the upper right panel of Fig. 19. The decay rates corresponding to these cases are also reported in Table 6. For n mag = 16 and critical or marginally critical initial values µ o /mu c = 1.0, 0.75, the evolution of the chemical potential is actually very similar to the (N, n mag ) = (16, 4) case (for the same initial critical or marginally critical initial condition), as the physical magnetic fields are the same. In particular, a double decay rate behavior is observed for the smallest gauge coupling e 2 = 0.5. The decays of the chemical potential for e 2 = 1.0 and e 2 = 2.0, exhibit however a single exponential damping rate both when µ o = 1.0µ c and µ o = 0.75µ c . For µ o = 0.50µ c and µ o = 0.25µ c , all gauge couplings e 2 = 0.5, 1, 2 exhibit also a single exponential damping rate. We report all decay rates of these cases in Table 6, based fitting the data with Eq. (5.5).
Finally we study the case when the chemical potential is either critical or sub-critical, for our larges magnetic seed n mag = 64, and for our largest lattice volume N = 32. This makes the external magnetic field weaker compared to the analogous (n mag , N ) = (64, 16) case. In Fig. 20 Figure 20. Evolution of the chemical potential when the initial value is critical (µ o = µ c ) or sub-critical (µ o < µ c ), in the presence of a magnetic field with seed n mag = 64. All plots correspond to a lattice volume N = 32, so the physical magnetic fields are weaker by a factor 0.25 compared to those in the analogous plots in Fig. 18. Continue lines correspond to our simulated data, whereas dashed lines correspond to fits.
All four plots correspond again to a N = 32 and n mag = 64, so the physical magnetic field is B = 0.393 T 2 in all cases. Let us first note that for the critical or marginally subcritical initial conditions, µ o /µ c = 1.0, 0.75, the chemical potential for the largest gauge coupling e 2 = 2, starts decaying initially faster than exponentially, until this behavior is stopped very rapidly (on its way down towards crossing zero) after a brief period of time (the shorter the smaller the initial chemical potential), while µ still remains positive. From then on the chemical potential simply decays following an exponentially damped behavior. For the subcritical initial conditions µ o /µ c = 0.5, 0.25, the chemical potential for e 2 = 2 simply follows an exponential decay from the beginning, with no initial fasterthan-exponential phase. All the rates of the exponential decays for e 2 = 2 are reported in Table X. For all µ o /µ c = 1.0, 0.75, 0.5, we observe that for the smallest gauge couplings e 2 = 0.5 or e 2 = 1, the chemical potential decays faster than exponential (similarly as in Fig. 17), turning negative in a brief interval of ∆t ∼ (50 − 100) T −1 . Once it becomes negative, it soon stops decaying, and starts approaching slowly and asymptotically to zero, while remaining negative. The smaller the initial chemical potential, the minimum (negative) value reached by µ is closer to zero , and the longer it takes the slow asymptotic approach to zero from negative values. In the case of µ o /µ c = 0.25, we observe that only  Table 7. Fits to magnetic field -driven decay rates of initially critical and sub-critical chemical potentials, in the presence of a physical magnetic field B = 0.393 T 2 , obtained for n mag = 64 and N = 32. The rates tabulated correspond to the different cases plotted in Fig. 20, obtained via fitting of the envelope function of the µ oscillations with Eq. (5.5). Upper bounds in the time interval of application of the fits is only due to the data becoming too much noisy.
for e 2 = 0.5 does the chemical potential decays faster than exponential (like in the cases with µ o /µ c = 1.0, 0.75, 0.5), whereas for e 2 = 1.0 µ follows an exponential decay from the beginning (see bottom right panel in Fig. 20).
To conclude this section, let us note that, qualitatively, the behavior of the system containing both a chemical potential and an external magnetic field, can be described by the system of Eqs. (2.10),(2.11), if we plug in the following ansatz for the vector potentials, which captures the development of the instability. The terms proportional to B describe the external magnetic field in direction of z-axis. The first terms in A 1 and A 2 correspond to maximally helical magnetic field with momentum k, which appears as a result of a instability in the presence of a chemical potential, whereas A 3 is associated with the electric field generated along the external B.
The field (5.6) leads to space-independent E · B = (Bḟ z −kfḟ ) and obeys the equations f = kf e 2 µ 4π 2 − k − σḟ (5.7) µ = 3e 2 T 2 π 2 (Bḟ z − kfḟ ) (5.8) For a sufficiently large chemical potential, the instability associated with the variable f drives the system into the state with chemical potential µ/T = 4π 2 /e 2 . For some time the system stays there, and eventually the chemical potential is diluted due to the generation of an electric field characterized by f z .

Conclusions
In this work we have studied the non-conservation of fermion/chiral number in Abelian gauge theories at finite temperatures. We have considered different mechanisms of fermionic charge disappearance in the plasma, and analyzed these with both analytical and numerical techniques. We use a bosonic effective theory where a non-zero fermionic charge density is characterized by a chemical potential µ, so that we add in the Hamiltonian a term H CS ∝ µQ, with Q ∝ d 4 xF µνF µν the Chern-Simons (CS) number.
In sections 1 and 2 we have provided a number of analytical results about the dynamical behavior of this system, both when µ = 0 and µ = 0, either in the presence or absence of an external magnetic field. In Section 3 we have discussed the details of our modeling, including its lattice formulation. We then dedicate Sections 4 and 5 to our numerical analysis, based on the outcome from the real time lattice simulations we performed.
In section 4 we have considered the diffusion of the CS number in the absence of fermions, hence setting µ = 0 as a fixed value at all times. We have studied initially the case with no external magnetic field in Sect. 4.1, where we obtain that the correlator Q 2 becomes constant, as expected. In Sect. 4.2 we have added an external magnetic field B, and studied the diffusion of Q 2 . We obtain that this correlator grows continuously, linearly in time. In both cases, without and with external magnetic field, based on the numerical outcome from the lattice simulations, we have obtained fits to the appropriate coefficients characterizing the behavior of Q 2 in each regime. We then compared these fits to the analytical estimates presented in sections 1,2.
The diffusive behavior of Q 2 in the presence of an external magnetic field B, indicates that the mechanism for fermionic number non-conservation for B = 0 is due to fluctuations of the gauge fields, similarly as in the case of non-Abelian gauge theories. Our numerical determination of such diffusion rate when B = 0, Γ κα 3 B 2 , is perhaps the most important result of the paper. We provide multiple fits to Q 2 in Sect. 4.2, based on different assumptions about the scaling of the correlator with α and B. Within the statistical errors, we obtain that the rate scales, parametrically, with the fine-structure constant as ∝ α 3 , and with magnetic field ∝ B 2 . We thus find that it scales with parameters as it was found long time ago from equations of MHD modified when accounting for anomaly effects. However, the numerical amplitude of the rate we find is much larger than previously estimated, as large as Γ (num) diff /Γ (th) diff 60 when compared to the rate derived from MHD. We argue that this is an indication that the small scale fluctuations of the electromagnetic fields with the size ∼ 1/(αT ), are important and must be accounted for. The higher rate of CS diffusion implies the higher rate of chirality (or fermion number) non-conservation in the presence of a magnetic field. Quite a number of works, e.g. [12][13][14][15][16] and those who used MHD for the chirality non-conservation rate [some of them written by the one of the authors (MS)] used the MHD approach for the study of combined evolution of the chirality/fermion number and the magnetic fields. We believe that the effects found in them may have to be reconsidered in the view of our present findings.
We have also performed numerical simulations in Sect. 5, for a system with a nonzero fermion charge density initially present, characterized by a chemical potential µ with initial value µ o . We have considered again both the absence and presence of an external magnetic field. Our numerical results clearly exhibit the expected instability of the system for B = 0 and non-zero µ o (when the latter exceeds some critical value which depends on the size of the system): the chemical potential behaves as µ(t) = Υ(t) · f (t) with f (t) a quasi-periodic oscillation, and Υ(t) ∼ e −Γt a damping envelope amplitude with decay rate Γ. When B = 0, we observe that the fluctuations of the bosonic fields lead to the transfer of chemical potential into Chern-Simons number for arbitrary µ.
In this paper, we have also provided a number of arguments about the fact that, even in the absence of the magnetic field, the fluctuations of bosonic fields may provide a mechanism for chirality non-conservation, unrelated to instabilities. This mechanism would lead to the damping of the fermion number with a rate independent on the chemical potential. Unfortunately we were not able to test this hypothesis in our lattice simulations due to the lack of computer resources.
There is quite a number of points where our analysis can be refined and improved. First, our study of the classical system with non-zero chemical potential was quite far from the physically interesting domain µ/T 1, and also from continuum limit. The reason for this is that one needs much larger lattices than we were able to used due to inequality (3.24). This also prevented us to verify the relation (2.14) which should be valid for small chemical potentials. Second, we studied only the classical theory without incorporating the hard thermal loops and Bödeker random force [38], the state of art that has been achieved for investigations of the non-Abelian sphaleron rate [41]. Yet another important question is related to applicability of our results for CS diffusion rate to high temperature electrodynamics without any scalar fields. While the Abelian Higgs model serves as a good approximation to the U(1) sector of the SM around the electroweak cross-over, this is not the case for smaller temperatures when all charged bosonic degrees of freedom are decoupled. We leave these points for future work.