Chiral charge dynamics in Abelian gauge theories at finite temperature

TL;DR

None

Contribution

None

Abstract

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

Tables

Table 1. Table 1 : Chiral dissipation rate Γ 5 subscript Γ 5 \Gamma_{5} . The errors are two standard deviations and are not rescaled by the χ 2 / d ​ o ​ f superscript 𝜒 2 𝑑 𝑜 𝑓 \chi^{2}/dof .

Fit quantity	Fit	${\chi }^2/\mathrm{dof}$
${\Gamma }_5$	${10}^{-2.094\pm 0.008}{\left(B\right)}^{2.001\pm 0.009}{\left(e^2\right)}^{2.76\pm 0.01}$	$0.88$
$\frac{{\Gamma }_5}{B^2}$	${10}^{-2.103\pm 0.002}{\left(e^2\right)}^{2.75\pm 0.01}$	$1.13$
$\frac{{\Gamma }_5}{e^{11/2}}$	${10}^{-2.098\pm 0.008}{\left(B\right)}^{2.005\pm 0.008}$	$1.13$
$\frac{{\Gamma }_5}{e^6}$	${10}^{-2.0\pm 0.1}{\left(B\right)}^{2.1\pm 0.1}$	$155$
$\frac{{\Gamma }_5}{e^{11/2}{B}^2}$	${10}^{-2.10\pm 0.02}$	$1.16$
$\frac{{\Gamma }_5}{e^6{B}^2}$	$\left(1.9\pm 0.1\right)\cdot {10}^{-3}\ln \left(\frac{62\pm 13}{e^2}\right)$	$1.80$

Table 2. Table 2 : Diffusion rate Γ diff subscript Γ diff \Gamma_{\rm diff} from Figueroa:2017hun , re-analysed. We take advantage of this re-analysis to take into account the statistical fluctuations in the fit, which was not done in Figueroa:2017hun . The errors represent two standard deviations. The relatively small χ 2 / d ​ o ​ f superscript 𝜒 2 𝑑 𝑜 𝑓 \chi^{2}/dof shows some degree of overfitting. As we do not use it to rescale the errors, they are probably overestimated.

Fit quantity	Fit	${\chi }^2/\mathrm{dof}$
$\frac{⟨Q^2⟩}{Vt}$	${10}^{-3.00\pm 0.23}\cdot {\left(e^2\right)}^{2.87\pm 0.12}\cdot B^{2.06\pm 0.07}\cdot V^{0.05\pm 0.09}$	$0.67$
$\frac{⟨Q^2⟩}{Vt{e}^6}$	${10}^{-3.06\pm 0.24}\cdot {10}^{-3}\cdot B^{2.05\pm 0.08}\cdot V^{0.06\pm 0.09}$	$0.74$
$\frac{⟨Q^2⟩}{Vt{e}^{11/2}}$	${10}^{-2.93\pm 0.24}\cdot {10}^{-3}\cdot B^{2.06\pm 0.08}\cdot V^{0.03\pm 0.06}$	$0.74$
$\frac{⟨Q^2⟩}{Vt{B}^2}$	${10}^{-2.91\pm 0.21}\cdot {10}^{-3}\cdot {\left(e^2\right)}^{2.87\pm 0.12}\cdot V^{0.03\pm 0.04}$	$0.70$
$\frac{⟨Q^2⟩}{Vt{e}^6{B}^2}$	${10}^{-2.90\pm 0.10}\cdot V^{0.06\pm 0.08}$	$0.76$
$\frac{⟨Q^2⟩}{Vt{e}^{11/2}{B}^2}$	${10}^{-2.84\pm 0.20}\cdot V^{0.00\pm 0.04}$	$0.77$

Equations

$$\mathcal{L}=-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\bar{\Psi }{\gamma }^{\mu }{D}_{\mu }\Psi -(D_{\mu }\varphi {)}^{\ast }(D^{\mu }\varphi )-V(\varphi )\text{leavevmode},$$

$$V(\varphi )=m^2|\varphi {|}^2+\lambda |\varphi {|}^4.$$

$${\partial }_{\mu }{J}_5^{\mu }$$

$$n_{{}_{\mathrm{CS}}}=\frac{e^2}{8{\pi }^2}\vec{A}\cdot \vec{B}.$$

$${\partial }_0{J}_5^0$$

$$J_5^0=\frac{1}{6}\mu T^2+\frac{{\mu }^3}{6{\pi }^2}.$$

$$\frac{d}{dt}\left(\frac{1}{6}\mu T^2+\frac{{\mu }^3}{6{\pi }^2}\right)=\frac{e^2}{8{\pi }^2}\frac{1}{V}\int d^{3}x{F}_{\mu \nu }{\tilde{F}}^{\mu \nu }.$$

$$\displaystyle S_{eff}=-\int d^{4}x\Big{(}$$

$$\displaystyle{1\over 2}(\partial_{0}a)^{2}-{1\over 4M^{4}}(\partial_{0}a)^{4}-{e^{2}\over(4\pi)^{2}}{a\over\Lambda}F_{\mu\nu}\tilde{F}^{\mu\nu}\Big{)}\,,$$

$$\displaystyle S_{eff}=-\int d^{4}x\Big{(}(D_{\mu}\phi)^{*}(D^{\mu}\phi)-$$

$$k<\frac{e^2}{4{\pi }^2}\mu =k_c.$$

$${\mu }_c=\frac{8{\pi }^3}{e^{2}N\mathrm{d}x}.$$

$${ϵ}_B(k)=\frac{1}{2{\pi }^2}{\left(\frac{\mathrm{d}x}{N}\right)}^3{⟨k^3|\tilde{B}(k){|}^2⟩}_{|k|},$$

$$⟨Q^2(t)⟩={\Gamma }_{\mathrm{diff}}Vt$$

$${\Gamma }_5=6\frac{{\Gamma }_{\mathrm{diff}}}{T^3}.$$

$$\displaystyle\begin{array}[]{lcl}\frac{\partial\vec{E}}{\partial t}=e\vec{j}-\vec{\nabla}\times\vec{B}-\frac{e^{2}}{4\pi^{2}}\mu\vec{B}\,,&&\frac{\partial\vec{B}}{\partial t}=\vec{\nabla}\times\vec{E}\,,\vspace*{0.3cm}\\ \frac{\mathrm{d}\mu}{\mathrm{d}t}=\frac{3e^{2}}{T^{2}\pi^{2}}\frac{1}{V}\int\mathrm{d}x^{3}\vec{E}\cdot\vec{B}\,,&&\vec{j}=-\sigma\vec{E}\,.\end{array}$$

$$\frac{\mathrm{d}\mu }{\mathrm{d}t}=-\frac{3e^2}{T^2{\pi }^2\sigma }\frac{1}{V}\int \mathrm{d}{x}^3\left(\frac{e^2}{4{\pi }^2}\mu \vec{B}+\vec{\nabla }\times \vec{B}\right)\vec{B}.$$

$$\mu \propto e^{-{\Gamma }_5^{(MHD)}t},$$

$$A_1=f(t)\sin (kz)-\frac{B}{2}y,A_2=f(t)\cos (kz)+\frac{B}{2}x,A_3=f_z(t),$$

$$\ddot{f}$$

$$\dot{\mu }$$

$$\ddot{f_z}$$

$$\ddot{f}=kf\left(\frac{e^2\mu }{4{\pi }^2}-k\right)-\sigma \dot{f},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mu =-k\frac{3e^2}{2T^2{\pi }^2}{f}^2+k\frac{3e^2}{2T^2{\pi }^2}{f}_0^2+{\mu }_0.$$

$$f_0{k}^2-\frac{e^2}{4{\pi }^2}{f}_{0}k{\mu }_0+k^2\delta f(t)+\frac{3e^4}{4{\pi }^4{T}^2}{f}_0^2{k}^2\delta f(t)-\frac{e^2}{4{\pi }^2}k{\mu }_0\delta f(t)+\sigma \dot{\delta f}(t)+\ddot{\delta f}(t)$$

$$+\frac{9f_0{k}^2{e}^4}{8{\pi }^4{T}^2}\delta f(t{)}^2+\frac{3k^2{e}^4}{8{\pi }^4{T}^2}\delta f(t{)}^3=0.$$

$$\left(k^2-\frac{e^2}{4{\pi }^2}k{\mu }_0\right)\left(f_0+\delta f(t)\right)+\sigma \dot{\delta f}(t)+\ddot{\delta f}(t)=0.$$

$$\delta f(t)\approx \frac{e^{-\frac{\sigma }{2}t}}{2{\omega }_{eff}}\left(2f_0{\omega }_{eff}\cosh \left({\omega }_{eff}t\right)+(2\dot{f_0}+f_0\sigma )\sinh \left({\omega }_{eff}t\right)\right)-f_0,$$

$${\omega }_{eff}^2=\frac{e^2}{4{\pi }^2}k{\mu }_0-k^2+{\left(\frac{\sigma }{2}\right)}^2.$$

$$\delta f(t{)}_{as}\approx \frac{\exp \left(\left(-\frac{\sigma }{2}+{\omega }_{eff}\right)t\right)}{2{\omega }_{eff}}\left(f_0\left({\omega }_{eff}+\left(\frac{\sigma }{2}\right)\right)+\dot{f_0}\right).$$

$$r_2(t)=\frac{k^2-\frac{e^2}{4{\pi }^2}k{\mu }_0}{\frac{9f_0{k}^2{e}^4}{8{\pi }^4{T}^2}}\frac{1}{\delta f(t{)}_{as}},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{r}_3(t)=\frac{k^2-\frac{e^2}{4{\pi }^2}k{\mu }_0}{\frac{3k^2{e}^4}{8{\pi }^4{T}^2}}\frac{1}{\delta f(t{)}_{as}^2}$$

$$t_b\equiv \min (t_2,t_3),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathrm{with}\text{leavevmode}\text{leavevmode}{r}_2(t_2)\equiv 1,\text{leavevmode}{r}_3(t_3)\equiv 1.$$

$$\sigma \frac{\partial \vec{B}}{\partial t}={\nabla }^2\vec{B}+\frac{e^2}{4{\pi }^2}\mu (\vec{\nabla }\times \vec{B}).$$

Full text

1 Introduction

Anomalous processes can be relevant in a large number of phenomena, from high energy particle physics to condensed matter. One of the most well-known applications is the explanation of the $\pi^{0}\rightarrow 2\gamma$ decay in quantum electrodynamics (QED), whereas in quantum chromodynamics (QCD) they play a decisive role in the resolution of the $U_{A}(1)$ problem tHooft:1976snw ; tHooft:1976rip . The rate at which anomalous processes occur is actually a relevant quantity whenever we are dealing with out-of-equilibrium processes. In highly energetic and dense matter environments like in the early universe, the fluctuations of gauge and scalar fields – sphalerons Klinkhamer:1984di – lead to rapid fermion number non-conservation in the Standard Model (SM) Kuzmin:1985mm , and to chirality non-conservation in QCD McLerran:1990de . The $SU(2)$ sphaleron rate is a crucial quantity to assess the viability of electroweak baryogenesis Kuzmin:1985mm , and has been extensively studied across decades, both from a purely analytical side and with the help of numerical simulations. Studies have been carried in the pure $SU(2)$ theory in Refs. Philipsen:1995sg ; Ambjorn:1995xm ; Arnold:1995bh ; Arnold:1996dy ; Arnold:1997yb ; Moore:1997sn ; Bodeker:1998hm ; Moore:1998zk ; Moore:1999fs ; Bodeker:1999gx ; Arnold:1999uy and in the electroweak theory e.g. Tang:1996qx ; Ambjorn:1997jz ; Moore:2000mx ; DOnofrio:2012phz , see DOnofrio:2015gop for the latest up-to-date prediction using the measured Higgs-mass.

Anomalous $U(1)$ processes have received some attention, especially in the cosmological context. In the electroweak theory of the SM, the anomaly in the fermionic and/or chiral current actually contains a $U(1)$ contribution, which is associated with the hypercharge field in the Higgs unbroken phase and to the photon field of QED in the Higgs broken phase. As in Abelian gauge theories there are no large gauge transformations, nor vacuum configurations with different Chern-Simons numbers, there is no irreversible fermion (or chiral) number non-conservation, contrary to the case of non-Abelian theories. This does not prevent however the fermion/chiral number in Abelian theories to be transferred into gauge configurations carrying Chern-Simons number, and to re-appear back again due to changes in the gauge field background. These processes may have an important impact on the problems of baryogenesis Giovannini:1997eg ; Kamada:2016eeb ; Kamada:2016cnb ; Fujita:2016igl ; Kamada:2018tcs , magnetic field generation in the early universe Joyce:1997uy , and chiral asymmetry evolution at $\sim\mathrm{MeV}$ temperatures Boyarsky:2011uy . Anomalous $U(1)$ processes have also received a renewed interest in the quark-gluon plasma community, where the chiral magnetic effect Vilenkin:1980fu ; Fukushima:2008xe and its potential experimental signatures are being studied, see Kharzeev:2015znc for a review.

The above $U(1)$ case have been studied mostly within the framework of magneto-hydrodynamics (MHD), which is an effective description accounting for distance scales exceeding the mean free path of the charged particles involved in the problem, see e.g. Brandenburg:2017rcb ; Rogachevskii:2017uyc , or Schober:2017cdw for recent numerical simulations. Despite the relevance of these processes, a full study beyond MHD, taking into account small scale fluctuations in detail, remains to be done. Some attempts in this direction were made in Buividovich:2015jfa ; Buividovich:2016ulp , where out-of-equilibrium techniques were implemented. The main limitations of these studies were the intrinsic numerical cost associated with a full-fledged treatment of fermions. Another approach was initiated in Figueroa:2017hun , of which this paper is a natural continuation. To explain the aim of the present research, let us set up a working model and fix notation. We are interested in the study of physics described by scalar electrodynamics coupled to a massless vector-like fermion field $\Psi$, so that our starting lagrangian is

[TABLE]

where $F_{\mu\nu}$ is the field strength tensor of the $U(1)$ gauge field $A_{\mu}$, $D_{\mu}=\partial_{\mu}-ieA_{\mu}$ and

[TABLE]

Taking a positive squared mass $m^{2}>0$, we can chose its value so that the lagrangian (1) becomes a toy-model for the SM hypercharge sector close to the electroweak phase transition. The chiral fermionic current $J^{\mu}={\bar{\Psi}}\gamma^{\mu}\gamma_{5}\Psi$ is not conserved at the quantum level, and satisfies the anomaly equation

[TABLE]

where $\tilde{F}_{{\mu\nu}}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ is the hodge-dual of $F_{\mu\nu}$, $N_{f}$ is the number of flavours and $K^{\mu}=\frac{e^{2}}{8\pi^{2}}\epsilon^{\mu\nu\rho\sigma}A_{\nu}\partial_{\rho}A_{\sigma}$ is the Chern-Simons current. It follows that the Chern-Simons number $n_{{}_{\rm CS}}=K^{0}$ is identified with the magnetic helicity density [we use a $(-,+,+,+)$ signature and $\epsilon^{0123}=-\epsilon_{0123}=1$]

[TABLE]

In the particular case of a homogeneous fermion distribution, the anomaly equation reduces to (we fix $N_{f}=1$ from now on)

[TABLE]

This allows us to write an effective description without fermions. Integrating them out, they can be represented by a (homogeneous) chemical potential $\mu$ sourced by the Chern-Simons number Redlich:1984md ; Niemi:1985ir . Furthermore, in the case of massless fermions, the relation between the chiral current and the chemical potential can be written in a closed form as111This is obtained for the ensemble average of the number density with the Fermi-Dirac distribution in the presence of a chemical potential: $J_{5}^{0}\equiv n_{+}-n_{-}$, with $n_{\pm}={1\over 2\pi^{2}}\int_{0}^{\infty}dE\,E^{2}\left(1+e^{E\pm\mu\over T}\right)^{-1}$.

[TABLE]

The anomaly, which represents a violation of chirality, can then be recast as an equation for $\mu$

[TABLE]

The equations of motion of the scalar and gauge fields, together with the anomaly equation (7), can be actually derived from an effective action

[TABLE]

upon identifying $\partial_{0}a={\Lambda\mu}$, $\Lambda^{2}=T^{2}/12$ and $M^{4}=\frac{\pi^{2}T^{4}}{9}$, where $T$ is the temperature of the system. The potential $V(\phi)$ was given in Eq. (2).

The main aim of this work is to extract the decay rate of the chemical potential in the presence of an external magnetic field. This quantity is related to Chern-Simons number diffusion rate in the absence of chemical potential through the fluctuation-dissipation theorem (we discuss this in detail in appendix B). It is of special interest as the corresponding diffusion rate has been measured from (independent) simulations in Ref. Figueroa:2017hun . The fluctuation-dissipation relation is a near-equilibrium relation, relevant for small chemical potentials $\mu\ll T$. Having this in mind, we will restrict our attention to the reduced system

[TABLE]

which corresponds to neglecting the $\mu^{3}$ term in the anomaly equation, hence leading to a linear dynamics in $\mu$. We have explored the dynamics of this system for different values of $\mu$, solving the linear equations of motion that follow from Eq. (LABEL:eq:action). In the case of large chemical potential $\mu\gg T$, the correct physical regime is rather described by the non-canonical kinetic term of $a(t)$ that we have neglected in Eq (LABEL:eq:action). Keeping that term as in the full action (8), would lead however to a non-linear set of equations of motion, requiring a more complicated algorithm for solving the evolution of the system in a lattice. In this work we focus for simplicity in the linearised description (at the level of the equations of motion) of the system given by Eq. (LABEL:eq:action), independently whether this corresponds to the correct description of the physical regime expected for a given $\mu$. Fortunately, in order to extract the chiral decay rate of the chemical potential in the presence of a magnetic field, it is enough to explore the regime $\mu<T$, for which Eq. (LABEL:eq:action) represents the correct description of the system.

We start nonetheless studying the system with large non-vanishing initial values of the chemical potential, $\mu>T$. Even though, as mentioned, Eq. (LABEL:eq:action) does not correspond to the correct description of this physical regime, we still expect to get a qualitative understanding of the system dynamics. In this case, the chemical potential is expected to be unstable, forced to decay down to lower values. Due to volume effects unavoidable in any lattice simulation, the chemical potential cannot decay completely, and rather relaxes into a volume-dependent finite critical value $\mu_{c}$. By reaching a large enough volume in the lattice, we have achieved sufficiently small values of $\mu_{c}$, so that we can characterise well the chemical potential decay. In particular, we observe that the decay proceeds through a self-similar/power-law regime, while at the same time the gauge field long-wavelengths develop an inverse cascade spectrum.

Next, we move into the main aim of our present work, studying the system in the presence of a background magnetic field. In this case, the theory becomes similar to its non-Abelian counterpart, as the vacuum becomes degenerate. In particular, $\mu$ decays exponentially all the way down to zero. We have measured the corresponding rate $\Gamma_{5}\equiv-{d\log\mu\over dt}$ in the regime of small $\mu\leq T$, where (LABEL:eq:action) describes the physically relevant theory. We provide in this way a new determination of the parametric dependence of the chiral rate $\Gamma_{5}$, which exhibits a best fit scaling as $\Gamma_{5}\propto e^{11/2}B^{2}$, with no residual volume dependence, where $B$ is the external magnetic strength. Furthermore, we have compared this prediction comparing the magnitude and parametric scaling of $\Gamma_{5}$ obtained from our present lattice simulations, against direct measurements of $\Gamma_{\rm diff}$ from (independent) simulations previously presented in Ref. Figueroa:2017hun , which according to the fluctuation-dissipation argument should be related as $\Gamma_{5}=6\Gamma_{\rm diff}/T^{3}$, with $T$ the temperature of the system.

On the technical side, the discretisation of this theory needs to be done with care, especially when considering the Chern-Simons number. An appropriate discretisation scheme for Abelian gauge theories, reproducing the continuum limit of the theory to quadratic order in the lattice spacing, was presented in Figueroa:2017qmv . It obeys the following properties on the lattice: $i)$ the system is exactly gauge-invariant, and $ii)$ shift symmetry of the axion is exact222An equivalent formulation for non-Abelian gauge theories was originally introduced in Moore:1996qs ; Moore:1996wn . In non-Abelian gauge theories the shift symmetry is however not preserved exactly on the lattice.. Property $i)$ implies that physical constraints such as the Gauss law or Bianchi identities, are exactly verified on the lattice (up to machine precision). Property $ii)$ implies that the lattice formulation naturally admits a construction of the topological number density with a total (lattice) derivative representation $\mathcal{K}\equiv F_{\mu\nu}\tilde{F}^{\mu\nu}=\Delta_{\mu}^{+}K^{\mu}$, which reproduces the continuum expression $\mathcal{K}=\partial_{\mu}K^{\mu}\propto\vec{E}\cdot\vec{B}$ up to $\mathcal{O}(dx_{\mu}^{2})$ corrections. Without this property, the interaction $a\,F_{\mu\nu}\tilde{F}^{\mu\nu}$ cannot be interpreted as a derivative coupling, and hence is not really shift symmetric. As it is precisely the shift symmetry which justifies the functional form of the interaction in first place, it is therefore relevant to preserve exactly such symmetry at the lattice level. Hence, in the present work, we obtain results from numerical simulations based on the discretisation scheme presented in Ref. Figueroa:2017qmv . In appendix, we provide a summary of the key equations of such lattice formulation. For further technical details we refer the interested reader directly to Ref. Figueroa:2017qmv .

Most of what is known about anomalous $U(1)$ dynamics comes from MHD predictions, which represent the long-wavelength effective description of our system. To get analytic insights, it is interesting to compare our findings with this effective theory. The phenomena we observed (e.g. self-similar decay of the chemical potential, gauge field inverse cascade dynamics, etc) are qualitatively well-modeled by MHD-inspired models, even though these do not take into account the full-fledged short-scale dynamics present in our lattice simulations.

The paper is organised as follows. In section 2, we present the outcome of our lattice simulations. Namely, we study the evolution of an initially non-vanishing large chiral chemical potential $\mu$ in the absence of an external magnetic field. This leads to a characterisation of the self-similar decay of $\mu$, and of the gauge field inverse cascade dynamics. We then add an external magnetic field, which introduce a vacuum degeneracy, and allows for the complete decay of the chemical potential. This gives us a way to measure in the physically relevant regime the amplitude and parametric dependence of the decay rate $\Gamma_{5}$, which we compare with $\Gamma_{\rm diff}$ as inferred from previous simulations. Section 3 is devoted to the comparison of our results to MHD. We provide there a simplified MHD model which allows us to describe the gross feature of the system. We then characterise the inverse cascade phenomenon and see that its dynamical evolution can be well fit to an MHD-like $ansatz$. We also discuss the expected parametric dependence of the chiral rate. In section 4 we discuss our results and present future outlooks. Appendix A recaps the lattice setup and in appendix B we review the derivation of the fluctuation-dissipation relation between the chiral rate and the Chern-Simons diffusion one.

2 Lattice results

In this section, we report the results obtained from our lattice simulations. We select from a thermal ensemble (with $\mu=0$) initial configurations, which we then evolve following the (lattice version of) the classical equations of motions derived from (LABEL:eq:action), with initial condition $\mu(0)=\mu_{0}\neq 0$. We follow the evolution of the chemical potential, which in the symmetric phase is expected to be unconditionally unstable, leading to the creation of long-range gauge fields Rubakov:1985nk ; Rubakov:1986am . We start by investigating the general features of the chemical potential evolution in the absence of a background magnetic field, making sure that the volume dependence of our observables is under control. Then, we study the self-similar behaviour of the decay for large initial chemical potential, and characterise the associated inverse cascade dynamics of the gauge field. We show that the magnetic power spectrum flows to the infrared modes. After that, we switch on a homogeneous background magnetic field. This is achieved by the use of twisted boundary conditions Kajantie:1998rz , which are already introduced in the Monte-Carlo process that generates the initial configuration (we describe this in appendix A). We observe different regimes, and in particular, we focus on the exponential decay induced by the presence of the external magnetic field. This allows us to extract the chiral magnetic decay rate $\Gamma_{5}=-{d\log\mu\over dt}$, and characterise its parametric dependence.

2.1 Chemical potential decay

The existence of an instability can be understood by looking at the free energy of the system under consideration. In momentum space, the usual magnetic term ${\mathcal{H}}_{\rm mag}={1\over 2}{\vec{B}}^{2}\sim{1\over 2}k^{2}A^{2}$ competes with the chemical potential term ${\mathcal{H}}_{{}_{\rm CS}}=\mu n_{{}_{\rm CS}}\sim\pm\frac{e^{2}}{8\pi^{2}}\mu kA^{2}$, where the different sign corresponds to the different polarisations of the gauge field, reflecting the chiral nature of the coupling. The chiral term may become dominant at small momenta, and for the gauge field polarisation for which ${\mathcal{H}}_{{}_{\rm CS}}$ and ${\mathcal{H}}_{\rm mag}$ have opposite sign, an instability occurs for sufficiently infrared modes

[TABLE]

In particular, if all $k$’s on the lattice are larger than the critical $k$, no instability can develop. As lattice momentum is discrete and has a minimum value $k_{min}$, it implies the existence of a critical chemical potential below which no instability can develop. Equation (10) can be rewritten to understand what is the largest chemical potential which is stable given a momentum resolution. Using the smallest momentum in a lattice $k_{min}=\frac{2\pi}{N\mathrm{d}x}$, with $N$ the number of lattice points in one direction and $\mathrm{d}x$ the lattice spacing, we find

[TABLE]

For $\mu<\mu_{c}(e,N,\mathrm{d}x)$ the instability is not captured any more for the given lattice. From now on we set $T\mathrm{d}x=1$, as any re-scaling of the lattice spacing can be translated into an inverse re-scaling of the gauge coupling constant $e^{2}$. In this way, relevant scales are always guaranteed to be captured in the lattice (see Figueroa:2017hun for more details). In the upper panels of figure 1 we show the behaviour of the chemical potential for different lattice sizes (in this figure and in the following, we indicate in brackets the corresponding units in powers of the temperature $T$). The critical threshold predicted by (11) is clearly appreciated, as indicated by the dotted lines in the figure. As long as the chemical potential remains much larger than the critical value for a given lattice size $N$, all simulations agree, as expected, independently of the volume. On the other hand, some finite volume effects are observed around the critical value. In particular, we see in the top right panel that for the smallest volume considered $N=56$, the chemical potential displays some oscillations just before approaching the critical value; such oscillations disappear for larger volumes. In the lower panel, we display the chemical potential evolution for different initial values, for the largest $N=448$ lattices we simulated. For large enough initial values, we observe a power law decay, which as we will explain in section 2.2, corresponds to a self-similar behaviour.

Another phenomenon we observe is that the smaller the initial chemical potential, the longer it takes for the decay to be triggered. In other words, we see the emergence of long $plateau$s, gradually longer the smaller the initial value of the chemical potential is. This phenomenon makes the regime of very small chemical potential difficult to be captured well on a lattice. Atop of requiring large volumes to decrease the critical chemical potential, one has to perform increasingly longer simulations.

2.2 Inverse cascade

In the previous subsection, we have seen that finite volume effects are under control for sufficiently large lattices $N\gg 50$, and that long initial $plateau$s appear for small chemical potentials. A remarkable feature of the chemical potential evolution in this context is that for large initial values, its decay corresponds to a universal power-law. We show this phenomenon in figure 2 for different charges and initial values of the chemical potential. After some transient behaviour, all decays enter a power-law regime well described by a $t^{-1/2}$ behaviour. All the dependence on the parameters comes from the charge and lies in the prefactor as we illustrate in the right-hand side part of the figure.

Actually, this phenomenon goes in pair with what is known as an inverse magnetic cascade. The decay of the chemical potential induces a transfer of magnetic energy from the ultraviolet (UV) into the infrared (IR). In other words, long-range magnetic fields are created. We consider the magnetic power-spectrum

[TABLE]

where the quantity $\langle\dots\rangle_{|k|}$ denotes an angular average over the spherical shell of radii $[|k|-{1\over 2}\Delta k,|k|+{1\over 2}\Delta k)$, with $\Delta k=k_{\rm min}\equiv{2\pi\over Ndx}$ the binning width. The quantity $\tilde{B}(k)=\sum_{x}e^{\vec{-ik}\cdot\vec{x}}\vec{B}(x)$ is the lattice discrete Fourier transform of $\vec{B}(x)$ and has the same units as $\vec{B}(x)$. This definition has the advantage of being volume independent and it is simply related to the real-space volume average as $\langle{\vec{B}}^{2}\rangle_{{}_{V}}=\sum_{k}\frac{k_{\rm min}}{k}\,\epsilon_{B}(k)$ mimicking the continuum relation $\langle{\vec{B}}^{2}\rangle=\int\frac{\mathrm{d}k}{k}\,\epsilon_{B}(k)$.

In figure 3, we plot the time evolution of the magnetic spectrum, where an energy flow from UV scales into IR scales, i.e. an inverse cascade, is clearly observed. As the chemical potential starts decaying, the magnetic energy is gradually transferred into the lower modes, so that at the end of the simulation most of it is peaked around the smallest lattice mode $k_{\rm min}$. Another way to display this information is by looking directly at the spatial distribution of the magnetic field, as we do in figure 4. There we show snapshots of the magnetic field arrows, together with some representative field lines emanating from the center. In the left panel of figure 4, we see that, immediately after thermalisation, no characteristic structures in the magnetic field distribution are observed. However, after the inverse cascade process occurs, we see long-range magnetic fields, as depicted in the right panel of figure 4, where the field lines reach out through the whole lattice.

2.3 Decay in the presence of an external magnetic field

The dynamics of a chiral charge is of special interest in the presence of a background magnetic field $B_{i}=\frac{1}{2}\epsilon_{ijk}F_{jk}$. In this case, as detailed in Figueroa:2017hun , the vacuum of the theory is degenerate and the situation is closely related to its non-abelian counterpart. As already mentioned, such a background can be introduced on the lattice through the use of twisted boundary conditions Kajantie:1998rz . The main difference from the case without magnetic field is that the chiral charge is now unconditionally unstable, and can decay all the way down to zero.

The early dynamics of the system depends crucially on the strength of the magnetic field. For large enough ones, it will drive the chiral charge decay from the beginning onward. For weak external magnetic fields, we expect the system to evolve initially like in the absence of magnetic field, and settle down to the critical value $\mu_{c}$. The chemical potential cannot stay however in the state $\mu=\mu_{c}$ in the presence of a background magnetic field, and eventually decays to zero.

These two limiting situations can be well observed in our simulations. We can also study the transition between them, see in particular figure 5. There we plot the contribution to the total energy from the chemical potential $\frac{\mu^{2}T^{2}}{24}$333Of course this represents a ”physical” energy only for $\mu<T$. . This allows us to compare the chemical potential energy with respect to that of the external magnetic field. In the upper panel, we plot configurations with large initial values of the chemical potential, whereas the analogous plots for weaker initial chemical potential values are presented in the lower panel. As a reference, we also display the outcome from simulations with the same initial chemical potential but no magnetic field. The dashed lines represent the magnetic energy carried by the external background. Looking first at the weakest magnetic fields, we see that its influence on the chemical potential is minimal. For larger initial chemical potentials, the simulations are almost not affected while in the case of smaller chemical potentials the length of the initial $plateau$ is slightly reduced. As expected, we see that the chemical potential relaxes first to its critical value, before undergoing a secondary decay due to the presence of the background magnetic field. For stronger background magnetic fields (top red curve, and bottom red and blue curves), we observe a transition regime where both effects are competing. In this case, the system is not sensitive any more to the finite number of infrared modes; the chemical potential does not stabilise to a critical value. In both cases, it happens when the external magnetic energy is roughly two orders of magnitude smaller than the initial chemical potential energy. For even larger magnetic energies, we enter into a different dynamical regime, where the chemical potential decays quickly to zero through damped oscillations.

It is also instructive to look at the spatial distribution of the magnetic field. In figure 6, we show the magnetic field lines corresponding to the orange and red curves in the top panel of Fig. 5, associated with large initial chemical potentials. In the left panels of Fig. 6, we show the field lines just after thermalisation, and in the right panels we show the field line distribution at the final stage of the decay. As expected, the magnetic field is initially oriented along the direction of the external magnetic field (which, without loss of generatily, was chosen to be oriented along the $z$-axis), and the stronger the external magnetic field the less important the thermal fluctuations are. In the final stage, we still see lattice-size structures. For weaker magnetic fields, the decay is essentially driven by the chiral instability of the chemical potential, so the field lines extend in all directions, as in the case without magnetic fields. Yet the effect it is still visible due to the presence of the external magnetic field, as more field lines expand along the $z$-direction. For a larger magnetic field background, the field lines which are not aligned to the $z$-axis are suppressed at the end of decay, and we essentially see only field strength lines along the $z$-axis.

2.4 Chiral magnetic rate

As presented in the previous section, the evolution of the system in the presence of an external magnetic field is subject to two competing dynamics. On the one hand, there is the intricate evolution related to the anomaly-induced creation of long-range gauge fields. On the other hand, the external magnetic field acts as a vacuum reservoir and induce an exponential decay of the chemical potential. In this section, we are interested in this dynamic.

To extract $\Gamma_{5}$ from our simulations, we proceed as follow. To disentangle the magnetic driven decay from the rest of the dynamics, we take our initial chemical potential to be the critical one. Thanks to our large volumes, it also allow us to reach the physically relevant regime $\mu<T$. As explained in section 2.1, the anomaly-related instability does not develop for sub-critical chemical potentials. Representative evolution of the critical chemical potential in the presence of an external field are shown in the upper-left plot of 7. After an initial relaxation, we see, as expected, a clear exponential decay. This allows us to measure $\Gamma_{5}$. We show its $B^{2}$ dependence for different volumes in the upper-right panel of 7. We observe a good convergence to the infinite volume limit, as the $N=448$ and $N=224$ cases differ only in $\sim 0.5\%$. We thus pursue the investigation only on the $N=448$ lattices, adding a $1\%$ systematic error to take into account any potential remaining volume dependence. We obtain that there is virtually no deviation from the expected $B^{2}$ dependence.

In the lower-left panel of Fig. 7, we look at the charge dependence of the rate. Fitting our data, we find that they are well modeled by a $\propto e^{11/2}$ dependence. This can also be seen in the lower-right panel of Fig. 7, where we show $\frac{\Gamma_{5}}{e^{11/2}B_{ext}^{2}}$. All the different determinations collapse to a constant; there is very little remaining dependence on the charge and magnetic field.

More quantitatively, we present in Table 1 the different fits we performed to our data. All the main fits agree with each other, giving an almost exact $B^{2}$ magnetic dependence, an effective $e^{11/2}$ charge dependence and a prefactor of $10^{-2.1}\approx 0.0079$. To check the robustness of the coefficients, we also to tried to fit $\frac{\Gamma_{5}}{e^{6}}$ with $f(B)=\alpha B^{\beta}$. Such a hypothesis is clearly excluded by the data as it leads to a $\chi^{2}/\mathrm{dof}$ of $155$. We also present in this table a fit in $\Gamma_{5}=aB^{2}e^{6}\ln\left(\frac{b}{e}\right)$ to show that the charge dependence can also potentially interpreted as logarithmic corrections to a leading $e^{6}$ scaling. As we will discuss in section 3.3, this is what can be expected from theoretical predictions.

These numbers may be compared to the values of the topological charge diffusion rate, which at large time reads

[TABLE]

with $Q(t)=N_{cs}(t)-N_{cs}(0)$, $N_{\rm cs}(t)=\frac{e^{2}}{8\pi^{2}}\int\mathrm{d}x^{3}\vec{A}\cdot\vec{B}$. This quantity was measured in Figueroa:2017hun , and it is related to $\Gamma_{5}$ through a standard fluctuation-dissipation argument444Note that the corresponding expression in Figueroa:2017hun contains a factor $2$ wrong, see appendix B here, for a clarification of this.

[TABLE]

For the sake of the comparison, we present in Table 2 a re-analysis of the data of Ref. Figueroa:2017hun . Unfortunately, the diffusion rate is a harder quantity to extract and leads to a less precise quantification. As the results of  Figueroa:2017hun were obtained at different small volumes, we keep an explicit volume dependence in our fits, even if we find it to be weak. Taking as a measured value $\Gamma_{5}=10^{-2.10\pm 0.02}B^{2}e^{11/2}/T^{3}$ leads to a prediction $\Gamma_{\rm diff}=10^{-2.88\pm 0.02}B^{2}e^{11/2}$. The prefactor and the magnetic field dependence are consistent with the measurements of $\Gamma_{\rm diff}$ reported in Figueroa:2017hun . At first sight, the scaling with $e^{2}$ seems to be in mild tension, but further analysis shows that the data from Ref. Figueroa:2017hun does simply not constrain well enough the charge dependence of $\Gamma_{\rm diff}$. This can be seen in the fits in Table 2, where either scalings $\propto e^{11/2}$ and $\propto e^{6}$ are compatible with the data, exhibiting a similar $\chi^{2}/dof$. We conclude that within the errors, the fluctuation-dissipation relation Eq. (14) is well verified.

3 Comparison to MHD

Different approaches can be followed to get some analytical understanding of a system in the presence of a chiral chemical potential. Recently, progresses have been made towards establishing a $chiral-$kinetic theory, which adds to the usual Boltzmann-like approach a self-consitent treatment of the anomaly, see Chen:2015gta and references therein for more information. Valuable knowledge may also be gained in the hydrodynamical regime by studying the effective field theory of the long wavelength field modes together with the electromagnetic interactions. In this approach, matter is a fluid, characterised by its velocity field. A variety of interaction terms going from the simple linear response electric conductivity to direct coupling of the velocity which can induce turbulence can be added. This kind of approach is referred to as MHD. The equations can be further modified to take into account anomalous processes Giovannini:1997eg ; Joyce:1997uy , to give anomalous MHD. Recently, an extensive study was carried out, both at the theoretical level Rogachevskii:2017uyc and using numerical simulations, Schober:2017cdw of full anomalous MHD models. We refer the interested reader to these references for more details.

To gain some understanding of our results,we will attempt to model the system using a truncated set of MHD equations. Namely, we will focus on the magneto-dynamics part, not considering any hydrodynamics evolution.

[TABLE]

where $\sigma$ is the linear-response electric conductivity.

This approximation will turn out to qualitatively work surprisingly well. Indeed, as our lattice simulations describes the full-dynamics of our matter field, there is no $apriori$ reasons for the matter field to contribute little to the global dynamics. It is even likely that in some regimes not explored in this work, matter effects which would be described as turbulence in a hydrodynamic language should manifest themselves.

In particular, this model predicts the existence of a chiral magnetic rate. Indeed, the evolution equation for $\mu$ can be rewritten as, neglecting the time derivative of the electric field,

[TABLE]

Taking $B$ to be constant, we then obtain

[TABLE]

with $\Gamma_{5}^{(MHD)}=\frac{3e^{4}}{4\pi^{4}\sigma T^{2}}$. We will compare it to our results in section 3.3.

Finally, note that from now on, when we refer to MHD or MHD-like predictions, we have in mind the specific model (17).

3.1 Qualitative behaviour and initial plateaus

As presented in previous sections, the chemical potential is transferred into long-range helical magnetic fields which carry a non-vanishing Chern-Simons number. To understand at a qualitative level some of the observed features, like the existence of $plateau$s in the initial stage, we can consider the evolution of a maximally helical field under anomalous MHD-like dynamics. In the presence of a background magnetic field in the $z$- direction, this corresponds making the $ansatz$ Rubakov:1985nk

[TABLE]

which leads to

[TABLE]

Let us first consider this toy-model without an external magnetic field ($B=0,f_{z}(t)=0$). Then the system of equations reduces to

[TABLE]

By plugging in the expression for $\mu$ into the first equation of (24), we obtain a non-linear equation for $f(t)$. Solving it numerically also gives us access to the evolution of $\mu(t)$. Here, we take $f_{0}$ and $\sigma$ to be some free parameters to be fitted. Examples of the time evolution of $f(t)$ and $\mu(t)$ obtained with this procedure, are shown in figure 8. Despite being in qualitative agreement, this simple modeling does not capture completely the fine details, especially for large initial chemical potentials. In figure 9, we compare the modeling based on the ansatz (20) with the numerical outcome from our lattice simulations. The lower panel shows the chemical potential whilst the upper panel shows the Chern-Simons number, which is useful to understand the initial dynamics. Solid lines were obtained from our simulations, while dashed and dotted ones are numerical solutions to equations (24) for different parameters (conductivity and initial conditions).

Let us first consider the case $\mu_{0}=8\mu_{c}$, as the model gives a better description for relatively small chemical potentials. The dotted curves were obtained by fitting the $ansatz$ model to the whole range of data. When looking at the chemical potential, we see that the model is able to describe reasonably well the initial $plateau$, but fails to describe well the decay. Actually, when looking at the Chern-Simons number, we see that even the initial phase is not very well described by this fit. A remedy to this is to restrict the fitting range to early times, where the system should be closer to the model. For example, the dashed curves are obtained by restricting the fitting range to times $t<10/T$ (delimited by a solid vertical black line). This gives a much better description of the initial phase; essentially the initial evolution is well captured by the simple $ansatz$. At a later stage, when a more complex dynamics has developed, such as the creation of long-range magnetic fields following a process of inverse cascade (or what can be described as a turbulent regime in MHD), the model is too simplistic to capture well the physics. The same can be said about large initial chemical potentials, as in the case of $\mu_{0}=60\mu_{c}$ depicted in the figure. There, we see explicitly that the power-law decay (expected due to the generation of long-range magnetic field) is completely missed by the $ansatz$ modeling. The red and orange vertical solid lines are a visual guide indicating the saturation of the Chern-Simons number and of the corresponding chemical potential.

We can get a better understanding of the early dynamics of $f(t)$ and $\mu(t)$ in the $ansatz$ modeling, by simply perturbing equation (24), i.e. considering $f(t)=f_{0}+\delta f(t)$, so that

[TABLE]

Keeping the linear order terms in $\delta f$ and neglecting the terms proportional to $f_{0}^{2}$ (we do not start in a state with an helical field), the equation reduces to

[TABLE]

This is a driven harmonic oscillator which admits a solution as

[TABLE]

where $\dot{f_{0}}$ is the initial time derivative of $f(t)$, and we have defined

[TABLE]

Eq. (27) indicates, first of all, that the solution corresponds to an IR instability for the modes $k<\frac{e^{2}}{4\pi^{2}}\mu_{0}$, which grow exponentially fast. Secondly, by estimating the range of validity of this solution, we can estimate the duration of the chemical potential $plateau$ until the onset of its decay. Indeed, the breakdown of this approximation corresponds to the end of the exponential growth of $f(t)$, which in turn triggers the chemical potential decay. The approximations cease to be valid when higher order perturbations become non-negligible. To see this in detail, we keep only the exponentially growing part of the solution and neglect the constant term,

[TABLE]

The relative weight between the linear to second-order perturbation terms in the equation of motion 25, and linear to third order terms, are defined by the ratios

[TABLE]

We expect the linear approximation to breakdown whenever either of these ratios become of order one. We thus define the breakdown time $t_{b}$ to be

[TABLE]

In figure 8, we show how this prediction performs, comparing the linearisation (dotted lines) to the numerical solution to Eq. (24). We see that $t_{b}$ in (31) gives a reliable prediction of the range of validity of the linearisation regime. In conclusion, the duration of the initial $plateau$ can be estimated well with our simple MHD-inspired $ansatz$ equation (20).

In the presence of an external magnetic field, the system of equations (21) does not admit a simple enough analytic treatment. We can solve them nonetheless numerically. We show the chemical potential obtained from numerical integration of these equations in figure 10. As in the previous case, we see that it captures reasonably well the dynamics of small chemical potentials. For weak magnetic fields, it is dominated by the chiral instability of the chemical potential. The system evolves as in the absence of magnetic fields, reaches a $plateau$ and only at later times it decays to zero because of the presence of the magnetic background field. For a sufficiently large background magnetic field, the decay is driven by the presence of such field, driving quickly the system into a $\mu=0$ state. For intermediate external magnetic fields, we have a superposition of the two effects.

3.2 Self-similarity and inverse cascade

The power law decay $\propto t^{-1/2}$ of the chemical potential reported in section 2.2 was also predicted in the context of MHD in several works Joyce:1997uy ; Boyarsky:2011uy ; Boyarsky:2015faa ; Tashiro:2012mf ; Hirono:2015rla . More specifically, various analysis of the time evolution of the magnetic spectrum, have been carried out in refs. Joyce:1997uy ; Tashiro:2012mf ; Hirono:2015rla . We will now compare these modelings with the outcome from our numerical simulations. First of all, we repeat the analysis of refs. Joyce:1997uy ; Tashiro:2012mf ; Hirono:2015rla . The starting point in MHD are equations  (17). Taking the curl of the first equation and neglecting the second time derivative of $\frac{\partial^{2}\vec{B}}{\partial t^{2}}\approx 0$, the equation for the magnetic field can be recast as

[TABLE]

We can now decompose the magnetic field in an orthonormal helicity basis $Q^{\pm}(\vec{k},\vec{x})$, which corresponds to a basis555The precise form of the basis depends on the base space the analysis is performed, e.g. Ref. Joyce:1997uy ; Tashiro:2012mf considers $\mathbb{R}^{3}$ whilst Ref. Hirono:2015rla used a sphere. In the analytical description above we will simply work in a Euclidean base space $\mathbb{R}^{3}$, as we are describing the continuum theory. Were we to write down a basis for the torus, we could construct it nonetheless from a superposition of the basis of Tashiro:2012mf . of eigenvectors of the curl operator,

[TABLE]

so that $\vec{k}\times Q_{\pm}(\vec{k})=\pm kQ_{\pm}(\vec{k})$. In that basis Eq. (32) reduces to

[TABLE]

which admits as a solution

[TABLE]

Furthermore, as soon as $\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\mu(t^{\prime})$ becomes sizable, $\alpha^{-}(t,\vec{k})$ becomes sub-dominant. Neglecting the latter, we can then write

[TABLE]

Equation (36) represents a prediction based on MHD, which we can use to compare against our data. A further simplification we make, following Tashiro:2012mf , is to neglect the $\vec{k}$ dependence666Note that we could also get rid of it by considering ratios of $|B(t,\vec{k})|^{2}$ at different times. An analysis using this method turned out to be less accurate, we believe because of the small number of momenta we can use for the fits. of the initial amplitude $\alpha^{+}(\vec{k})\approx\alpha_{0}^{+}$. This leads us to write down the following function

[TABLE]

which will be used as a fit to $|B(k,t)|^{2}$ obtained from our simulated data at different times. This is precisely what we show in figure 11. In the left panel we plot $|B(k,t)|^{2}$ for three different times together with a best fit of the form (37). Despite the limited number of modes available, we see that the data is well fit by a functional form as in Eq. (37).

Let us look now at the time dependence of the fit coefficients $c_{1},c_{2}$ and $c_{3}$. From comparing (37) and (36), we expect $c_{1}$ to be a linear function of $t$. This is exactly what we observe in figure 11, where the top right panel displays a linear fit from our data as $c_{1}(t)=(18.0\pm 0.1)t+(256\pm 20)$. The crossing of the fit with the y-axis encodes actually various effects. First of all, we learn that $t_{0}$ is not zero, and rather corresponds to an initial time at which $|B(t,\vec{k})|^{2}$ is equal to $|\alpha^{+}_{0}|^{2}$ (we will estimate it shortly, in a self-consistent manner, from the fit to the second coefficient $c_{2}$). Secondly, a non-zero crossing may also reintroduce some $k^{2}$ dependence in the initial condition, which will partially compensate for the approximation  $\alpha^{+}(\vec{k})\approx\alpha_{0}^{+}$.

The coefficient $c_{2}$ is expected to go as $M(t_{0},t)\equiv\frac{1}{t}\frac{e^{2}}{4\pi^{2}}\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\mu(t^{\prime})$. Its dependence on $t_{0}$ allows us to estimate this initial time, by comparing it to $M(t_{0},t)$ for different values of $t_{0}$, and looking for the best match. This is shown in the middle right panel of figure 11. This procedure leads to an estimated value of as $t_{0}\approx 9.2/T$. A last check is provided by the lower right panel of figure 11, which shows how close $c_{3}$ is to a constant. The fact that this MHD-inspired model achieve such a good description of the data is a bit surprising, as the chemical potential probed here lie outside of the range of validity of (17). The relevant length scales which are to be compared are the typical instability length scale associated to the chemical potential decay, $l_{inst}\sim\frac{1}{k_{inst}}=\frac{4\pi^{2}}{e^{2}\mu}$ to the typical mean free path in the plasma $\lambda\sim\frac{1}{\alpha^{2}T}$. MHD is expected to be valid for $l_{inst}>\lambda$, which happens for $\mu<e^{2}T$. The only place where the ”breaking” of the effective theory can be seen is if we try to extract a conductivity from the coefficient $c_{1}$. This gives an effective conductivity $\sigma_{eff}=(0.11\pm 0.01)T$, which differs from the kinetic theory prediction, see section 3.3 for more details.

An interesting question would be to compare $\sigma_{eff}$ to the intrinsic measurement of our system’s conductivity with the use of the Kubo formula Kubo:1957mj . For classical field theories on a lattice with finite volume $V$, for a homogeneous and isotropic plasma, $\sigma$ can be written as Figueroa:2017hun

[TABLE]

where $j_{i}$ are the spatial components of the electric current, and $\langle...\rangle$ is an ensemble average over different realisations of our thermal initial conditions. This formula is readily to be implemented as an observable in our lattice simulations. The expression is however UV sensitive (in the free theory it is even linearly divergent with momentum), so the results are expected to be sensitive to the natural UV cutoff imposed by the lattice $k_{\max}\sim 1/dx$. The correlator in (38) is expected to decay exponentially in time, $\Sigma(t)\sim\exp{(-\gamma t)}$, and to oscillate with the plasma frequency. When obtaining numerically $\Sigma(t)$ from our simulations, it exhibits an oscillatory pattern with frequency considerably smaller than $1/T$, which is presumably attributed to the plasma frequency (plots of this can be found in figure 7 (right panel) of Ref. Figueroa:2017hun , so we do not reproduce them here again). After some time, oscillations occur around zero, indicating the dumping. However, the presence of short time oscillations, associated with the lattice UV cutoff, is also very noticeable. These contribute to make the behaviour of $\Sigma(t)$ very ’noisy’, making difficult to obtain a trustworthy attempt to extract $\sigma$ by this procedure.

Note that even if we were to extract a conductivity from this method, there is no reason to expect it to be equal to the quantum field theory computation. Indeed, as it is a UV-dominated quantity, it is sensitive to the distribution of large-momentum modes, which are intrinsically different between the classical and quantum field theory, see Arnold:1997yb for a careful discussion on this issue.

3.3 Chiral magnetic rate and MHD

As mentioned at the beginning of the section, one can infer a chiral magnetic-rate from MHD

[TABLE]

To be able to use this prediction, we need to use the kinetic theory prediction for the conductivity. The computation consists in finding the coefficient $\sigma$ in the relation $\vec{j}=-\sigma\vec{E}$ where $\vec{E}$ is an external uniform and time-independent electric field and $\vec{j}$ is the induced electric current. The main process to be taken into account is the mutual scattering of charged particles with an energy of the order of temperature. At leading-log order, for one charged fermion species, it can be estimated as Arnold:2000dr

[TABLE]

For one charged scalar species, the conductivity has been also obtained recently as Sobol:2019kiq

[TABLE]

The coefficients $q_{\rm F},q_{\rm S}$ are $\sim\mathcal{O}(1)$, but can only be fixed from matching a computation at full leading order. The leading behaviour of the chiral rate is then expected to scale, independently of whether we use the fermion or scalar computation, as

[TABLE]

with $q$ some $\mathcal{O}(1)$ parameter.

In the case of a charged fermion species, the coefficient $q_{\rm F}$ in the log can be determined using the full leading numerical result, for instance evaluated at $e^{2}=\frac{4\pi}{137.04}$, as can be read out from table 2 of Ref. Arnold:2003zc . Matching such value to Eq. (40), leads to $q_{\rm F}=4.2$, and hence to

[TABLE]

Unfortunately, there is no analogous full leading computation available for the scalar field case, so we cannot determine exactly the value of $q_{\rm S}$. Assuming that $q_{\rm S}=q_{\rm F}$, then $\sigma_{\rm S}$ remains $\sim 14.5\%$ smaller than $\sigma_{\rm F}$. If we consider, however, values within the range e.g. $q_{\rm F}-2\leq q_{\rm S}\leq q_{\rm F}+2$, then $\sigma_{\rm S}$ ranges, for $e=1$, between $\sim 56\%$ larger and $\sim 33\%$ smaller, compared to $\sigma_{\rm F}$ given by Eq (43) evaluated at $e=1$. The difference becomes however smaller for smaller values of $e$.

In section 2.4, were we extracted this rate in the physical region $\mu<T$ , we observed an almost perfect $B^{2}$ dependence of $\Gamma_{5}$, but a certain deviation from an exact $\propto e^{6}$ scaling. In particular, fitting to a simple power law $\propto e^{p}$, we obtain a best fit with $p=11/2$, whereas fitting to a form $e^{6}\ln\frac{q}{e}$ we obtain $q=7.9$, see Table 1777Note that we present there the result of the logarithmic fit as a function of $q^{2}\approx 62$.. Both of these fits describe well our data, although the logarithmic fit leads to a slightly worse $\chi^{2}/dof$. Given our limited range of values of $e^{2}$, we cannot clearly discriminate between these two options.

Using as a reference $\sigma_{F}$, we can now estimate a chiral magnetic decay rate from MHD and compare it to our simulations. Plugging in numbers, we get

[TABLE]

which for $e^{2}=1$ gives $\sigma_{\rm F}\simeq 10.95\,T$. In the presence of a background magnetic field $B$, the chiral rate is actually inversely proportional to the conductivity, and it is given by

[TABLE]

which for $e^{2}=1$ is $\Gamma_{5}^{(MHD)}(e^{2}=1)\simeq 7.0\cdot 10^{-4}\cdot B^{2}/T^{3}$. In general then, $\Gamma_{5}^{(MHD)}\sim 7\times 10^{-4}\cdot e^{6}\cdot B^{2}$, modulo electric charge logarithmic corrections. Alternatively, from the fluctuation-dissipation relation $\Gamma_{\rm diff}={1\over 6}\Gamma_{5}T^{3}$, the effective diffusion rate is expected in MHD as

[TABLE]

which e.g. for $e^{2}=1$ is $\Gamma_{\rm diff}^{(MHD)}\simeq 1.1\cdot 10^{-4}\cdot B^{2}$, or again in general $\Gamma_{\rm diff}^{(MHD)}\sim 10^{-4}\cdot e^{6}B^{2}$ modulo logarithmic corrections on the electric charge.

Comparing these rates against our numerical fits to $\Gamma_{5}$ (c.f. Table 1) or against the fits to $\Gamma_{\rm diff}$ from Ref. Figueroa:2017hun (c.f. Table 2), we observe that the numerically extracted rates are a factor $\sim 10$ larger than the MHD counterparts. In particular for $e^{2}=1$ and imposing an exact scaling $\propto B^{2}$ over our data (something that it is very well verified, recall the discussion in section 2.4), we obtain888The ratio ${\Gamma_{5}^{(num)}/\Gamma_{5}^{(MHD)}}$ of Eq. (47) was obtained with the power law fit $\Gamma_{5}^{(num)}\propto e^{11/2}$. If the ratio is obtained instead with the logarithmic fit $\Gamma_{5}^{(num)}\propto e^{6}\ln\left(\frac{62}{e^{2}}\right)$, this yields (for $e^{2}=1$) ${\Gamma_{5}^{(num)}/\Gamma_{5}^{(MHD)}}=11.3\begin{subarray}{c}+1.0\\ -1.4\end{subarray}$, which is equally consistent with the other ratios.

[TABLE]

The two ratios in Eq. (47) are, as expected, consistent with each other, even if these are numbers independently obtained. This is simply due to the fact that the fluctuation-dissipation argument relating $\Gamma_{5}$ and $\Gamma_{\rm diff}$, is actually well verified in our data (within the errors of the diffusion rate), recall section 2.4. In summary, both the chiral decay rate or the Chern-Simons diffusion rate, obtained by completely independent simulations, lead to rates much larger than predicted in the MHD picture. This is the main quantitative result of our work.999Note that this result can also be interpreted as having an effective conductivity $10$ times smaller than that coming from kinetic theory.

The equations 47 lead us to conclude that formula (39) is not describing correctly the rate of fermion number non-conservation in the presence of magnetic fields in our system. We take this as an evidence of the impact of short-scales fluctuations on the system dynamics. Assuming it this to be correct, it calls for a revision of the cosmological implications of fermion number and chirality non-conservation in finite temperature Abelian gauge theories.

In fact, the failure of the MHD to predict correctly the CS diffusion rate in the presence of a magnetic field is not surprising. The MHD is only valid at distances larger than the mean free path of particles in the plasma and thus it is not accounting correctly for the shorter scale electromagnetic fluctuations that may change the CS number and which are automatically included in our simulations. The similar phenomena exists in the diffusion of the CS number in non-Abelian theories, where the relevant distance scale $\sim\frac{1}{\alpha T}$ is much smaller than the scattering length $\sim\frac{1}{\alpha^{2}T}$. Though our lattice simulations, being classical in nature, do not include several effects - hard thermal loops (HTL) and scattering of energetic particles with energy $\sim\mathcal{O}(T)$, we do believe that they are reliable. The HTLs introduce a Debye-screening of the electric fields, plasmon masses for thermal excitations, and the Landau damping Braaten:1989mz ; Andersen:1999sf . All these effects have the energy scale $\sim eT$. In our discretised system, the similar effects are generated through lattice artifacts Bodeker:1995pp . With our parametrisation, we have a UV cutoff of order $\sim T$, which will induce corrections mimicking HTLs Bodeker:1995pp . Regarding the effect of collisions on the CS diffusion rate, we expect them to be negligible, as the longest scales which enters our problem is smaller than the mean free path. As detailed in section 2.2 of Figueroa:2017hun , the length scale associated to the chemical potential decay is the typical size of a configuration carrying CS-number one $l_{CS}\sim\frac{1}{\alpha T}$, while the mean free path in the plasma is of order $l_{mf}\sim\frac{1}{\alpha^{2}T}$. These conjectures can be checked by including HTLs and effective collision terms into classical equations of motion along the lines of ref. Bodeker:1999gx ; Rajantie:1999mp ; Hindmarsh:2001vp . We leave this for the future work.

4 Conclusions

In this work, we studied the evolution of the fermionic charge in an Abelian gauge theory at finite temperature, which can be a proxy either for fermion non-conservation or for chirality breaking. In section 2.1, we studied the dependence of the evolution on the initial value of the chemical potential. We observed that for small chemical potential, it develops long initial $plateau$s, before the decay is triggered. We saw in section 2.2 that for larger initial chemical potential, the decay happens through a self-similar process. This leads to a phenomenon of inverse cascade in the gauge field sector, with power transferred from UV to IR scales. We observed and quantified both of these phenomena, measuring the critical exponent of the chemical potential self-similar decay to be $-\frac{1}{2}$.

In section 2.3, we moved on to the study of the effect of an external magnetic field on the chemical potential evolution. Both for large and small initial chemical potential, three different situations can happen. First, when the external magnetic field is small, the chemical potential relaxes to its critical value, in the same way that in the absence of chemical potential. Then, from the critical value, it eventually decays to zero. For large external magnetic field, the dynamic is fully determined by the magnetic field, and the decay happens through damped oscillations. For moderate magnetic fields, both effects can be observed at the same time.

A way to isolate the effect of the external magnetic field influence is to study configurations with initial chemical potential equal to the critical one, as presented in 2.4. Then, only the exponential decay is visible. Following this procedure, we could extract a chiral dissipation rate and study its dependence on the parameters of the theory, finding $\Gamma_{5}=10^{-2.10\pm 0.02}B^{2}e^{11/2}/T^{3}$. We also compared this rate to the fluctuation rate of the topological charge in the absence of chemical potential, which was measured in Ref. Figueroa:2017hun . These rates are related through a fluctuation-dissipation theorem. We found a good agreement between them, providing a solid self-consistency check for our framework.

In section 3 we analysed our results in light of MHD, which is the long-wavelength effective theory of our system. Despite being out of its range of validity, qualitative features of the chemical potential can be well described. For example, in section 3.1, we provide a simplified MHD modeling which reproduces the time-scale of the observed $plateau$s in the chemical potential decay. In section 3.2 we also showed that the spectrum evolution of the inverse cascade is well-modeled by a MHD-inspired description. In section 3.3, we studied the chiral rate in the physically relevant range of chemical potentials. We showed that our numerical rates are an order of magnitude larger than the MHD predictions. We interpreted this as evidence of the effect of short-scale fluctuations on the system dynamics. This can potentially shed a new light on the role of the abelian contribution to the anomaly on fermion number/chirality violating processes.

The main outlook of this work will be to confirm or refute these results. We intend to verify numerically the effects of hard thermal loops and collisions, using some effective theory, see e.g. Bodeker:1999gx ; Rajantie:1999mp ; Hindmarsh:2001vp .

Another direction which remains to be explored is the regime of extremely small chemical potential. This is however challenging from a technical point of view, as the smaller the initial values of the chemical potential, the larger the volume needs to be, given the existence of a critical value $\mu_{c}$. Moreover, one also need longer simulations, as the initial $plateau$ gets longer. Our computer resources thus limits us to explore systematically this regime.

Acknowledgements

The authors want to thank Oleksandr Sobol for relevant comments and Kohei Kamada for useful discussions. A.F. wishes to thank Hauke Sandmeyer for giving him access to his numerical fitting routines and Pavel Buividovich for an insightful discussion. This work was supported by the ERC-AdG-2015 grant 694896 and the Swiss National Science Foundation. The numerical simulations were performed on the Intel Broadwell based cluster Fidis, provided by the EPFL HPC center (SCITAS). The discrete Fourier transforms were performed with the AccFFT library 2015arXiv150607933G .

Appendices

Appendix A Lattice set-up

The results presented in this work were obtained using the lattice discretisation presented in Ref. Figueroa:2017qmv . Initial conditions for the gauge fields and the scalar fields are drawn from a thermal ensemble, generated by a simple Metropolis algorithm. Gauss law is then enforced on the thermalised configurations. After that, the system is evolved along the classical trajectory specified by the set of discretised equations of motion

[TABLE]

with

[TABLE]

and $\Delta^{\pm}_{\mu}f=\pm\frac{1}{\rm dx}(f_{\pm\mu}-f)$, $D_{\mu}^{\pm}f=\pm\frac{1}{\rm dx}(e^{\mp ie{\rm dx^{\mu}}A_{\mu}(n\pm\frac{1}{2})}f_{\pm\mu}-f)$ the forward/backward finite difference operator and covariant derivatives. The notation $f_{a,\mu}$ means that the component $a$ of the vector field $f$ is to be evaluated at the point $\vec{n}+\hat{\mu}$, with $\hat{\mu}$ a unit displacement in the direction $\mu$, $f_{a,\mu}=f_{a}(n+\hat{\mu})$. Notice that the equations are built out of composite field so that all the fields can be expanded consistently about the same point to reproduce the continuum equations to order $O(\mathrm{d}x^{2})$. The scalar fields live on lattice edges while the gauge fields are link variables; they live on half-integer sites. The time differential operators evolve by half a step in times. More details are to be found in Figueroa:2017qmv .

The constant background magnetic field is introduced through twisted boundary conditions Kajantie:1998rz , which imposes a constant flux. To specify, we modify the periodic boundary conditions of the first component of our gauge field as follow

[TABLE]

This corresponds to introducing a constant background magnetic field in the $z$-direction of magnitude $B=\frac{2\pi n_{mag}}{(\mathrm{d}xN)^{2}}$.

For large volumes, it is essential to initiate the Monte-Carlo with configurations which already satisfy the twisted-boundary conditions (54). This can be achieved by taking as seed

[TABLE]

where the non-trivial dependence of $A_{2}$ ensures that $A$ is consistent. In more details, a gauge field which lives on a periodic manifold must transform as

[TABLE]

with potentially three different gauge transformations $\alpha_{i}$, one by direction. These gauge transformations cannot be completely arbitrary, they satisfy so-called compatibility conditions. These are the requirement that the relation of the gauge field value at a given point $x+N\cdot\hat{i}+N\cdot\hat{j}$ to its value at $x$ better not depend on whether one first goes through the $i$ or through the $j$ boundary. In other words, we have

[TABLE]

which implies

[TABLE]

The twisted boundary conditions (54) tells us that $\partial_{1}\alpha_{2}(1,j,k)=-\frac{2\pi n_{mag}}{\mathrm{d}x}$. The choice of initial field (55) - (57) is made to satisfy the compatibility conditions (59).

Appendix B Fluctuation-dissipation theorem

The dissipative chiral decay rate $\Gamma_{5}$ is related to the diffusion rate of the topological charge $\langle Q^{2}(t)\rangle=\Gamma Vt$. To derive this relation, we will use Zubarev’s formalism zubarev1974nonequilibrium , in the spirit of Ref. Khlebnikov:1988sr . We treat $\mu$ as a dynamical variable and want to understand its out of equilibrium properties. To do so, we introduce a ”local equilibrium” partition function

[TABLE]

with $\chi(t)$ some Lagrange parameter which drives $\mu(t)$ locally out of equilibrium and the factor $VT^{2}$ was introduced to make $\chi(t)$ of the same dimensions than $\mu$. The function $H$ is the Hamiltonian associated with our system. The relevant $\mu$-dependent part is

[TABLE]

This function is used to compute $\langle\mu\rangle$ at any time zubarev1974nonequilibrium and thus determines $\chi(t)$ as a function of $\langle\mu\rangle$.

[TABLE]

However, it is not stationnary, i.e. it does not satisfy Liouville’s equation

[TABLE]

with $\langle\cdots\rangle$ the usual Poisson brackets. A way to fix this is to introduce a second density matrix

[TABLE]

with $N$ such that $\rho$ integrates to one. Integrating by part and neglecting $\dot{\chi}(t)$, i.e. considering slow processes, this can be recast as

[TABLE]

with $h=VT\chi\mu-VT\int_{-\infty}^{t}e^{T\epsilon(t-t^{\prime})}\chi(t^{\prime})\dot{\mu}(t^{\prime})\mathrm{d}t^{\prime}$. Working at linear order in $h$, it becomes

[TABLE]

where $\rho_{0}=N_{0}e^{-\frac{H}{T}}$ and $\langle h\rangle_{0}$ is the average with respect to $\rho_{0}$.

To obtain a relation between the chiral decay rate and the diffusion rate of the topological charge, we compute $\langle\dot{\mu}\rangle$ using eq. (69). Using the fact that in equilibrium we have $\langle\dot{\mu}\rangle_{0}=0$ and $\langle\dot{\mu}\mu\rangle_{0}=0$, we obtain

[TABLE]

where again we assumed that $\chi(t)$ as a weak dependence on time. Now we can make use of the anomaly equation (7)

[TABLE]

with $q(t)=\frac{e^{2}}{16\pi^{2}}\int d^{3}x\,F_{\mu\nu}\tilde{F}^{\mu\nu}(x,t)$ the topological charge density. With this notation and setting taking $t=0$ as a reference time, we can define the topological charge as

[TABLE]

Inserting Eq. (63) and setting $\epsilon=0$, we then have

[TABLE]

This last quantity also appears in the diffusion rate of the topological charge. In the absence of chemical potential, $Q(t)$ follows a random walk and we can define its diffusion rate by, following the convention of Figueroa:2017hun

[TABLE]

for times $t$ much larger than a reference time $t_{0}$101010In equation (13), we implicitly set $t_{0}=0$.. We can differentiate this expression by $t_{0}$ to get

[TABLE]

Finally, plugging this in (78), we get

[TABLE]

Leading to the relation

[TABLE]