The Chiral Separation Effect from lattice QCD at the physical point

: In this paper we study the Chiral Separation Effect by means of first-principles lattice QCD simulations. For the first time in the literature, we determine the continuum limit of the associated conductivity using 2+1 flavors of dynamical staggered quarks at physical masses. The results reveal a suppression of the conductivity in the confined phase and a gradual enhancement toward the perturbative value for high temperatures. In addition to our dynamical setup, we also investigate the impact of the quenched approximation on the conductivity, using both staggered and Wilson quarks. Finally, we highlight the relevance of employing conserved vector and anomalous axial currents in the lattice simulations.


Introduction
The non-trivial topological structure of the QCD vacuum manifests itself in a variety of fascinating aspects that can be studied theoretically and probed experimentally.In this context, anomalous transport phenomena have attracted the most attention recently.These effects appear due to the interplay between quantum anomalies and electromagnetic fields or vorticity, explaining their intimate relation to the topological nature of QCD and representing the impact of event-by-event CP-violation in QCD [1].The prime example among these phenomena is the Chiral Magnetic Effect (CME): the generation of a vector current due to magnetic fields and a chiral imbalance [2].In the last decade, major experimental campaigns have been launched in order to measure the CME in different setups ranging from condensed matter systems [3] to heavy ion collision experiments [4][5][6].
Another anomalous transport phenomenon is the Chiral Separation Effect (CSE): the generation of an axial current in a dense and magnetized system [7,8].Although highly analogous to the CME in its formulation, the equilibrium interpretation of the two effects turns out to be quite different.Together with the CME, it is expected to form the Chiral Magnetic Wave (CMW), a collective excitation emerging in a magnetized environment [9],  which is also the subject of intense experimental searches [10].The key parameter, describing the strength of the effect is the conductivity coefficient C CSE , to be defined in detail below.
The analytical approaches to assess the CSE span a broad range and include for example chiral kinetic theory [11], effective models of QCD [12], holography [13] or the field correlator method [14].The first calculations were performed for massless quarks in the absence of gluonic interactions, giving C CSE = 1/(2π 2 ) [7].It was later recognized that this value is not fixed by the axial anomaly but rather depends on the quark mass m as well as the temperature T .In fact, as we show in App.A, for non-interacting quarks, an analytical treatment of this problem gives in agreement with [7,8].The resulting curve is shown in Fig. 1, revealing a suppression of the conductivity coefficient for heavy quarks (or, equivalently, low temperatures), see also Refs.[15,16].Note that this result corresponds to the linear behavior of the CSE current, valid for small chemical potentials, and this is the domain that we aim to study in this paper.
Given the sensitivity of the CSE to the mass and the temperature, it is natural to expect that the impact of gluonic interactions on C CSE will also be substantial (see also Ref. [17]).In the phenomenologically interesting region around the finite temperature QCD crossover [18], these interactions are non-perturbative.The most successful nonperturbative tool, based on the first principles of the underlying quantum field theory, is lattice QCD.However, at nonzero quark density -which is required for the discussion of the CSE -standard lattice simulations break down due to the complex action problem, see e.g. the recent review [19].For this reason, many of the existing lattice investigations [20,21] were not performed in full QCD, but either in the two-color version of it, or in the so-called quenched approximation -both of which are free of the complex action problem.Even with this caveat in mind, some of the findings in the literature gave different predictions about whether QCD interactions at intermediate temperatures affect the value of C CSE [21] or not [20].The lattice regularization has also been employed to study the CSE for free fermions [22,23] and with classical-statistical real-time simulations [24,25].
In this work, we demonstrate for the first time how the CSE conductivity can be constructed using a Taylor-expansion around zero density, thereby enabling its treatment in full QCD.We determine C CSE using dynamical staggered quarks with physical masses and carry out its continuum extrapolation for a broad range of temperatures in the transition region as well as deep in the confined phase.We shed further light on the impact of the quenched approximation and determine C CSE in this setup using both staggered and Wilson quarks.In this setup, we also discuss the effect of increasing the quark masses away from the physical point.Finally, we also determine the CSE conductivity in the absence of gluonic interactions for both the staggered and the Wilson fermion discretizations -this provides a useful benchmark of the appropriate lattice currents.
The present effort will be useful for studying further anomalous transport effects like the CME, for which most lattice simulations so far [23,[26][27][28][29][30] were either based on indirect approaches or are yet to be performed in full QCD at the physical point.

Conductivities and lattice methods
As introduced above, the chiral separation effect amounts to the generation of an axial current J ν5 in the presence of a background magnetic field and nonzero quark density.The latter is parameterized by a chemical potential µ that couples to the corresponding conserved number density.The magnetic field B is considered to be constant and homogeneous and, without loss of generality, pointing along the third spatial direction.Moreover, the magnetic field is measured in units of the elementary electric charge e > 0, so that we can work with the renormalization group invariant combination eB.In this paper, we are interested in the leading-order behavior of the current for weak chemical potentials and weak magnetic fields.This response is described by the linear Taylor-coefficient C CSE that we refer to as the CSE conductivity coefficient, or simply as CSE conductivity.

Currents and chemical potentials
For a single quark flavor, the definition of the axial current and the chemical potential is unambiguous.For several quark flavors, like in full QCD, there are different options which -as we will find below -give different results for the strength of the effect.One may consider the axial current with different quantum numbers i = B, Q or L in flavor space, where f = u, d, s, . . .labels the quark flavors and q f are the corresponding electric charges.The currents (2.1) are rendered intensive by dividing by the Euclidean four-volume V /T .The above choices correspond to the currents coupling to baryon number B, electric charge Q or the so-called light baryon number L [31].The latter serves to mimic heavy-ion collision setups.
Analogously to the axial currents, the quark density can also be controlled either by a baryon chemical potential µ B , an electric charge chemical potential µ Q or a light baryonic chemical potential µ L , which couples to the corresponding vector currents, 3) The chemical potentials for the individual quark flavors are µ f = µ B /3 in the baryonic and µ f = µ Q q f /e in the charge case, while they are set as for the light baryon chemical potential.
The most common choice in the flavor space corresponds to the quantum number Q for the current and B for the chemical potential.In this case, the CSE conductivity coefficient is defined via the leading order-behavior for the expectation value Here we use the short notation C CSE ≡ C QB CSE and we analogously define C ij CSE , where the first superscript i ∈ {B, Q, L} corresponds to the current and the second one j ∈ {B, Q, L} to the chemical potential.Note that the conductivity of Eq. (2.4) -as well as the related conductivities C ij CSE -are defined through the leading-order response of the current to the chemical potential and the magnetic field at physical quark masses and are, generically, functions of the temperature.For free quarks at zero temperature, for instance, C CSE vanishes, as discussed in the introduction.
Full lattice QCD simulations suffer from the complex action problem at µ B ̸ = 0, µ Q ̸ = 0 or µ L ̸ = 0. To circumvent this issue, we consider the Taylor-expansion of the current expectation value in the chemical potential.To leading order, this gives the first chemical potential derivative of the CSE current, which, using Eq.(2.4), yields and analogously for the other flavor combinations.This, leading-order Taylor expansion of the current only requires simulations at zero density, free of the sign problem.Finally, C CSE can be extracted via numerical differentiation with respect to eB.We note that, alternatively, one could also perform the derivative with respect to the magnetic field analytically.On the lattice, such derivatives can be implemented using the methodology developed in Ref. [32], see also Refs.[33,34].In the present case, this would lead to an expression including a three-point function with two vector currents and one axial current, revealing the relation between the conductivity and the anomalous triangle diagram.
For each choice of the flavor structures, there is an overall proportionality constant involving the quark baryon numbers or quark charges and the number of colors N c = 3, which can be factored out of the conductivity.Since the coupling to the magnetic field always occurs via the quark charges, these overall factors read and similarly for the remaining ones.For massless, non-interacting quarks, these factors appear as overall normalizations in C CSE .(Note that C BB dof = 0 for three light quarks.)Note, however, that for massive, non-degenerate quarks each flavor contributes differently, depending on the choice of the quantum numbers.In the interacting case, disconnected diagrams are also different for B, Q or L. This implies that the individual C ij CSE are not equivalent despite scaling out the above factors.Since these overall factors can always be restored, from this point we rescale all our results by C dof unless explicitly stated otherwise.
Throughout this paper, we are working in Euclidean space-time.Taking into account the relation between the Minkowski (indicated by the superscript M) and Euclidean Dirac matrices (γ M 0 = γ 4 , γ M i = iγ i ), we find that the Minkowski-space observable reads Therefore we consider, for the rest of the text, the imaginary part of the Euclidean twopoint function calculated on the lattice.The real part of the Euclidean observable was checked to be consistent with zero.

Lattice setup -staggered quarks
We begin by describing the setup with rooted staggered quarks, which we used to study the CSE in full dynamical QCD with up, down and strange quark flavors.Here, the partition function Z is written using the Euclidean path integral over the gluon links U as with β = 6/g 2 denotes the inverse gauge coupling and m f the quark masses for each flavor f = u, d, s.In Eq. (2.8), S g is the gluonic action, which in our setup is the treelevel Symanzik action, and M f is the massive staggered Dirac operator with twice stoutsmeared links.The Dirac operator contains the quark charges q u /2 = −q d = −q s = e/3.The quark masses are tuned to the physical point as a function of the lattice spacing a [35].The simulations are carried out using four-dimensional lattice geometries with N s spatial and N t temporal points.The physical spatial volume is given by V = L 3 = (aN s ) 3  and the temperature by T = (aN t ) −1 .The lattice sites are labeled by the four integers n = {n 1 , n 2 , n 3 , n 4 } and ν denotes the unit vector in the ν direction.
The magnetic field B is included as a classical background field, so we do not consider dynamical photons in our setup.The electromagnetic potential enters the Dirac operator for the quark flavor f in the same way as the SU(3) links, i.e. as a parallel transporter between two lattice points u f ν = exp(iaq f A ν ).These electromagnetic links are chosen such that they represent a homogeneous magnetic field in the x 3 direction.The flux of the field is quantized as eB = 6πN b /(aN s ) 2 with the integer flux quantum N b ∈ Z [36].Moreover, the chemical potential is also included in the exponential form, multiplying the temporal links by exp(±aµ f ).Even though our simulations are at µ f = 0, we will need this form to calculate the Taylor-expansion coefficients.
Within the staggered formalism, the quark fields ψ f are transformed in coordinate and Dirac space to the so-called staggered fields χ f in order to partially diagonalize the Dirac operator.To proceed, we need to express the vector and axial currents in terms of χ f .The prescription that maintains a conserved vector current and an anomalous (flavor-singlet) axial current are the point-split bilinears [37], involving the staggered counterparts of the Dirac matrices [38], (2.10) Here, η ν (n) = (−1) ρ<ν nρ are the staggered phases and ϵ νραβ the totally antisymmetric four-index tensor with the convention ϵ 1234 = +1.We mention that these Dirac structures depend explicitly on the links as well as on the magnetic field and the chemical potential.
In particular, note that Γ f 35 involves hoppings in the temporal direction and thus depends on µ f (but not on the chemical potentials for the other quark flavors With these definitions, we can now perform the Grassmann integral over the staggered fields, giving the expectation value for the (volume-averaged) axial current, where Tr refers to a trace in color space and a summation over the lattice coordinates.The factor 1/4 results from rooting.Next, we can calculate the derivative (2.5) required for the CSE conductivity.This derivative generates the usual disconnected and connected terms, together with an additional tadpole term arising due to the derivative of Γ f 35 with respect to µ B , where the flavor coefficients (2.2) entered.We emphasize that the tadpole term -the last line in Eq. (2.12) -appears due to the exponential form of introducing the chemical potential.This form is required to maintain gauge invariance on the lattice and to avoid chemical potential-dependent ultraviolet divergences [39].The tadpole contribution is analogous to the one that appears in the calculation of quark number susceptibilities for staggered quarks, see e.g.Ref. [40].We also mention that in deriving (2.12) we exploited the charge conjugation symmetry of the µ B = 0 system, i.e. that the expectation value of the baryon density and that of the axial current vanish at µ B = 0 (we refer back to this point in Sec.3.2).The expectation values in Eq. (2.12) are to be evaluated at µ B = 0 but nonzero magnetic field eB.The conductivities C ij CSE for the other flavor quantum numbers can be calculated analogously, and differ from Eq. (2.12) by factors of c i f c j f ′ in the disconnected and c i f c j f in the connected and tadpole terms under the flavor sums.To enable a direct comparison to existing lattice results in the literature, we also consider QCD in the quenched approximation.This amounts to dropping the fermion determinant from the partition function (2.8), while leaving the fermionic operators in the observables unchanged.This results in a non-unitary theory where quarks behave differently in operators as in loop diagrams.It is a commonly employed approximation that simplifies simulation algorithms significantly.For completeness, in the quenched case we employ both the staggered and the Wilson discretization of fermions.
Finally, we also calculate C CSE in the absence of gluonic interactions, both for Wilson and for staggered quarks.To this end, we calculated the eigenvalues and eigenvectors of the staggered Dirac operator exactly and constructed the necessary traces from these (see App. B), while for Wilson fermions we relied on stochastic techniques to estimate the traces.

Lattice setup -Wilson quarks
Next we describe our setup involving Wilson quarks.This discretization we only consider for free fermions and in the quenched approximation -simulations with non-degenerate light Wilson quarks (differing in their electric charges) would be a computationally much more challenging setup.In this case, the currents are constructed from the bispinor fields Here the Dirac structure is the same as in the continuum, (apart from the Wilson term), but the conserved vector and anomalous axial currents again involve a point-splitting, just as for staggered quarks [41] Γ where r is the coefficient of the usual Wilson term in the Dirac operator (taken to be one in our setup).Notice the presence of the Wilson term in the vector current and its absence in the axial current1 .Using the above definitions, the expectation value of the axial current reads, Since Γ f 35 this time only involves hoppings in the 3 direction, it is independent of the chemical potential.Thus, in comparison to the staggered case, our observable takes a simpler form, involving only disconnected and connected terms.The conductivities for the other flavor quantum numbers follow similarly.
In the literature, it is also customary to consider a local vector current instead.Even though it is not conserved, it has the same quantum numbers as J f ν and is often employed in Wilson fermion simulations, such as for the study of various aspects of hadron physics at zero and non-zero temperature, for instance.It was also used to study the CME in quenched and dynamical QCD [28].It is also possible to introduce a chemical potential µ loc B that couples to these local currents in the action.The analogue of (2.15) now reads This introduction of a chemical potential is a naive extension of the continuum formulation to the discretized theory and it is well known that this leads to new ultraviolet divergences [39].Therefore we use this setup with the non-conserved, local current merely for comparison, and we will test the impact of these ultraviolet divergences on the observable (2.17) below.

Results
To estimate the traces appearing in Eqs.(2.12), (2.15) and (2.17), we use the standard noisy estimator technique.
We performed our numerical simulations using 100 Gaussian noise vectors, which was found to be sufficient to reliably calculate the observables.For free staggered fermions, we calculated the traces using the exact eigensystem instead, see App.B for details.
To determine C CSE , we calculated the current derivative for different magnetic fields, and extracted the coefficient from its slope with respect to eB.We considered a linear fit function with no offset term: a one parameter fit whose optimal value is found via a usual χ 2 minimization method.The error analysis for each simulation is performed using the jackknife method with 10 bins, which is also used to propagate the error to the fit.This is what we refer to as the statistical error of C CSE .Furthermore, we repeat the fit, successively eliminating data points at the largest value of eB, until we are left with only one point.The largest difference in the slopes between the original fit and the different repetitions is what we consider the systematical error of C CSE .
In Fig. 2 we show an example of the obtained results and the fits, both for free fermions and full QCD with staggered quarks.In these plots, we can see that the expected linear behavior is confirmed.Now we turn to the precise analysis of C CSE in different situations.

Free quarks
Since the free case is analytically solvable, we can use this simple situation as a check of our setup.We consider one color-singlet fermion with charge q and mass m -this implies that the flavor quantum numbers are irrelevant and it is sufficient to treat a single chemical potential µ and axial current J 35 .The overall proportionality constant is thus C dof = (q/e) 2 , which is used to normalize the results.
In Fig. 3 we show the results for free staggered fermions at different m/T values.To compare to the analytic formula, both the continuum limit (a → 0), as well as the thermodynamic limit (L → ∞) need to be taken.The different panels of Fig. 3  The dashed line represents the optimal fit, while the band is our estimate for its uncertainty, obtained by adding the statistical and systematical errors in quadrature.
to the continuum extrapolations for different aspect ratios LT by increasing N s and N t (with LT and m/T kept constant).For small values of m/T , finite-size effects are expected to be sizeable, and the results confirm this expectation.In particular, LT = 4 is already found to agree with the infinite volume limit for all masses.A cross-check with Wilson fermions is also shown for the largest value of m/T .The main message that we get from this calculation is that C CSE approaches the value given by Eq. (1.1) when the volume is sufficiently large and the continuum limit is taken.In the case of Wilson fermions, we also learn that the correct setup is to consider a conserved vector current and the anomalous axial current.If instead we use the local vector current (2.16), the analytical result is no longer recovered and the continuum limit shows a possibly divergent behavior.A very similar phenomenon occurs for staggered fermions if we exclude the tadpole term of (2.12), again leading to a non-conserved (in fact ultraviolet divergent) vector current.Having cross-checked the consistency of our setup with the analytical result, we can proceed toward the physical setting by turning on the gluonic interactions.

Quenched QCD
Before moving on to full QCD, there is an intermediate step where one can analyze the impact of gluons on the CSE: the quenched approximation.Furthermore, as we will see, the quenched theory reveals an interesting feature of the CSE due to the presence of an exact center symmetry.To appreciate this, next we briefly introduce the notion of the center symmetry and the Polyakov loop.
In the absence of dynamical fermions, the system undergoes a first-order phase transition from confined to deconfined matter at around T q c ≈ 270 MeV [42], for which the

Polyakov loop
acts as the order parameter.This behavior arises from the Z(3) center symmetry of the gauge action, which is invariant under transformations U 4 (n) → zU 4 (n), with z ∈ Z(3), while the Polyakov loop is not, enforcing a vanishing expectation value of P .Below the transition temperature, center symmetry is intact and on typical configurations P = 0.In turn, in the deconfined phase center symmetry is spontaneously broken, and the theory chooses one of the three vacua with arg P = 0, ±2π/3 with equal probability.
The Polyakov loop sectors at high temperature (where these are spatially approximately homogeneous) can be mimicked by an imaginary baryon chemical potential iµ B /T = 0, ±2π/3.
As we show in App.A, nonzero imaginary chemical potentials give a non-trivial contribution2 to the observable (2.5).In Fig. 4 In full QCD, the quarks explicitly break the Z(3) symmetry and, consequently, the theory is always in the real Polyakov loop sector -corresponding to a vanishing imaginary baryon chemical potential.This indicates that the relevant QCD contribution to the CSE originates from configurations with real Polyakov loops.Therefore, for the quenched analysis of C CSE below, we rotate all our gauge configurations to this sector by the appropriate center transformation.We note that while this procedure has a clear impact at high temperatures, where the Polyakov loop background is homogeneous (and thus equivalent to a nonzero iµ B ), around T q c the non-trivial spatial distribution of center clusters (see e.g.[43]) complicates the interpretation.This is a shortcoming of the quenched approximation for the CSE.
Having this caveat in mind, we present the results for C CSE in the quenched approximation in Fig. 5.We use configurations generated with the plaquette gauge action at β = 5.845, 5.9, 6.0, 6.2, 6.26 and 6.47, and use both staggered and Wilson fermions in the valence sector.We employed these gauge configurations already in Refs.[44,45].For the staggered discretization, the quark masses were tuned to a pion mass of M π ≈ 415 MeV, while for Wilson fermions we use M π ≈ 710 MeV.The latter choice is motivated by the giving rise to an additional term ∝ Tr(Γ (2.12) for example.This term is included for the comparison in Fig. 4.
3 That is to say, the imaginary part of the Euclidean observable, see Eq. (2.7).
CME study [28].Although here we only present results at a few temperatures, we can clearly observe two distinct regimes: in the confined phase the CSE is severely suppressed, reaching zero toward T ≈ 0, while at temperatures well above T q c , the result approaches the massless free fermion value.One may compare this with the findings of the quenched study [20], where no corrections due to QCD interactions were found.However, in that case the quenched configurations were probed with a massless overlap Dirac operator in the valence sector, complicating the interpretation of the results.We emphasize moreover that according to our results, the transition between the regime where CSE is suppressed and the one where the massless case is approached, appears to occur in the vicinity of T q c ≈ 270 MeV.All these hints given by the quenched results will be confirmed by the full QCD simulations, which we present next.

Full QCD at physical quark masses
Finally, we present the main result of this study: the conductivity coefficient C CSE in full QCD, in particular for N f = 2 + 1 flavors of staggered fermions at physical quark masses.This is the first fully non-perturbative result for C CSE at the physical point.The measurements were performed on an already existing ensemble of configurations for different magnetic fields [36,46].In Fig. 6, we present the dependence of the conductivity on the temperature using several finite-temperature lattice ensembles 24 3 × 6, 24 3 × 8, 28 3 × 10, 36 3 × 12 as well as two zero-temperature ensembles 24 3 × 32 and 32 3 × 48.The impact of the temperature on C CSE is very pronounced, as we have already seen in the quenched results.Above the QCD crossover temperature T c ≈ 150 MeV, C CSE is found to approach the value corresponding to free massless quarks, in accordance with asymptotic freedom.In turn, at temperatures below T c , the conductivity decreases until reaching zero.This indicates a suppression of the CSE at low temperatures, which is consistent with a previous study in two-color QCD [21].
In the low temperature region, we consider a simple model to describe the system: a non-interacting gas of hadronic degrees of freedom.For this particular observable, only electrically charged hadrons contribute that couple to chirality (i.e. the Dirac indices of γ 5 ) 4 .Thus, in our model we only consider a gas of protons and Σ ± baryons.The contribution of heavier charged spin-1/2 baryons is negligible in the relevant temperature range, and we do not include spin-3/2 baryons either.The so constructed model also features a strong suppression of C CSE in the confined regime and is found to agree with the lattice results for T ≲ 120 MeV.It is also in qualitative agreement with other approaches to the CSE [11,14].
The continuum limit is carried out using the lattice data in the range of temperatures 90 MeV ≲ T ≲ 400 MeV.To reliably extrapolate to the continuum5 , we consider a spline fit procedure combined with the continuum limit.For this, we consider a T -dependent spline fit of all lattice ensembles with a-dependent coefficients.The best fitting surface in the a − T plane is found by minimizing χ 2 /dof.The details of the spline fit can be found in [50].The statistical error of this procedure is calculated using the jackknife samples, while we estimate the systematic error by repeating the spline fit removing the coarsest lattice.The maximum difference between the limits obtained with these two data sets is taken as the systematic error, which is added in quadrature to the statistical one.
We further provide a complete parameterization of C CSE , to be used for comparisons with effective theories of QCD and in (anomalous) hydrodynamic descriptions of heavy-ion collisions.The parameterization, which is also shown in Fig. 6, smoothly connects the continuum extrapolation of the lattice results with the perturbative limit and the lowenergy model.The details of the parameterization are given in App. C. In summary, we conclude that the CSE conductivity is an observable very sensitive to the finite temperature QCD crossover, and it may even be used to define its characteristic temperature.In particular, the inflection point of the parameterization is found at T c ≈ 136(1) MeV.As an alternative definition, we also consider the temperature, where the conductivity assumes half its asymptotic value.This gives T c ≈ 149(4) MeV.
The results above correspond to  charged axial current to the baryon chemical potential.Next, we compare results for other combinations of the flavor quantum numbers Q, B and L. The corresponding proportionality factor C dof is different for each (see Sec. 2.1), but even after normalizing the conductivities by these factors, the results are not equivalent, since the different quark flavors interact with each other via gluons.Out of the nine possibilities, in Fig. 7 we present the continuum extrapolated results for the six combinations that are most interesting from a phenomenological point of view.
Except for the BB case, all combinations show similar trends, being zero at low T and approaching the massless free fermion result at high temperature.Still, differences at intermediate temperatures are visible, in particular in the cases with a baryonic axial current.The BB setup is exceptional, because the three-flavor perturbative result in this case vanishes (see the discussion in Sec.2.1).Our full QCD result for C BB CSE also approaches zero for high temperatures.While at intermediate T , the perturbatively expected cancellation between the contributions of the three flavors is only partial, the result is still found to be much smaller than for the other quantum numbers.Moreover, the data for C BB CSE is rather noisy due to the disconnected contribution in Eq. (2.12), especially around T c .While our current results do not enable us to resolve this temperature dependence, they do show that the continuum limit of this conductivity lies within the yellow band shown in the right panel of Fig. 7.

Summary and outlook
In this paper, we have presented a comprehensive study of the Chiral Separation Effect using first-principles lattice QCD simulations.In particular, we have calculated the conductivity coefficient C CSE using dynamical staggered quarks with physical quark masses as well as quenched staggered and quenched Wilson quarks.In addition, we calculated C CSE in the absence of gluonic interactions both in the continuum and on the lattice with staggered and Wilson quarks.
In the free case, our analytic formula (1.1) reproduces the known results in the literature [7,8].Moreover, we find that both lattice discretizations reproduce the analytical result once the continuum limit is taken and the volume is sufficiently large.We highlight that the use of conserved vector currents and anomalous axial currents on the lattice is crucial -non-conserved vector currents are demonstrated to lead to incorrect results.
Having cross-checked our setup, we moved on to full QCD, where we have determined the continuum limit of C CSE for a broad range of temperatures at physical quark masses.At temperatures higher than T c , the conductivity coefficient approaches the prediction for massless free fermions, as expected, while at temperatures below T c , the conductivity goes to zero, revealing a strong suppression of the CSE in the confined phase.The result depends on the flavor quantum numbers used for the axial current and the chemical potential.Even after rescaling the conductivity by the multiplicative factor C dof containing the degrees of freedom, for example a charge chemical potential induces slightly different baryonic or charge axial currents.
A widely employed approximation to take into account gluonic interactions is quenched QCD, which we also explored in this paper.In the quenched theory, we highlighted the impact of the Polyakov loop sector on the CSE at high temperatures, an important subtlety that needs to be treated carefully when working in this approximation.Only the real Polyakov loop sector delivers the physical result -configurations in imaginary Polyakov loop sectors deviate from it and also come with drastically enhanced fluctuations.This issue can be overcome by rotating the configurations to the real sector by a center transformation.The results are in qualitative agreement with the picture based on our full QCD results: a strong suppression of the CSE at low temperatures and a tendency to the massless free fermion result at temperatures much higher than T c .
Our final results in full QCD are interpolated by a simple parameterization, explained in App.C, for use in hydrodynamic descriptions of heavy-ion collisions.Finally, we note that our technique can be generalized to study the CME, which we plan to carry out in an upcoming publication.Performing the analogous simulations for that case will contribute to a better theoretical understanding of anomalous transport phenomena, as a counterpart to the large-scale experimental efforts being made to detect this effect.

A Analytical treatment of the CSE coefficient for free quarks
We will now shortly summarize the continuum calculation leading to Eq. (1.1) for noninteracting fermions.As mentioned in the main text, it is sufficient to consider one colorless fermion flavor (with mass m and charge q), a single chemical potential µ and axial current J 35 .To enable a direct comparison to the lattice results, we will evaluate the µ-derivative of the axial current in a homogeneous magnetic field background, in which case free fermion propagators are known exactly [51,52].We will carry out the calculation starting out in Minkowski metric and extending to finite temperature during the calculation using the Matsubara formalism.Accordingly, the Dirac matrices are the Minkowski ones fulfilling {γ µ , γ ν } = 2η µν = 2 diag(1, −1, −1, −1) contrary to the main text where Euclidean Dirac matrices are used.
Since the chemical potential, µ couples to the four-volume integral of the zeroth component of the vector current J µ = ψγ µ ψ, the µ-derivative of the z-component of the chirality current can be written as We regularize in the ultraviolet using the Pauli-Villars (PV) scheme, meaning that all diagrams are replicated by the regulator fields with coefficients c s and masses m s .The

B CSE with free staggered fermions
In this appendix, we illustrate how to calculate C CSE for free staggered quarks using explicitly the eigenvalues and eigenvectors of the Dirac operator.Therefore, here we turn off the gluon fields and consider one color component only.Just as in App.A, we discuss a single quark flavor with mass m, one chemical potential µ and axial current J 35 .Below we work on an N 3 s × N t lattice and write everything in lattice units, setting a = 1.In the non-interacting theory, the disconnected term of Eq. (2.12) vanishes, so we only need to calculate two terms and the expectation value indicating the fermion path integral was omitted for brevity.Furthermore, for convenience we redefine the staggered phases as The massless staggered Dirac operator is antihermitian, / D † = − / D, therefore its eigenvalues are purely imaginary.Moreover, due to staggered chiral symmetry, { / D, η 5 } = 0 (with η 5 = (−1) n 1 +n 2 +n 3 +n 4 ), the eigenvalues come in complex conjugate pairs, In addition, we will need the analogous eigensystem for where each eigenvalue is doubly degenerate due to (B.3).Using this eigensystem as basis, the traces in Eq. (B.1) can be written as where we inserted a complete set of eigenstates j Ψ j Ψ † j = 1 in the first term and used Next, we make use of the separability of the problem and reduce Eq.(B.4) to one twodimensional and two one-dimensional eigenproblems.This will enable us to determine the complete spectrum on much larger lattices than with a direct, four-dimensional approach.To this end, we write the free Dirac operator as / Similarly, the staggered Dirac matrices (2.10) simplify to Notice that the hop operators satisfy the property due to the way that the chemical potential enters the temporal hoppings.
The squared operator M M † = M M † 12 + M M † 3 + M M † 4 separates into three terms that act in the respective subspaces and only depend on the respective coordinates -thus they commute with each other.Therefore, the eigenvectors in (B.4) factorize, 6 Below we use a shorthand notation and simply write Ψ i = ρ i ϕ i ξ i .
We proceed by expanding the operators appearing in the matrix elements in Eq. (B.5) into separate contributions that depend only on n 1,2 , n 3 or n 4 , Γ 4 = (−1) n 1 +n 2 S 4 , (B.11) with We employ the LAPACK library to separately solve the two-dimensional eigenvalue problem in the 12 plane in the presence of the magnetic field, as well as two one-dimensional problems in the 3 and 4 directions, both without electromagnetic phases.The eigensystem of the former problem gives Hofstadter's butterfly [54], the spectrum of a well-known solid-state physics model.Its relevance for lattice QCD has been pointed out in Ref. [55] and it has been generalized to the presence of gluonic interactions [56,57] as well as for inhomogeneous magnetic fields [58].The final value for the current derivative can be reconstructed using Eq.(B.18), and C CSE can be extracted in the same way as explained in the main text by calculating the observable at different values of the magnetic field.

C Parameterization of the CSE conductivity coefficient
In this appendix we provide a parameterization of the CSE conductivity in QCD as a function of temperature.We consider the continuum extrapolated data for C QB CSE in Fig. 6 and we smoothly interpolate to the prediction given by our model of a gas of protons and Σ ± baryons at low T , while at high temperatures we impose an asymptotic approach to the result for a gas of free massless fermions (m/T = 0 in Eq. (1.1)).We found the following parameterization to adequately describe the data,

Figure 1 :
Figure 1: The CSE conductivity coefficient as a function of m/T for non-interacting quarks.The marked points correspond to m/T = 0, 1, 4, which will be the values we compare to in Sec.3.1.

Figure 2 :
Figure2: Derivative of the axial current with respect to the chemical potential as a function of the magnetic field for a 24 3 × 6 lattice with staggered quarks for free fermions (left) and in full QCD with 2 + 1 flavors and physical quark masses at T = 305 MeV (right).The dashed line represents the optimal fit, while the band is our estimate for its uncertainty, obtained by adding the statistical and systematical errors in quadrature.

Figure 3 :
Figure 3: CSE conductivity coefficient for free fermions at different values of m/T using the staggered discretization.For m/T = 4, Wilson fermions as well as discretizations with a non-conserved vector current are also included.The dashed lines correspond to the marked values from Fig 2.

3 Figure 4 :
Figure4: Results for the Polyakov loop (left) and the derivative of the current with staggered fermions (right) in an ensemble of 100 configurations on a 32 3 × 8 lattice at T ≈ 400 MeV and at a magnetic field of eB = 0.74 GeV 2 .The three sectors can be clearly distinguished at this temperature, and the contribution of configurations with imaginary Polyakov loops is observed to be very different as compared to those in the real sector.

Figure 5 :
Figure 5: Results for C CSE with Wilson staggered fermions in the quenched approximation, measured on configurations generated using the plaquette gauge action.The configurations were rotated to the real Polyakov loop sector (for details see the text).The valence pion mass is set to M π ≈ 415 MeV for staggered quarks and M π ≈ 710 MeV for Wilson quarks.The quenched critical temperature T q c is indicated by the dashed vertical line.

Figure 6 :
Figure6: C CSE for a broad of temperatures in full QCD with 2 + 1 flavors of staggered fermions.The continuum limit extracted via the spline fit procedure described in the text, together with its uncertainty, is indicated by the orange band.The result of our low-energy model involving p and Σ ± baryons is shown by the continuous line at low T , while the parameterization for C CSE in the whole range of temperatures is displayed as a dashed-dotted line.

Figure 7 :
Figure 7: Continuum limits of C CSE for different flavor quantum numbers (as defined in Eq. (2.2)) for the axial current and the chemical potential.We present the phenomenologically relevant cases involving an electrically charged axial current (left panel) and a baryon axial current (right panel).

2 µ 3 Figure 8 :
Figure 8: The CSE conductivity as a function of m/T for non-interacting quarks, at different values of the imaginary chemical potential µ I .The averaged result corresponds to averaging over the three center sectors, that is µ I /T = 0, 2π/3, −2π/3.

Table 1 :
2• exp −h/t 4 • 1 + a 1 /t + a 2 /t 2 + a 3 /t 3 1 + a 4 /t + a 5 /t 2 + a 6 /t 3 , with the parameters included in Table.1.As visible in Fig.6, this function captures all details of our results -more specifically, the confidence interval of the parameterization is about 95%.Parameters of the function (C.1).