Lattice study of electromagnetic conductivity of quark-gluon plasma in external magnetic field

We study the electromagnetic (e.m.) conductivity of QGP in a magnetic background by lattice simulations with $N_f = 2+1$ dynamical rooted staggered fermions at the physical point. We study the correlation functions of the e.m.~currents at $T=200,\,250$\,MeV and use the Tikhonov approach to extract the conductivity. This is found to rise with the magnetic field in the direction parallel to it and to decrease in the transverse direction, giving evidence for both the Chiral Magnetic Effect and the magnetoresistance phenomenon in QGP. We also estimate the chiral charge relaxation time in QGP.

The Chiral Magnetic Effect (CME) is a well known anomaly-based phenomenon which can be realized in different systems with relativistic fermionic degrees of freedom [1][2][3]. The CME is the generation of a nondissipative electric current along the external magnetic field in systems with a net imbalance between the number of right-handed and left-handed fermions or nonzero chiral density.
The nonzero chiral density should be generated in order to experimentally observe the CME. In heavyion experiments the chiral density might be generated due to sphaleron transitions in the quark-gluon plasma (QGP) [4,5]. In condensed matter systems the chiral density can be generated as the result of lattice deformations [6]. Another way to generate the chiral density and to observe the CME is to apply parallel electric and magnetic fields. In this case the chiral anomaly generates the imbalance between the right-handed and left-handed fermions which leads to the CME which manifests itself through the rise of electric conductivity along the magnetic field. This CME current has already been observed experimentally in condensed matter systems [7][8][9].
Similarly to condensed matter systems, the latter mechanism can be realized in heavy-ion experiments, where colliding ions create hot QGP with deconfined relativistic quarks. In addition, in non-central collisions the QGP is affected by huge magnetic fields generated by the motion of colliding heavy ions [4]. As a result the electromagnetic (e.m.) conductivity of QGP along the magnetic field might be significantly enhanced.
Let us consider QGP in parallel electric E and magnetic B fields. Due to the axial anomaly these fields lead to the generation of a chiral density with the rate [7] where C = N c f q 2 f . The first term in Eq. (1) describes the production of chiral charge due to the chiral anomaly, while the second term stands for the decrease of chirality due to the chirality-changing processes with the relaxation time τ . Note that Eq. (1) has the stationary solution which describes the balance between anomaly based production rate and chirality relaxation processes. The chiral charge density can be parameterized by the chiral chemical potential µ 5 through the equation of state (EoS) ρ 5 = ρ 5 (µ 5 ). We use the linear response theory and consider the electric field E as a perturbation. In this limit the generated chiral chemical potential is small and the EoS reads where the χ(T, B) is a function of magnetic field and temperature. We mostly consider large magnetic fields (q f eB T 2 ), thus the chiral density is governed by the lowest Landau level degeneracy, χ ∝ eB (χ = N c f |q f | eB/2π 2 in the non-interacting approximation). The CME generates the electric current It is assumed that the magnetic field is applied along the z axis.
In addition to the CME current there is also Ohmic current in the system. The total conductivity is the sum of Ohmic and the CME conductivities σ = σ O + σ CM E . If the electric field is applied along the x axis, the Lorentz force reduces the transverse conductivity σ O xx . The σ CM E xx component is zero in this case. The decrease of σ xx in external magnetic field is called magnetoresistance. On the other hand, if electric field is applied along the magnetic field, there is no Lorentz force and magnetoresistance. At the same time σ CM E zz is a rising function of the magnetic field which can be a manifestation of the CME 1 . These facts allow one to expect that the transport properties of QGP in heavy ion collision experiments can be considerably modified by the external magnetic field. Since the transport properties of QGP are particularly important for understanding of heavy ion collision phenomenology, in this paper we are going to study the conductivity of QGP in external magnetic field.
It should be noted that the e.m. conductivity of QCD was calculated in a number of lattice studies (see for instance [10][11][12][13]). At the same time, some e.m. properties of the QGP in the presence of a magnetic background, like its magnetic susceptibility, have been already explored [14][15][16][17][18][19][20][21], as well as the emergence of anisotropies related to the magnetic background in other relevant quantities [17,[22][23][24][25]. Quenched lattice study of the e.m. conductivity of QCD in external magnetic field was carried out in [26], where no sign of neither CME nor magnetoresistance in QGP was found. We would like also to mention the lattice study of the e.m. conductivity with the external magnetic field in the Dirac semimetals [27] where the CME and magnetoresistance were observed in the semimetal phase which is similar to QGP in some properties. Finally we would like to mention lattice study of the CME in thermodynamic equilibrium [28].
In this paper we carry out the first lattice study of the e.m. conductivity of QGP in external magnetic field with N f = 2+1 dynamical staggered quarks at physical quark masses, more details about the lattice discretization and algorithms are provided in Appendix. We consider temperatures T = 200, 250 MeV for several values of the external magnetic field. Most simulations are carried out on a 16 × 64 3 lattice, with spacings a = 0.0618 fm and a = 0.0493 fm correspondingly. To check lattice spacing dependence we also consider a 10 × 48 3 lattice with a = 0.0988 fm. To study the ultraviolet (UV) properties of the correlator of two e.m. currents we consider simulations on a 96 × 48 3 lattice at a = 0.0988 fm, which corresponds to approximately zero temperature.
To study the conductivity we apply the following strategy. We first calculate the lattice correlation function where τ is the Euclidean time and J i (τ ) is the conserved current x , χ f x are staggered fermion fields of f = u, d, s flavours, and U x,i is the gauge field matrix.
The well known property of the staggered fermions is that the correlator (6) corresponds to two different operators for the even τ = 2n × a and odd τ = (2n + 1) × a slices. In the continuum limit C ij (τ ) reads where s e, o = ±1 is the timeslice parity and and ψ f is Dirac spinor of the flavour f . Notice that the operator A i corresponds to e.m. current in the continuum whereas we would like to remove the B i contribution. Next let us recall that the current-current Euclidean correlators both for even and odd slices C e, o ij are related to its spectral functions ρ e, o ij (ω) as where K(τ, ω) = cosh ω(β−τ /2) sinh ωβ/2 . The e.m. conductivity σ ij is related to the spectral densities ρ e, o ij (ω) through the Kubo formulas Notice that in last formula the contribution of the correlator B i (τ )B j (0) to the sum ρ e ij + ρ o ij cancels out and in the continuum limit the e.m. conductivity is reproduced. It is important to notice that similarly to [10][11][12][13]26] in this calculation of the correlation function (6) only connected diagrams are accounted.
Given the correlation functions C e, o ij one needs to invert the integral equation (9) and determine the spectral functions ρ e, o ij to find the conductivity. To do this one can apply the Backus-Gilbert (BG) [29] or Tikhonov regularization (TR) [30] approaches. A detailed description of these approaches can be found in Appendix. Our calculation shows that both approaches give similar results but the TR resolution function for the conductivity is a little narrower (see below). For this reason we calculate the conductivity within the TR approach.
The TR method is based on the calculation of the estimator of the spectral functionρ(ω) instead of the spectral function ρ(ω) itself. The estimator is defined as where δ(ω,ω) is the resolution function peaked around ω. If δ(ω,ω) = δ(ω −ω) the estimator of the spectral function exactly reproduces the spectral functioñ ρ(ω) = ρ(ω). However, in real calculations the resolution function has a finite width of few T . In particular, at thē ω = 0 the width of the resolution function is ∼ 3.5 T . The estimator averages the spectral function over the width of the resolution function. The TR method can be used to reconstruct the spectral function at ω = 0 if the width of the resolution function δ(ω = 0, ω) is of the order of or smaller than the characteristic variation scale of the spectral function around ω = 0, otherwise the TR method might underestimate it. Lattice data for the correlation functions of the e.m. currents are well described by either the anzats combining the transport peak at small frequencies and UV contribution at large frequencies [11][12][13] or by the AdS/CFT spectral function [13]. In the temperature interval under consideration the widths of the resolution functions are close or smaller than the variation scale of the spectral functions from [11][12][13]. For this reason we believe that the TR method can be used to calculate the e.m. conductivity in QGP. Notice also that our results for the conductivity at zero magnetic field shown in Fig. 1 are in agreement with the previous studies which hints at a correct reconstruction of the conductivity.
Another important issue is the UV contribution to the reconstructed conductivity. For instance, in the studies of shear and bulk viscosities of gluon plasma [31,32] the UV spectral density scales as ρ ∝ ω 4 , which results in a large UV contribution to the estimator (11). This contribution should be subtracted in order to obtain reliable results. For the conductivity the UV contribution scales as ρ ∝ ω 2 and our calculation shows that the UV gives ∼ 20 − 30 % contribution atω = 0. In Appendix we give a detailed description of the UV subtraction procedure.
To summarize, the calculation is done in the following steps. First we measure lattice correlation functions C e, o ij (τ ). Then we calculate the estimatorsρ e, o (ω)/ω at ω = 0 within the TR approach and subtract the UV contribution. Finally, using Eq. (10), we calculate the e.m. conductivity.
The e.m. conductivities normalized to the factor T C em (C em = e 2 f q 2 f ) at zero magnetic field and temperatures T = 200, 250 MeV are shown in Fig. 1. In addition we plot the results of [10,12]. Notably our results are in agreement with previous lattice studies within the uncertainties. Let us now consider the e.m. conductivity of QGP in the presence of the external magnetic field. From a technical point of view, the magnetic field affects directly the path-integral measure and the fermion propagators entering the construction of the e.m. currents; moreover, the e. m. U (1) phases enter the gauge links in the definition of the split current in Eq. (7). Apart from this, the problem turns out to be easier than at eB = 0. In particular, instead of the correlation functions C e, o eB we consider the difference ∆C e, o = C e, o eB −C e, o eB=0 . Since, for the chosen values of the lattice spacing, the UV regime starts at ω 0 ∼ 2 GeV, we note that q f eB ω 2 for all frequencies in the UV regime and magnetic fields. Thus, one can consider the UV spectral function magnetic fieldindependent and assume that the differences ∆C e, o do not contain the UV contribution. The results for ∆C e, o turn out to be more accurate since the UV-estimation uncertainty is absent in this case. The correlator ∆C e, o is related to additional conductivity due to the presence of the magnetic field. In our further study we apply the TR approach to the differences ∆C e, o .
The e.m. conductivity due to the external magnetic field ∆σ = σ eB − σ eB=0 normalized to T C em at temperatures T = 200, 250 MeV is shown in Fig. 2. It is seen that ∆σ rises with magnetic field for both temperatures, which is the observation of the CME in QGP on the lattice. Notice also that the rise of the ∆σ becomes linear for sufficiently large magnetic field, in agreement with Eq. (5). A similar linear growth of parallel conductivity σ with large magnetic fields was obtained in [33] by studying kinetic equations. In turn, ∆σ ⊥ is negative and decreases with the magnetic field, which is the observation of the magnetoresistance in QGP. Note also that the slopes on both functions ∆σ (eB), ∆σ ⊥ (eB) decrease with temperature. We believe that this can be explained by a decrease of the relaxation time with temperature because of the increased thermal activity. In order to estimate the finite N t effects we calculate the conductivity at T = 200 MeV, eB = 0, 0.52, 0.79, 1.12 GeV 2 on the lattice 10 × 48 3 in addition to the lattice 16 × 64 3 at hand. The results of this calculation are shown in Fig. 1 and Fig. 2. It is seen that within the uncertainties of the calculation the conductivities calculated at different lattice spacings are in agreement with each other. From this we conclude that discretization effects are under control in our study.
From equations (6), (7) it is seen that there are two contributions to the conductivity. The first one is the valence quarks contribution to the current operator (7) and the other results from the sea quarks in the fermionic determinant of the partition function. The valence quarks' contribution can be separated into u-, d-and s-quark contributions to the conductivity which can be calculated from the quark loop of the corresponding flavour 2 . In Fig. 3 we plot the u-, d-and s-quark contributions to the conductivities ∆σ , ∆σ ⊥ normalized to the factor T e 2 q 3 f at temperature T = 250 MeV. The normalization factor was chosen so as to reduce the dependence of the corresponding contribution on the quark flavour: the q 2 f results from the correlation function (6) while the additional q f results from the leading order coupling of the magnetic field to the quark q f eH. From Fig. 3 it is seen that within the uncertainty we do not see the dependence of ∆σ ⊥ /q 3 f on the quark flavour. In turn, within the uncertainty the contributions of the d-and s-quarks to ∆σ /q 3 f agree, while the contribution of the u-quark is slightly larger. This can be explained by the larger charge of the u-quark. We thus conclude that the leading dependence of ∆σ , ∆σ ⊥ on the quark flavour is proportional to q 3 f . In addition the relatively heavy s-quark mass does not influence ∆σ , ∆σ ⊥ within the uncertainties. The TR method also allows to reconstruct the spectral function, for instance ∆ρ . The reconstructed ∆ρ (ω)/f (ω), where f (ω) = ω 2 / tanh ω/2T , at T = 200 MeV is shown in Fig. 4. It is seen that the infrared part of the spectral function (ω/T < 10) has a transport peak with the height rising with the magnetic field. The ultraviolet part of the spectral function weakly depends on the magnetic field and remains close to zero as expected.
The dependence of ∆σ on the magnetic field which is responsible for the CME allows us to estimate the relaxation time of the chiral charge (see Eq. (5)). The relaxation times turn out to be τ The e.m. conductivity of QCD in external magnetic field was studied in [26] in the quenched approximation, reporting no evidence of either CME or magnetoresistance in QGP. A possible source of the disagreement is the small magnetic field used in [26], where the largest field used in the deconfinement phase is eB = 0.36 GeV 2 : at the same value our signal is quite small (see Fig. 2), so probably the signal was hardly detectable in [26]. The authors of [26] have also conducted simulations in the confinement phase and observed the rise of ∆σ and drop of the ∆σ ⊥ . Similarly, we have also calculated the conductivities in the confinement phase using the approach developed in this paper and obtained similar results: ∆σ rises while ∆σ ⊥ drops with magnetic field. However, we would like to stress that contrary to the deconfinement phase the structure of the spectral function in the confinement phase is rather complicated. For instance, it contains the contribution of the intermediate π + π − mesons or the ρ meson peak which has a large spectral weight in the confinement [12]. It is reasonable to expect that the external magnetic field modifies the spectral function, for instance through the light meson masses modification [35][36][37]. Thus, in order to check the presence of CME in the confinement phase one has to separate the contribution to the spectral function due to the conductivity ω ∼ 0 from the contribution of the light mesons ω ∼ 2m π , m ρ , .... This is a difficult task which can not be done with the data used in this paper; that might be the case for Ref. [26] as well. Note also that in [27] ∆σ was studied in Dirac semimetals both in semimetal and insulator phases. Due to chiral symmetry breaking in the insulator phase it was found that ∆σ = 0 and there is no CME in this phase. In the confinement phase there is also chiral symmetry breaking. For this reason one can expect that there is no CME in the confinement phase at sufficiently low temperature.
In conclusion, this paper is devoted to a lattice study of the e.m. conductivity of QGP in the presence of a magnetic background field. It is found that the conductivity along the magnetic field rises with the magnetic field, which is a possible manifestation of the CME. On the contrary, the conductivity in the transverse direction is decreasing with the magnetic field, which is the magnetoresistance phenomenon. Thus we observe evidence for the CME and magnetoresistance in QGP. Finally, we have also computed the relaxation time of the chiral charge in QGP for the explored temperature range. Maximal Entropy Method (MEM) is a popular method for the reconstruction of the spectral functions [38].
For the e.m. conductivity MEM was applied in papers [11,26]. It is rather difficult to carry out our study with MEM. This is because for staggered fermions we have N t /2 = 8 points(due to the periodicity of the correlator) in temporal direction which are splitted into 4 points for even time slices and 4 points for odd time slices. To conduct the reconstruction in this case you have to reconstruct separately even and odd spectral functions. We believe that for MEM this is very complicated task. Notice also that MEM can be applied only for positive spectral functions. However, this is not the case for odd branch of the spectral function in magnetic field. For these reasons we decided to apply Backus-Gilbert(BG) and Tikhonov regularization (TR) methods.
The BG and TR methods are non-parametric approaches which can be used to study the spectral function 3 . These methods are aimed at the solution of the equation where K(ω, τ ) = cosh ω(τ −β/2) sinh ωβ/2 f (ω) and f (ω) is an arbitrary function. Within the BG and TR methods instead of ρ(ω) one reconstructs the estimatorρ(ω) expressed as where δ(ω, ω) is the resolution function that has a peak aroundω and normalized. The BG and TR are the linear methods and the resolution function is taken in the form thus the estimator is a linear combination of the correlation function values Accurate reconstruction of ρ(ω) requires minimization of the width of δ(ω, ω). However, too small values of the estimator make the method unstable and susceptible to noise in the data. Thus, the method requires regularization that should be properly adjusted.
Within the BG method one minimizes the functional H(ρ(ω)) = λA(ρ(ω)) + (1 − λ)B(ρ(ω)). The term A represents the width of the resolution function: regularizes ρ(ω) making it less susceptible to noise. In terms of the covariance matrix and functions q i (ω) used to defineρ(ω) in Eq. (14), it reads B( q) = q TŜ q. Thus, statistical uncertainties are reduced at cost of increasing the width of the resolution function through decrease of λ.
The minimization of H gives the following linear functions on the form (14) The TR method is another way of the regularization of the same problem. While in the BG method the regularization is performed as W ij → λS ij + (1 − λ)W ij , in the TR scheme the SVD decomposition of W −1 = V DU T is regularized.
We remark that both in the BG and TR methods, the resolution function is an outcome of the method itself, and cannot be chosen a priori. This makes it difficult to well define a continuum limit, since there is no guarantee that the measured quantity is defined in the same way across different gauge ensemble. However, in our calculation we empirically observe that the dependence of the resolution functions on the parameters of the calculation is very weak. This is reflected in the good agreement of the results of N t = 10 and N t = 16 lattices. In the future, if a continuum limit has to be carried out, a fixed resolution function must be employed, for example following the approach suggested in [40].
As mentioned, the calculation of the electromagnetic conductivity is carried out in the following steps. Firstly, we measure the lattice correlation functions C e, o ij (τ ) (8). Then we calculate the estimatorsρ e, o (ω)/ω atω = 0 within the TR approach. For the eB = 0 case we subtract the UV contribution. Finally using Eq. (10) we calculate the electromagnetic conductivity.
In the reconstruction procedure one has to choose the value of the regularizer γ. We found that in the region γ < 1 the width of the resolution function is ∼ 3T , but the method becomes unstable what leads to large uncertainties of the calculation. At the same time in the region γ > 10 the method is stable with small uncertainties but the resolution function is rather wide (width > 4T ). This, in the region 1 < γ < 10 the method is stable and the resolution is sufficiently narrow ∼ 3.5T . For this reason we vary the regularizer in the region γ ∈ (1, 10) and use f (ω) = ω.
For zero magnetic field we subtract the ultraviolet contributions from even and odd spectral densities. To this end, we use the model for the spectral densities at large frequencies. Taking into account asymptotic freedom in QCD it is reasonable to assume that real spectral densities at ω Λ QCD do not deviate considerably from their tree level expressions. This assumption will be confirmed below. It allows us to propose the following forms of the spectral densities at large frequencies where Z e , Z o are the coefficients for the even and odd branches. At the tree level approximation Z e = 1 2 , Z o = 3 2 , but these coefficients can be renormalized by the interactions. In [32,41] it was shown that the BG method with proper scaling can be used to determine the UV coefficient such as Z e , Z o . Following [32,41], we apply BG approach with the rescaling function 4 f (ω) = 3 4π 2 ω 2 tanh ωβ 4 and use the lattice data for the correlators calculated on the lattice 96 × 48 3 . In Fig. 6 we plot the UV behavior of the rescaled reconstructed spectral functionρ(ω) for different values of λ = 10 −4 , 10 −5 , 10 −6 at the tree level case. From Fig. 6 one can see that the reconstructed asymptotic values match the tree-level values 3/2, 1/2, what confirms the validity of the method.
In Fig. 7 we perform the same procedure in the interacting case for λ = 10 −4 , 10 −5 , 10 −6 . From Fig. 7 it is 4 Notice that in order to account discretization uncertainties we use lattice expressions for the function f (ω). seen that in the UV the spectral functions indeed correspond to the models (18). Notice also that the coefficients Z e , Z o are considerably renormalized as compared to their tree level values, but their mean value is only slightly renormalized, A = (Z e + Z o )/2 ∼ 1.0. Taking into account the uncertainties of the calculation we obtain A = 1.05 ± 0.05. Notice that this result agrees with previous calculations [12,13], where the renormalization of A was shown to be small. Finally, to perform the subtraction of the UV contribution, we take the mean value of two branches conductivity ρ(ω) = (ρ e (ω) + ρ o (ω))/2. Then the UV contribution is subtracted in the form ∆ρ(ω) = A +∞ ω0 dωδ(ω, ω) 3 4π 2 ω 2 tanh where the ω 0 is the frequency which represents the asymptotic freedom region (19) onset. Unfortunately we are not able to determine the value of the ω 0 within the BG method. In the calculation of the conductivity we vary ω 0 ∈ (1.5 GeV, 3.0 GeV) and account this as the systematic uncertainty. Note that this range of ω 0 is in good agreement with the one obtained in [12] within the fitting procedure. Using formula (19) it is not difficult to find that before the subtraction the UV contribution gives 20-30% to the reconstructed conductivity.