Local CP-violation and electric charge separation by magnetic fields from lattice QCD

We study local CP-violation on the lattice by measuring the local correlation between the topological charge density and the electric dipole moment of quarks, induced by a constant external magnetic field. This correlator is found to increase linearly with the external field, with the coefficient of proportionality depending only weakly on temperature. Results are obtained on lattices with various spacings, and are extrapolated to the continuum limit after the renormalization of the observables is carried out. This renormalization utilizes the gradient flow for the quark and gluon fields. Our findings suggest that the strength of local CP-violation in QCD with physical quark masses is about an order of magnitude smaller than a model prediction based on nearly massless quarks in domains of constant gluon backgrounds with topological charge. We also show numerical evidence that the observed local CP-violation correlates with spatially extended electric dipole structures in the QCD vacuum.


Introduction
Quantum Chromodynamics (QCD) is the theory of the strong interactions. At low temperatures QCD is confining, implying that the elementary particles of the theory -quarks and gluons -only exist as components of bound states (hadrons). The asymptotic freedom property of QCD ensures that at high temperatures the interaction between quarks and gluons weakens, and a transition to the quark-gluon plasma (QGP) phase occurs, where the dominant degrees of freedom are no longer colorless bound states but colored objects. According to lattice simulations, this transition is no real phase transition but an analytic crossover [1] and takes place at around T c ∼ 150 MeV, see e.g. refs. [2,3].
The high-temperature QGP phase is routinely produced in contemporary high energy heavy-ion collisions, for example at the Relativistic Heavy Ion Collider (RHIC), where temperatures exceeding T c can be reached [4]. Besides extreme temperatures, another interesting feature of such a heavy-ion collision is the presence of strong magnetic fields generated by the spectator particles in non-central events. This magnetic field is perpendicular to the reaction plane and may reach values up to √ eB ∼ 0.1 GeV for RHIC and √ eB ∼ 0.5 GeV for the Large Hadron Collider (LHC) [5], depending on the beam energy and centrality. Even though the generated magnetic field has a very short lifetime, of the order of 1 fm/c, this magnetic 'pulse' coincides with the formation of the quark-gluon plasma and thus may play an important role in the description of the collision. Strong magnetic fields also represent an important concept for cosmology [6] and for the description JHEP04(2014)129 of dense neutron stars called magnetars [7]. Therefore, a clear theoretical understanding of the response of QCD matter to external magnetic fields is desirable.
An important characteristic of the QCD vacuum is its transformation property under parity (P) and charge conjugation (C). In the absence of a θ-parameter, the theory prohibits violation of both the P-and CP-symmetries. Indeed, experimental bounds -mostly coming from measurements of the electric dipole moment of the neutron -on the degree of this violation turn out to be extremely tiny. Nevertheless, CP-violation could still be realized in the local sense, through fluctuations of CP-odd observables. One manifestation of this in the QGP phase created in heavy-ion collisions might be through the presence of domains with a non-trivial topological structure of the gluon fields (see, e.g., ref. [8]). Such a nonzero topology is indicated by the non-vanishing value of the topological charge Q top (defined below) within that particular domain. Since the magnetic field is odd under CP transformation, it is natural to expect that it can be used to effectively probe the CP-odd domains of the quark-gluon plasma and, thus, the CP-violating fluctuations in the QCD vacuum.
A possible realization of the coupling between the strong magnetic field and the nontrivial topological structure of the QGP is the so-called chiral magnetic effect (CME) [9,10]. For close to massless quarks, helicity is an approximately conserved quantity, and in strong magnetic fields the quark spins tend to align themselves either parallel (for positive charges) or antiparallel (for negative charges) to the external field. Therefore, right-handed, positively charged quarks and left handed, negatively charged quarks will tend to have their momenta parallel to the direction of the magnetic field. In a domain of the quarkgluon plasma with nonzero topological charge density, there is an imbalance between the number of left-and right-handed quarks, due to the Atiyah-Singer index theorem. As a consequence, a net current of quarks can be produced (anti)parallel to the external magnetic field, or, equivalently, the domain in question will be electrically polarized in the direction of the magnetic field. An alternative formulation of the effect is in terms of a chiral chemical potential [10], which couples to the anomalous axial current and creates a chiral imbalance by preferring right-handed over left-handed quarks.
The effects of the electric polarization of the plasma domains may persist at later stages of the collision. After hadronization takes place, this can result in a preferential emission of charged particles above and below the reaction plane [11,12]. Indications for such a charge asymmetry were observed in the STAR experiment at RHIC [13,14] and in the ALICE experiment at the LHC [15]. However, to access observables related to the CME, certain parity-even experimental backgrounds have to be taken into account, which complicates the interpretation of the observed data. Thus, the exact meaning of these results is still debated, see, e.g., refs. [16][17][18][19]. For recent reviews on the subject see, e.g., refs. [20,21].
The CME and topology-induced CP-violation have been studied in various approaches, ranging from effective theory/model calculations to Euclidean lattice simulations. The former include among others settings like the Nambu-Jona-Lasinio model with an additional coupling to the Polyakov loop (PNJL model) [22], the holographic approach [23,24], hydrodynamics (see, e.g., refs. [25,26]) or using a chiral effective action [27]. On the lattice, the CME was first studied by measuring current-and chirality fluctuations in quenched SU(2) [28] and quenched SU(3) gauge theory [29]. Surprisingly, around the transition JHEP04(2014)129 temperature, the fluctuations of the current parallel to the magnetic field were found to decrease with growing B in the small magnetic field region [28], a result which still lacks a qualitative understanding. Another approach to investigate the CME on the lattice is using the chiral chemical potential, see refs. [30,31]. Finally, the interplay between magnetic fields and topology was also studied by discretizing a continuum instanton configuration, and measuring the electric polarization in the presence of the magnetic field [32], see the illustration in the left panel of figure 1.
In the present paper, we pursue a different approach and measure the extent to which the topological charge and the electric polarization of the quarks correlate locally, when exposed to external magnetic fields. Instead of having to consider classical instanton configurations, this approach enables us to use real QCD gauge backgrounds and to consider the local fluctuations of the topological charge on them, see the right panel of figure 1. Moreover, while there is no need to introduce any anomalous current or chemical potential, the method still gives a handle on relating the topological and the electromagnetic properties of the QCD vacuum in a Lorentz invariant manner. This approach is similar to that of ref. [33], where chirality-electric polarization correlators were measured in quenched SU (2) gauge theory to detect the induced electric dipole moment of valence quarks. However, in our case the quarks and the external magnetic field are introduced dynamically, which allows us to observe spatially extended electric dipole structures in the QCD vacuum.
We indeed find that in local domains with nonzero topological charge density, an electric dipole moment is induced parallel to the external field. The strength of this effect is determined for various magnetic fields and temperatures around T c for several different lattice spacings. A scheme for defining the continuum limit of the results, utilizing the gradient flow [34] of the fields, is also introduced and used to perform the continuum extrapolation. Finally we compare the lattice results to a model calculation that employs nearly massless quarks and constant (anti)selfdual gluon backgrounds -a setting in which the problem can be treated analytically [35]. This comparison reveals that the numerical result found for full non-perturbative QCD with physical quark masses is by an order of magnitude smaller than the model prediction.

JHEP04(2014)129 2 Formulation
In our setup, we consider the local correlation of the quark electric dipole moment with the topological charge density, where G µν (x) is the SU(3) field strength at the point x, and tr denotes the trace in color space. The space-time integral of q top gives the total topological charge Q top . In order to define the local electric dipole moment operator, let us consider the spin polarization of the quark of flavor f (represented by the field ψ f ), where γ µ are the Euclidean Dirac matrices. In the presence of a constant Abelian external field F µν , the spin polarization develops a nonzero expectation value [36], where q f is the charge of the quark of flavor f , and the factor of proportionality τ f is conventionally written as the product τ f = ψ f ψ f χ f of the quark condensate and the magnetic susceptibility. Note that the expectation value in eq. (2.3) involves an integral over space-time and a normalization by the four-volume, to exploit translational invariance. The xy component of eq. (2.3) is induced by an external magnetic field F xy = B z , whereas the zt component by an external (Euclidean) electric field F zt = E z . Accordingly, the polarizations correspond to a magnetic and an electric dipole moment of the quark, respectively. 1 The electric dipole moment is a parity-odd quantity, just as the topological charge density of eq. (2.1); their product is therefore parity-even, and can have a nonzero expectation value in the presence of the external magnetic field. We consider this product locally and write it, similarly to eq. (2.3), as where we factored out the magnitude of the topological fluctuations to define the correlator of the two quantities. A similar combination was studied in ref. [33]. Since q top = 0 we use the square root of the expectation value of q 2 top for the normalization. In this way, similarly to eq. (2.3), we obtain an observable with mass dimension 3, and introduce the proportionality factorτ f . We emphasize that we consider the product of the two densities on the left hand side locally, in order to see the local correlation of topology and the polarization (this local product contains contact terms which need to be removed by an JHEP04(2014)129 adequate renormalization prescription, see section 3 below). Eq. (2.4) expresses the fact that there is a local correlation between the topological charge density of the non-Abelian vacuum and the induced electric dipole moment, and that this correlation is proportional to the external magnetic field. 2 Consider now the ratio of eq. (2.4) and the xy component of eq. (2.3). Here the external field cancels to leading order, giving directly the ratioτ f /τ f , which has dimension zero, and is particularly suited for the lattice determination. Note that in this ratio all multiplicative renormalization factors cancel.

Observables and renormalization
We calculate the expectation values appearing in eq. (2.5) on the lattice with an external magnetic field in the positive z direction, B z ≡ B. The lattice geometry is N 3 s × N t , and the lattice spacing is denoted by a, such that the spatial volume of the system is given by V ≡ (aN s ) 3 and the temperature by T = (aN t ) −1 . We consider the three lightest quark flavors u, d and s, for which the charges are q u /2 = −q d = −q s = e/3 (here e > 0 is the elementary charge). We derive our observables from the QCD partition function, which, in the staggered discretization of the fermionic action reads where β = 6/g 2 is the inverse gauge coupling, S g the gauge action and M f = M f (U, q f B, m f ) = / D(U, q f B) + m f 1 the fermion matrix, for which we apply two steps of stout smearing on the gluonic links U . The quark masses are tuned along the line of constant physics (LCP) as m u = m d < m s , ensuring that the isospin averaged zerotemperature hadron masses equal their experimental values [38] (for the present action the most precise LCP can be read off from figure 1 of ref. [39]). Further details of the action and the simulation setup can be found in refs. [38,40,41]. Since the external field couples directly only to quarks, B enters only through the fermion determinants. Note that the dependence on B is always of the form of the renormalization group invariant combination For the gauge action S g , we use the tree-level improved Symanzik action, which contains the product of links along closed loops of size 1 × 1 (the plaquettes P µν ) and of size 2 × 1. The topological charge (2.1) at the space-time point x can be calculated via the field strength G µν (x), which can be discretized as the sum of the antihermitian part of the four JHEP04(2014)129 plaquettes touching the site x, and the product in the four plaquettes starts at point x and advances counter-clockwise. To suppress the noise originating from short-range fluctuations, the links used in eq. (3.2) are the twice stout smeared links that we also use in fermionic observables. We find that this choice for the definition of G µν -and, thus, of q top -reduces the noise in the correlation between q top and the electric dipole moment, necessary for the coefficient C f of eq. (2.5).
Note that the continuum limit of C f is unaffected by this choice. Let us add here that it is customary to use improved definitions of q top (see, e.g., ref. [42]) or much more extensive smearing of the gluonic links in order to obtain an integer value for the total topological charge Q top . Here we do not aim to determine the total charge, or its susceptibility, but concentrate on local fluctuations in q top and its correlation with fluctuations of the electric dipole moment, for which we carefully checked that our setup is appropriate. The expectation value of the spin polarization with respect to the partition function (3.1) reads where the trace (in color and coordinate space) is determined using noisy estimators η i , such that the polarization at point x is (color indices are suppressed here) with no summation over x. Here, N v is the number of estimators, which we set, depending on the ensemble, in the range 40 . . . 80. Furthermore, σ µν stands for the staggered representation of the tensor operator, see ref. [43] for the implementation we use. Using the expressions (2.1), (3.2), (3.3) and (3.4), the expectation values appearing in the ratio C f of eq. (2.5) are determined. The so obtained C f is yet to be renormalized, since both its numerator and denominator contain divergent terms. The magnetic dipole moment, for example, contains a logarithmic additive divergence, which may be eliminated using the operator m f ∂/∂m f , see ref. [43]. The square of q top is also subject to renormalization, as it contains the contact term, see, e.g., ref. [44]. Similarly, one expects the numerator to contain terms that are infinite in the continuum limit. These divergences are related to the fact that two densities are multiplied at the same space-time point. To remove these unphysical contributions, we use the gradient flow [34] for the fields contained in q top and in Σ f µν . The gradient flow was shown to eliminate additive divergences in fermionic observables like the condensate or the pseudoscalar correlator [45]. Likewise, we find that evolving the fields up to a fixed physical flow time t ph -or, equivalently, applying a nonzero smearing range R s = √ 8 t ph -renormalizes the observable C f and, at the same time, suppresses noise considerably. Our implementation of the gradient flow is detailed in appendix A. Finally, the operator Σ f µν is also subject to multiplicative renormalization by the tensor renormalization constant Z T , which was calculated in perturbation theory for the present action in ref. [43]. However, this factor cancels in the ratio C f . Altogether, C f is ultraviolet finite, if the continuum limit is approached along a fixed nonzero smearing range R s . On the lattice, this corresponds to taking the limit N t → ∞ at a fixed temperature T = (N t a) −1 and tuning the smearing range in lattice units as R lat s = R s /a. We repeat the continuum extrapolation for several ranges R s > 0 and subsequently extrapolate the results to R s = 0.

JHEP04(2014)129
Let us point out that in the present study smearing is applied in two different contexts. First, stout link smearing is employed in the fermionic action in order to suppress lattice discretization errors and, thus, to improve the convergence towards the continuum limit. Second, the fields in certain observables are evolved according to the gradient flow, which is equivalent to performing infinitesimal smearing steps. The latter reduces unphysical ultraviolet contributions in some observables, allowing for a clean definition of the continuum limit.

Results
We first analyze the response of q top (x) · Σ f zt (x) to the external magnetic field. Together with the results for Σ f xy (x) , this is plotted for the down quark in the upper left panel of figure 2. The ratio of the two expectation values is expected to be independent of the magnetic field, up to corrections of O((q f B) 2 ), in accordance with Lorentz invariance. Within the range of the applied magnetic fields, these corrections are found to be small, and thus the ratio is to a good approximation constant, see the lower left panel of figure 2. In order to determine the leading order B-dependence of the ratio, in the following we fit the data either to a constant, or consider corrections of O((q f B) 2 ). Our strategy for the determination of the systematic error of the result will be discussed below. The next step to obtain the coefficient C f of eq. (2.5) is to measure the local fluctuations 3 in q top . We find that q 2 top (x) depends quadratically on eB (again in accordance with Lorentz invariance), however, with a coefficient that changes sign as the temperature is increased across the transition temperature T c . This behavior is reminiscent of that of the chiral condensate [41,46,47] as well as of the gluonic action [48], which undergo magnetic catalysis at low temperatures and inverse catalysis in the transition region. The change in the local fluctuations due to the magnetic field, ∆ q 2

JHEP04(2014)129
is shown in the right panel of figure 2 for different temperatures. We note that although this change is significant, its magnitude is negligible compared to q 2 top (x) at B = 0 for the magnetic fields under study, in accordance with the findings for the two-point function of the topological charge density in ref. [48].
We proceed with the renormalization, and investigate the effect of the gradient flow on the coefficient C f . According to our expectations, C f is unphysical for a → 0 at vanishing flow time (vanishing smearing range), whereas for any nonzero t ph ∝ R 2 s , it has a finite continuum limit. We demonstrate this in figure 3, where C u (R 2 s ) is shown for four lattice spacings at a fixed temperature T = 113 MeV. While a power-type divergence is clearly absent from the R s = 0 data points, a logarithmic divergence cannot be excluded. At finite smearing ranges, we observe the convergence of the results to improve drastically -at R 2 s ≈ 0.5 fm 2 , the data points for all lattice spacings lie essentially on top of each other. Moreover, we also observe that the signal to noise ratio improves by up to an order of magnitude as the smearing range is increased beyond 1 fm.
For each R s > 0 dataset, we extrapolate the results to the continuum limit by a JHEP04(2014)129 113 MeV 0.132(10) 0.130 (14) 0.096 (7) 163 MeV 0.14(2) 0.12 (3) 0.09(2) quadratic fit in the lattice spacing (motivated by the O(a 2 ) scaling properties of the action we use). For this extrapolation we use the three finest lattices and only include N t = 6 in the fit to estimate the systematic error. We find that the so obtained extrapolations are very well described by a linear function in R 2 s (i.e., linear in the physical flow time t ph ), which we use to extrapolate to R s = 0, see the left panel of figure 3 for the results for the up quark at T = 113 MeV. We also consider a quadratic dependence on R 2 s , which we do not find to improve the fit qualities. Altogether, we take into account 2 × 3 × 2 × 2 different fits (constant or quadratic fit in eB; including or excluding the point with the largest or the smallest eB; continuum extrapolation including or excluding N t = 6; linear or quadratic extrapolation in R 2 s to R s = 0). The a → 0, R s → 0 limits are used to build a weighted histogram, and the average value and systematic error is estimated -following refs. [49,50] -by the mean and width of the obtained distribution, respectively. ( figure 3 shows one representative fit out of the many.) The central values and the total (systematic and statistical) errors obtained from this procedure are given in table 1 and indicated by the gray regions at R s = 0 in figure 3.
We perform a similar analysis in the deconfined phase, at T = 163 MeV using three ensembles with N t = 6, 8 and 10. The coefficient τ f of the magnetic dipole moment quickly approaches zero as the temperature is increased, see ref. [43]. At the same time, the coefficientτ f of the topological charge density-electric dipole moment correlator is also found to drop, which lowers the signal-to-noise ratio in C f . Moreover, we also observe that the continuum extrapolated data at R s > 0 show a much less pronounced dependence on R 2 s , as compared to the case at T = 113 MeV, see the right panel of figure 3. Motivated by this, in addition to the linear fits we also fit the data to a constant to extrapolate to R s = 0. The systematic error is again found by considering the width of the histogram built from results obtained by the various fit procedures.
For the down quark -again as a consequence of the q f B-independence to leading order -the results are within errors consistent with those obtained for the up quark. We find C s to be somewhat suppressed compared to the light quark coefficients, due to the larger mass of the strange quark. Our final results in the continuum limit at R s = 0, for the two temperatures under consideration, are shown in table 1. Note that the values for the two temperatures agree within errors for all flavors. Finally we remark that within our range of magnetic fields (eB < 0.5 GeV 2 ), the behavior shown in the left panel of figure 2 persists also at nonzero smearing ranges R s > 0 in the gradient flow, and the ratio of polarizations q top (x) · Σ zt (x) / Σ xy (x) shows no significant dependence on B.
Interpreting Σ f µν as the electric dipole moment of the quark, it might seem that the

JHEP04(2014)129
induced polarization is point-like and is not related to spatial charge separation. However, due to the fluctuations in q top (x) and their interaction with dynamical sea quarks, the local electric dipole moment correlates with spatially extended dipole structures and, thus, with the spatial separation of the electric charge. To show that these extended structures exist, let us consider the electric current operator and compose the observable 4 where we employed the same normalization as in the definition (2.5) of C f . The ratio D f (∆) represents the correlation between the topological charge density and the electric charge density at two distinct points separated by a four-vector ∆. We remark that in our Euclidean setting, the correlator in the numerator of eq. (4.2) is imaginary. Since the observable contains no dependence on the (imaginary) time, its analytic continuation simply amounts to a multiplication by i, giving a real observable in Minkowski space-time. 5 In the left panel of figure 4 we show this correlator for the up quark in the xz plane. The figure reveals an excess of positive charge above (∆ z > 0) the topological 'source' and an excess of negative charge below (∆ z < 0) it. Thus, we indeed observe an electric dipole structure aligned with the magnetic field.
To show that this spatially separated electric charge is not a lattice artefact, in the right panel of figure 4 we plot D u (∆) for ∆ = (0, 0, ∆ z , 0). To approach the continuum limit in a well-defined manner, we again make use of the gradient flow and consider a nonzero smearing range. The results using three lattice spacings N t = 6, 8 and 10 lie almost perfectly on top of each other, showing small discretization errors and a fast scaling towards a → 0 -similarly as we observed for C f , compare figure 3. For the strange quark (which exhibits a better signal to noise ratio) we also considered the dependence of the spatial integral d 3 ∆ D s (∆) ∆ z -corresponding to the electric dipole moment of the configuration -on the smearing range. The results indicate that this integral remains nonzero even in the limit R s → 0. Altogether, we conclude that the spatial separation of the electric charge remains a well-defined concept in the continuum limit. 4 Note that in accordance with Lorentz-symmetry -namely that D f should be antisymmetric in the two indices appearing in its definition -the correlator involving J f t along ∆z equals minus the correlator involving J f z along ∆t. 5 To see this, note that the density of topological charge qtop(x) receives a factor of i via the continuation.

Comparison to a model
Let us now interpret our result for C f in a model and in the context of heavy-ion collisions.
It is instructive to think of the quark-gluon plasma as depicted in the right panel of figure 1, with small independent domains containing gluon backgrounds of topological charge Q top . In each domain an electric polarization Σ zt is induced by the magnetic field and by the local Q top . This can be compared to the magnetic polarization Σ xy , which is uniform in the whole volume. Let us further assume that the topological charge in each domain is created by constant selfdual or antiselfdual non-Abelian fields of strength G, such that Q top ∼ ±G 2 . This is a generalization of the approach in ref. [35] that allows to describe the case B > G as well as B < G. Like in ref. [35], both polarizations in the local domain can be calculated analytically, when other gluonic interactions are neglected. The calculation (for technical details, see appendix B) simplifies tremendously in the lowest-Landau-level (LLL) approximation, which amounts to neglecting the quark masses compared to the field strengths. The magnetic field-dependence of the polarizations then reads LLL, m 2 G, |B − G|, |B + G|: 1) where an identical proportionality factor has been neglected. This result is valid for one quark flavor (whose electric charge is set to unity) and gauge group SU(2) for spatially aligned Abelian and non-Abelian fields, and agrees with the calculation of ref. [35]. (Note that eq. (5.1) does not hold for G = 0 where the LLL approximation is invalid. In fact, in this limit Σ zt vanishes but Σ xy remains finite.) In appendix B we also discuss the case of non-aligned fields and gauge group SU(3) resulting in similar formulae.
Based on our lattice results, we can make two important statements about this model. First, the ratio q top (x) · Σ zt (x) / Σ xy (x) is found to be B-independent for QCD with physical quark masses in the relevant range of magnetic fields, cf. the left panel of figure 2.

JHEP04(2014)129
The equivalent of this quantity in one domain in the model treatment is Q top · Σ zt /Σ xy , which is independent of B -and, thus, reproduces the lattice findings -only if the non-Abelian scale G exceeds the external field B. Second, in this regime (B < G), we may compare the model prediction to the lattice results quantitatively. In order to compute the coefficient C f in the model, we need to assume a distribution of the topological charge among the local domains. A reasonable approximation is a Gaussian average 6 over Q top . In addition, we also consider an arbitrary angle ϑ between the non-Abelian and Abelian fields and integrate over ϑ. This averaging over Q top and over ϑ is denoted by . . . . Using eq. (5.1) for B < G and its generalization to non-aligned fields, eq. (B.24) -which we derived for B G -we obtain is a Q top -independent factor describing the dependence of the polarizations on the angle. Its average over ϑ cancels in the ratio C f . In fact, the above obtained number is independent of the width of the Gaussian distribution of Q top (due to the matching powers of Q top in the numerator and the denominator). However, it differs from our lattice determination C f ∼ 0.1 by an order of magnitude. Put differently, the strong interaction between quarks prevents their full polarization predicted by such a model.
The above comparison reveals that an effective description of QCD with magnetic fields has to take the strong interaction into account non-perturbatively and beyond the simple assumptions of this model. In the same spirit one can question the lowest-Landau-level approximation used in the model setting. It corresponds to the idealized situation where the quark mass vanishes, and all quarks which are spin polarized by the magnetic field interact with the gluonic background and contribute to the electric polarization. However, heavier quarks are less sensitive to topology, and, accordingly, we expect the ratio Σ zt /Σ xy to decrease as m grows. This is consistent with our results C u,d > C s .

JHEP04(2014)129
We also estimated this coefficient using a model calculation employing nearly massless quarks, the lowest-Landau-level approximation and constant selfdual gluon backgrounds. This model was found to overestimate C f by an order of magnitude. In other words, there is a substantial quantitative difference of the strength of local CP-violation for quasifree quarks used in model approaches and fully interacting quarks in realistic physical situations. Whether the electric current in the formulation of the CME with a chiral chemical potential [10] is also subject to a similar suppression due to non-perturbative QCD effects (first lattice results indicate a suppression by a factor of 3-4 [30]), does not follow directly from our results. However, we take the results as a hint that effects due to local CP-violation in general contain similar suppression factors.
Let us finally add that we employed the staggered discretization of the QCD quark action in the lattice simulation, which in some topology-related aspects gives rise to large systematic/discretization errors. The topological susceptibility, Q 2 top /V , for example, shows a rather slow scaling towards the continuum limit, see, e.g., ref. [52]. We find that for our particular observable, C f , the continuum extrapolation is much flatter. This may have to do with the fact that C f is a local observable whereas the susceptibility is not. Nevertheless, it would be desirable to confirm our numerical findings with chiral fermions that have nicer topological properties.

A Details of the gradient flow
The smearing of the gluonic and fermionic fields is performed by evolving these fields in flow time t (t ph = t · a 2 is the physical flow time). The evolution in flow time amounts to finding the solution of the flow equations for the gluonic links [34], and for the quark fields [45], and the corresponding equation forψ t f . Here, Z(U t ) is the (algebra-valued) derivative of the plaquette action with respect to the link variable, and ∆ is the lattice discretization of the Laplace operator (see below). The solution of the flow equations can be found by numerical integration, which is done using the third-order Runge-Kutta integrator described in refs. [34,45] (a stepsize of ∆t = 0.02 was found to be optimal here, see also ref. [53]). Integrating the flow equations up to a fixed physical time t ph = t · a 2 corresponds to a smearing of the fields over a range of R s = √ 8 t ph [34].

JHEP04(2014)129
The definition of the quark condensate -or, of the fermionic bilinears appearing in eq. (2.5) -at nonzero flow time requires the use of the adjoint flow for the noisy estimators η i of eq. (3.4) from flow time t back to flow time 0, see ref. [45]. For this, an optimal scheme for the storage of the evolved links U t for 0 ≤ t < t is implemented. The evolution along the gradient flow is started from our original gauge action, thus with unsmeared links. The stout smearing is then applied only for the measurement of the operators, see eqs. (3.2)-(3.4).
We remark that there is a peculiar issue that arises if one applies the fermionic gradient flow for staggered quarks in a naive way. In the staggered fermionic discretization, the Dirac components of the quark field ψ at site x are distributed over vertices of the fourdimensional hypercube touching x. This distribution of the components is devised in a manner such that the staggered action becomes diagonal in Dirac space, and the only remnant of the original Dirac structure is through space-dependent real numbers, the so-called staggered phases. In particular, the mass termψψ and the Dirac operatorψ / Dψ are diagonal in Dirac space, therefore they can be represented in terms of the staggered quark fields χ in the same form, e.g.χχ. However, the naive discretization of the Laplace operator is not diagonal after the staggered transformation, giving no straightforward correspondence between the representation with the original fields and that with the staggered fields.
To construct the Laplacian, let us define the forward and backward covariant difference operators, where U µ ∈ SU(3) are the gluonic links and u µ ∈ U(1) the phases corresponding to the magnetic field. The naive one-step discretization of the Laplace operator, ∆ naive = ∇ † µ ∇ µ indeed becomes off-diagonal as it mixes the tastes distributed over the hypercube in a non-trivial way. One possibility to avoid this mixing of the tastes is to use the two-step discretization of the covariant differences, µ . This two-step discretization was used in the flow equation eq. (A.2). The non-diagonal nature of ∆ naive results in an explicit Lorentzsymmetry breaking of the evolved fermionic fields, even at B = 0. This is indicated by asymmetric expectation values of the bilinear structuresψγ µ ψ. Using the two-step Laplacian ∆ diag , (the lattice discretized version of) Lorentz-symmetry is maintained, and ψ γ µ ψ = 0 for all µ.
Finally we remark that we also attempted to use the square of the staggered Dirac operator in place of the Laplacian for the evolution of the fermionic fields in eq. (A.2). The results obtained for the coefficient C f after the flow with / D 2 , however, showed an inferior

JHEP04(2014)129
scaling towards a → 0, as compared to the case with ∆ diag . Performing the extrapolation to the continuum limit was only feasible for the latter choice.

B Polarizations in topological backgrounds
In order to evaluate C f in a topological background, we consider one quark flavor in constant commuting selfdual or antiselfdual non-Abelian fields, which exist in a finite Euclidean box with quantized fluxes [54] (they can also be thought of as fields deep inside instantons or antiinstantons [35]) plus an Abelian magnetic field B. In these backgrounds, both the topological charge Q top and the polarizations Σ zt and Σ xy are constant in space. The quark mass is denoted by m and the electric charge is set to unity for simplicity. Moreover, our notation is such that the QCD coupling does not enter the covariant derivative. We follow two equivalent approaches to determine the electric and magnetic polarizations in this model setting. First we employ a spectral representation of the observables using Landau-levels. Second we write down the polarizations using the exact quark propagator in the specific background.

B.1 Polarizations using the spectral representation
Let us first consider the case where the non-Abelian field G is (anti)parallel to the Abelian one B. Without loss of generality we can assume that B points in the z direction. Taking SU(2) for the non-Abelian group, the xy and zt components of the total field strength f read where we diagonalized the field strengths via a gauge transformation (for constant field strengths this is always possible). We also inserted the sign of the topological charge Q top = ±G 2 /(2π 2 ) in the electric components to account for both the selfdual and the anti-selfdual cases. Let us first discuss the upper color component and denote b ≡ B + G and e ≡ sign(Q top )G. The Dirac eigenvalues of this system are obtained through two independent Landau-level problems in the (x, y)-and (z, t)-planes (the arrows indicate the eigenvalues of the corresponding operators),

JHEP04(2014)129
where |b||e|/(4π 2 ) is the degeneracy of all Landau-levels. The spin-dependence is such that only the corresponding lowest Landau-levels contribute: {n b = 0, s b = −sign(b)} for Σ xy , whereas {n e = 0, s e = −sign(e)} for Σ zt , cf. appendix B in ref. [43], giving Note that the polarizations change sign when their corresponding field strengths e or b are reversed, as they should. The sum in h contains an m-and |f |-independent divergence, arriving at eq. (5.1) used in section 5. The first and third lines agree with eqs. (81) and (82) of ref. [35], while the second line can also be obtained from the number of zero modes, eq. (47) of that reference. Note that at |G| = B, where Σ zt would have a cusp, the lowest-Landau-level approximation breaks down in the color sector with field strength −|G| + B.
We now turn to the gauge group SU(3). One can again diagonalize the field strength, now it has two independent amplitudes in the fields (G 1 , G 2 , −G 1 − G 2 ) in three color sectors, and |Q top | ∼ [G 2 1 + G 2 2 + (G 1 + G 2 ) 2 ]. This slightly complicates the calculations. For the simplest case of space-parallel fields in the lowest Landau-level approximation one gets, in analogy to (B.6)-(B.7)

JHEP04(2014)129
We have found that the ratio Σ zt /Σ xy is B-independent and equals sign(Q top ) when all three fields |G 1 |, |G 2 | and |G 1 + G 2 | are large compared to B, and that it is smaller and becomes B-dependent if one of them is not.

B.2 Polarizations using the exact propagator
We proceed by generalizing the above calculation and allow for an arbitrary polar angle ϑ between the non-Abelian and Abelian fields, This case should be equally relevant for estimating C f in realistic QCD configurations. We again considered the selfdual and the antiselfdual cases simultaneously by inserting sign(Q top ) in the electric fields.
It is now advantageous to represent the polarizations (first we discuss a single color sector) by

JHEP04(2014)129
Let us denote the invariants of f (proportional to 'action' and 'topological charge' density) as Then the eigenvalues are given by [55] .
(B.15) By explicit comparison we found that the other factor appearing in S(x, x)σ αβ can be represented as The second line confirms that Σ zt changes sign if the topological charge does so. Plugging all this into the proper time integral (B.11) shows that the integral diverges as t → 0. This is the same divergence that we encountered in eq. (B.5). Here we eliminate it by dividing the observable by m and differentiating it with respect to m 2 , cf. eq. (B.4). This indeed renders the integral finite and also reveals that the divergence is independent JHEP04(2014)129 of m and of the fields. Setting ϑ = 0 (and consequently z = ±1 etc.) reproduces the finite part of eq. (B.5).
For ϑ = 0, two hyperbolic sine functions are left in the denominator of eq. (B.15), such that the proper time integral cannot be performed easily. Since we argued that the region |G| > B (i.e. |G| > B cos ϑ for all ϑ) is the relevant one for comparison with the lattice data, we now specialize to this case. Similarly as in eq. (B.6), we now resort to the approximation m 2 |G|, B. Moreover, in order to simplify the integral, we also assume B |G|. Expanding the fraction in B/|G| and m 2 /|G|, we can perform the t-integral to arrive at Here the leading term ∼ |G|G vanishes upon adding the second color sector of SU(2), which amounts to the same expression with G → −G, cf. eq. (B.7). Adding the contributions of both sectors and integrating in m 2 we finally get (−4π 2 m) · Σ xy = 2B|G| · g(ϑ), (−4π 2 m) · Σ zt = sign(Q top ) · 2B|G| · g(ϑ), (B.24) which, at ϑ = 0, reproduces eq. (B.7) for the case B < G. This expression was inserted in eq. (5.2). Note that the average over the polar angle factorizes and gives