Numerical Study of the Chiral Separation Effect in Two-Color QCD at Finite Density

We study the Chiral Separation Effect (CSE) in finite-density SU(2) lattice gauge theory with dynamical quarks. We find that the CSE is well described by the free quark result in the high-temperature quark-gluon plasma phase. As one enters the confinement regime with broken chiral symmetry at chemical potential smaller than half of the pion mass, the CSE response is gradually suppressed towards low temperatures in comparison to the free quark result. This suppression can be approximately described by assuming that the CSE current is proportional to the charge density, rather than the chemical potential, as suggested in the literature (ArXiv:1712.01256, Phys.Rev.D97(2018)085020). We also provide an upper bound on the contribution of disconnected fermionic diagrams to the CSE, which is consistent with zero within our statistical errors and small compared to that of the connected diagrams. Our results are obtained mainly in the QCD-like regime of SU(2) gauge theory at low densities, and hence should be at least qualitatively applicable to QCD as well.


I. INTRODUCTION
Anomalous transport phenomena are transport responses of quantum matter that originate in quantum anomalies, inevitable violations of classical symmetries upon quantization [1,2]. In particular, in strongly interacting matter described by Quantum Chromodynamics (QCD) the classical symmetry between left-handed and right-handed fermions is violated by the Adler-Bell-Jackiw axial anomaly [3]. This violation manifests itself in the infamous Chiral Magnetic Effect [4] -the generation of an electric current along a magnetic field in chirally imbalanced matter -as well as the closely related Chiral Separation Effect (CSE) [5,6] -the generation of an axial current along an external magnetic field in a dense medium (see Fig. 1).
In the last decade, anomalous transport phenomena in QCD matter were systematically and intensely studied in heavy-ion collision experiments at the RHIC [7] and LHC [8] colliders, and will also be studied at the NICA collider [9]. These studies are not conclusive yet due to large background effects, which contaminate the signatures of anomalous transport [8,[10][11][12][13][14]. A dedicated run with isobar nuclei has been recently completed at RHIC in order to disentangle these background effects [7,15], and the produced experimental data is currently being analyzed [16]. Just as the viscosity of the quark-gluon plasma is related to hadronic elliptic flow [17], anomalous transport coefficients characterizing the strengths of the Chiral Magnetic and Chiral Separation Effects can be related to correlations of angular distributions of oppositely charged hadrons in heavy-ion collisions [18,19].
One of the most popular ways to interpret experimental data on these correlations relies on the anomalous-viscous fluid dynamics (AVFD) framework [20][21][22]. AVFD is based on anomalous hydrodynamics [23][24][25] which incorporates anomalous transport along with more conventional transport responses such as viscosity and electric conductivity. Hydrodynamic simulation codes which include anomalous transport phenomena on an event-by-event basis are currently being ac-tively developed [21,[26][27][28] and are becoming more and more realistic.
The anomalous hydrodynamic description of QCD matter requires the values of anomalous transport coefficients as an input. For a single-component chiral fluid, anomalous transport coefficients are fixed by thermodynamic consistency [24,29]. On the other hand, in a quark-gluon plasma (nearly) chiral quarks interact with dynamical non-Abelian gauge fields, which themselves behave as a viscous fluid and in fact dominate the hydrodynamic flow. When interactions with dynamical gauge fields are present, all anomalous transport coefficients might receive both perturbative [30] (see Fig. 1) as well as non-perturbative [6, [31][32][33][34][35] corrections. However, at present not much is known about the magnitude of these corrections.
In a few lattice gauge theory simulations, the Chiral Magnetic [36] and Chiral Vortical [37] Effects were studied by measuring the responses of naively discretized axial currents and energy momentum tensors to constant external magnetic or axial magnetic fields (the study [37] of the Chiral Vortical Effect used a trick to replace rotation by a background axial gauge field). In both works [36,37] the CME and CVE transport coefficients were found to be 5 to 20 times smaller than those obtained for free quarks, both in high-and in low-temperature phases. If the corrections in the full gauge theory should indeed have the effect to make the anomalous transport responses so small, they might well be unobservable in current heavy-ion collision experiments.
However, the works [36,37] used non-chiral lattice fermions with non-conserved currents and an energymomentum tensor without proper renormalization. On the other hand, a numerical study of the CSE in quenched SU (3) lattice gauge theory with exactly chiral overlap fermions and a properly defined axial current [38] found no noticeable corrections to the free quark result in both confinement and deconfinement phases. For free quarks, the axial current induced by the CSE is where j A i = fq f γ 5 γ i q f is the axial current density with quark fieldsq f , q f of flavor f and N c colors, µ is the quark chemical potential, B i is the magnetic field, and C em = f Q f = Tr (Q) is the electromagnetic charge factor in which Q f denotes the electric charge of quark flavor f . The result (1) corresponds to the triangular diagram in Fig. 1 a). To simplify notation, in what follows we assume that C em , which appears in all formulae as a simple prefactor, is equal to unity: C em = 1. The correct value of C em can be restored in all results in an obvious way.
In this paper we study the Chiral Separation Effect in the full gauge theory with dynamical quarks, taking into account the contributions of virtual fermion loops and disconnected fermionic diagrams like the one in Fig. 1 b). These contributions are expected to modify the free quark result (1) and are thus important to estimate the detectability of the Chiral Separation Effect in heavy-ion collisions. Rather than studying the theoretically clean, but rather academic, limit of exactly chiral quarks, we address the fate of the CSE in a more realistic setup with finite quark and pion masses and at finite temperatures in the vicinity of the chiral crossover. While giving some general insight into the magnitude of non-perturbative corrections to anomalous transport coefficients, this might also help to estimate the observable consequences of the CSE, such as the electric quadrupole moment of the quark-gluon plasma [39] due to Chiral Magnetic Waves [40].
Since the CSE is a feature of finite-density fermions, studying it in QCD would require simulations at finite baryon density, which are complicated by the infamous fermion-sign problem [41]. With current simulation methods one could only obtain first-principle lattice QCD results for µ/T ≪ 1.
In this work we circumvent the fermion sign problem by using two-color QCD, i.e. the SU (2) gauge theory with N f = 2 light quark flavours instead of QCD. The path integral weight is manifestly positive in this case, thus the sign problem is absent and the theory can be simulated at finite density. The SU (2) gauge theory is expected to be qualitatively similar to QCD at small densities µ < m π /2 [42,43]. In this regime there is a conventional QCD-like chiral crossover at some finite temperature, which separates the quark-gluon plasma regime and the hadron gas regime dominated by light pion states [44][45][46][47][48][49]. We expect that due to this qualitative similarity our numerical study of the CSE in SU (2) gauge theory is also at least qualitatively relevant for real QCD at small densities.
At larger densities, for µ > m π /2, SU (2) gauge theory is no longer similar to QCD, since the chiral condensate qq is rotated into the diquark condensate qq . Diquarks are bound states of two quarks which are color singlets and hence "bosonic baryons" in the SU (2) gauge theory. Instead of the first-order liquid-gas transition of nuclear matter, one therefore observes Bose-Einstein condensation together with a BEC-BCS crossover inside the diquark condensation phase [50,51]. Similarity to QCD, although at a different conceptual level, can be again expected at very large densities and low temperatures, in the conjectured quarkyonic and color-superconducting phases [52,53].

II. LINEAR RESPONSE APPROXIMATION FOR THE CHIRAL SEPARATION EFFECT
Within the linear response theory, the Chiral Separation Effect is characterized by the spatial correlator of vector and axial-vector currents j A i (k) j V j (−k) , where k is the space-like momentum. Assuming that the only nonzero momentum component is the spatial component k 3 [54], at small momenta this correlator can be written as where σ CSE is the anomalous transport coefficient in (1) characterizing the strength of the CSE. For free quarks, It is therefore convenient to define a momentumdependent transport coefficient σ CSE (k) as In the low-momentum hydrodynamic regime, the anomalous transport coefficient σ CSE in (1) is given by the zeromomentum limit of σ CSE (k). For exactly chiral fermions which may interact with each other, but not with other dynamical degrees of freedom (like dynamical gauge fields), σ CSE is expected to be universal and equal to the free fermion result due to the relation with the Adler-Bell-Jackiw axial anomaly [5]. However, corrections are still possible for nonzero quark mass [6] and due to fermionic disconnected diagrams like 1b) in Fig. 1. As calculations of [6] suggest, corrections to the CSE can be related to the amplitude g π 0 γγ for the π 0 → γγ decay: Within the linear sigma model where m is the constituent quark mass. In this case, σ CSE is still approximately linear in µ.
Another calculation of the flavor non-singlet CSE axial current j a A generated by a finite isospin chemical potential was carried out within chiral effective field theory in [35]. It suggests that in the low-temperature phase, where the CSE current is saturated by pions, σ CSE is proportional to the isospin charge density ρ a V rather than the isospin chemical potential, While this calculation is not directly applicable to the flavor-singlet axial current in (1) and (5), at least in the large-N c limit the flavor-singlet axial current should behave similarly to the flavor non-singlet one [55]. Our numerical results presented in Section IV below suggest that a parametrization similar to (6) might also work at finite density in the SU (2) gauge theory with dynamical quarks.

III. LATTICE SETUP
In this work we use the same lattice setup and the same ensembles of gauge field configurations as in our recent works [44,56], so here we will provide only a brief summary. We use the standard Hybrid Monte-Carlo algorithm with a tree-level improved Symanzik gauge action and N f = 2 flavours of mass-degenerate rooted staggered fermions with bare lattice quark mass am stag q = 5 · 10 −3 .
Our lattices have spatial sizes L s = 24 (am π L s = 3.8) and L s = 30 (am π L s = 4.7) and varying temporal extent L t = 4, 6, . . . , 22 to control the temperature T = 1/aL t . We use a single value of the lattice gauge coupling β = 1.7, thus working in a fixed-scale approach, which significantly simplifies the analysis of renormalization of lattice observables.
Most of our low-temperature ensembles with L t ≥ 12 were generated with a small diquark source λqq in the action with aλ = 5 · 10 −4 in order to facilitate diquark condensation, which would otherwise be impossible in a finite volume. This diquark source has very little effect on current-current correlators outside of the diquark condensation phase [44], see also Fig. 4. Estimates of phase boundaries based on our data sets are shown in Fig. 2 (see [44] for full details). We use Domain Wall (DW) and Wilson-Dirac (WD) valence fermions to measure the correlators of axial and vector currents in (2). On the one hand, for DW fermions the renormalization factor Z A for the flavour-singlet axial current is expected to deviate from unity by at most few percent [57], which is below our the statistical uncertainty of our Monte-Carlo simulations (and also well below experimental uncertainties). On the other hand, DW fermions are computationally very expensive, and we use the cheaper WD fermions to produce results with better precision covering more points on the phase diagram. A comparison between the results obtained with DW and WD fermions further demonstrates the smallness of axial current renormalization. We do not use staggered valence fermions in order to avoid artifacts related to the unphysical taste symmetry. Such a mixed lattice action with staggered sea fermions and DW valence fermions has already been used in a number of studies of the nucleon axial charge [57,58].
We tune the bare quark masses a m DW q = 0.01 and a m W D q = −0.21 in the DW/WD Dirac operators to match the pion mass am stag π = 0.158 ± 0.002 obtained with staggered valence quarks. The ratio of pion to rhomeson mass is m π /m ρ ≈ 0.4. To improve the chiral properties of DW and WD fermions without using much finer and larger lattices, we follow [58] and use HYP smearing [59] for gauge links in the DW and WD Dirac operators. For DW fermions the lattice size in the fifth dimension is L 5 = 16, which is typically sufficient to suppress additive mass renormalization [57,58].
For WD fermions, we use the conserved vector current where D xy is the Dirac operator, x, y, z, . . . label lattice sites, γ µ are the Euclidean gamma-matrices, with P ± µ = (1 ± γ µ )/2, U z,µ are the SU (2)-valued link variables,μ denotes the unit lattice vector in the direction µ, and θ z,µ is an external U (1) gauge field. We also use the conventional point-split definition of the axial current for WD fermions [60], j A z,µ x,y = = iγ µ γ 5 U z,µ δ x,z δ y,z+μ − iγ µ γ 5 U † z,µ δ x,z+μ δ y,z . (8) For DW fermions, the four-dimensional vector and axial currents are defined in the standard way by summing the five-dimensional conserved current over the fifth dimension. For the vector current a unit weight is used, for the axial current the summation weight changes from +1 to −1 in the middle of the lattice extending in fifth dimension [61]. The five-dimensional conserved current has a form similar to (7), except that the index µ takes five values and x, y, z live on the five-dimensional lattice with open boundary conditions along the fifth dimension.
We measure the contributions of both connected and disconnected fermionic diagrams to the axial-vector current-current correlator in (2). In coordinate space these contributions are: j A x,µ j V y,ν disc = Tr j A x,µ D −1 Tr j V y,ν D −1 , (10) where the traces are taken over the lattice site, spinor and color indices of the quark fieldsq, q. The disconnected contribution is measured using standard stochastic estimator techniques. After measuring j A x,µ j V y,ν conn and j A x,µ j V y,ν disc in coordinate space, we perform a discrete Fourier transform to obtain the momentum-space correlators which enter the linear response relations (2).

IV. NUMERICAL RESULTS
In Fig. 3 we present our lattice results for the momentum-dependent CSE transport coefficient σ CSE (k) defined in (4). For comparison, we combine the results obtained with DW and WD fermions, and with the spatial lattice sizes L s = 24 and L s = 30. For WD fermions on the L s = 24 lattices we show the contributions (9) and (10) of both connected and disconnected fermionic diagrams, for other data sets only the connected contributions are shown.
We also compare the gauge theory results with the results obtained for free WD quarks on same lattices. For the free WD quarks, we use a bare quark mass of am W D q = 0.01 (as compared to am W D q = −0.21 in the full gauge theory), since for free quarks there is obviously no mass renormalization. Therefore, in this case we choose the same bare quark mass as for the DW fermions, for which mass renormalization is expected to be weak.
In our calculations we also combine results obtained with zero diquark source λ at high temperatures (L t < 14) and with aλ = 5·10 −4 at low temperatures (L t ≥ 14). In Fig. 4 we demonstrate that for these two values of λ the CSE transport coefficients σ CSE (k) are practically indistinguishable.
We see from Fig. 3 that for most values of temperature and chemical potential the momentum-dependent CSE transport coefficient σ CSE (k) is very close to the corresponding free quark result. An explicit calculation of σ CSE (k) for free quarks at finite temperature in the continuum is sketched in Appendix A. The results of this continuum calculation are shown in all plots of Fig. 3 as solid black lines.
The gauge theory result for σ CSE (k) only becomes noticeably smaller than the free quark result at small values of the chemical potential aµ 0.10 and low temperatures L t 16 (see e.g. the plot for L t = 20 and aµ = 0.05 corresponding to µ = 0.32 m π in Fig. 3). In this regime SU (2) gauge theory is expected to be qualitatively similar to real QCD, thus the observed suppression of the CSE in the confined and chirally broken phase is also likely to happen in low-temperature, low-density QCD.
For aµ = 0.05 and aµ = 0.20 we have also calculated σ CSE (k → 0) for two different spatial lattice sizes, L s = 24 and L s = 30. Since the discrete momentum values for both lattices do not coincide, the low-momentum data for both lattices cannot be compared in a direct way. To overcome this difficulty, we construct a thirdorder spline interpolation of the data for both lattices, which is shown on Fig. 3 as solid lines going through the corresponding data points. Spline interpolations for both lattice sizes coincide with a very good precision for almost all available data sets. We only observe a noticeable difference between the interpolations of the L s = 24 and L s = 30 data on the plot for L t = 12 and aµ = 0.20, with L s = 30 data being much closer to the free quark results. Our data therefore suggest that the suppression of CSE that we observe at low temperatures and densi- ties is not a finite-volume artifact. In fact, the observable finite-volume effects are somewhat stronger in the high-temperature and high-density regime than at low temperatures. The contribution of disconnected fermionic diagrams is consistent with zero within our statistical errors for all values of chemical potential and temperature. The upper bound which we are able to set on these disconnected contributions appears to be least strict for low temperatures and small µ -that is, exactly in the corner of the phase diagram where also the connected contributions deviate most strongly from the free quark result (see Fig. 3, plot for L t = 20 and aµ = 0.05). We note, however, that the calculation of axial-vector current-current correlators is most difficult precisely in this regime, since because of the small chemical potential the CSE signal is also small compared to statistical fluctuations.
The values of σ CSE (k) calculated with Wilson-Dirac (WD) and Domain Wall (DW) fermions appear to be very close to each other. This suggests that the effect of multiplicative renormalization of the axial-current operator is small and plays a minor role in comparison with our statistical errors. Indeed, the renormalization factor for the axial singlet current typically appears to be close to unity on fine lattices with sufficiently light pions, especially for Domain Wall fermions [57,58]. For these reasons, we have not determined the precise value of Z A in this work. Small deviations of Z A from unity can also be expected to be significantly smaller than any systematic and statistical uncertainties in the experimental detection of anomalous transport phenomena.
Let us now investigate the CSE transport coefficient σ CSE (k → 0) in the low-momentum limit, which is most relevant for anomalous hydrodynamics [54]. To this end in Fig. 5 we illustrate the temperature dependence of σ CSE (k) at the smallest nonzero value of the lattice momentum, a k min = 2π/L s , and at different values of the chemical potential µ, and compare it with the corresponding results for free quarks on lattices of the same size and for the same momentum.
Again, we observe that σ CSE (k min ) becomes significantly smaller than the free quark result only at low temperatures L t 14 and small values of the chemical potential µ 0.10. This region of the phase diagram almost coincides with the QCD-like regime with spontaneously broken chiral symmetry, which should be dominated by pions. In this regime, the free quark result for σ CSE (k min ) has a rather weak temperature dependence, but the gauge theory result is strongly suppressed towards lower temperatures.
A comparison of the values of σ CSE (k min ) for L s = 24 and L s = 30 might make an impression that the suppression of σ CSE (k min ) becomes weaker for larger volumes. However, the apparent deviation of the results for L s = 24 and L s = 30 is caused simply by the difference of the minimal nonzero momenta a k min = 2π Ls for different lattice sizes. As one can see from Fig. 3, the data interpolated to continuum momentum values does not show significant volume dependence.
Let us now try to describe the observed CSE suppres-sion at low temperatures and densities in terms of some phenomenological formula. Let us try a formula of the form (6), but with a flavour-singlet current and chemical potential: To this end, in Fig. 5 we also show the rescaled charge density α ρ V (µ) /µ, tuning the coefficient α to achieve the best coincidence between the data points for σ CSE (k min ) /µ and α ρ V (µ) /µ. To achieve this we minimize the mean squared deviation of α ρ V (µ) /µ from σ CSE (k min ) /µ at a µ = 0.05 for the lattice size L s = 30.
We find that at low temperatures L t 14 the dependence of σ CSE (k min ) /µ for a µ = 0.05 on L t can be indeed well described by the formula (11) with α ≈ 14 a 2 (here we simply express α in lattice units, without implying that α scales as a 2 ). The same value of α also describes the data for a µ = 0.10 and L t 16 reasonably well.
For larger values of the chemical potential, a µ = 0.20 and a µ = 0.50, the formula (11) does not seem to work. Interestingly, Fig. 5 suggests that for these values of µ the CSE coefficient σ CSE (k min ) appear to be even slightly larger than the free quark result on the same lattice, although the deviations do not exceed 20%. According to Fig. 2, the data points for a µ = 0.20 and a µ = 0.50 are within the quark-gluon plasma regime or at the boundary of the diquark condensation phase, except for the point with aµ = 0.50 and L t = 20 that is in the diquark condensation regime. While one can see a drop of σ CSE (k min )/µ by around 20% for this single data point, this change is comparable to statistical errors, and we cannot make definite conclusions on the behavior of σ CSE (k min ) in the diquark condensation phase. We can only rule out a significant suppression of the CSE around the boundary of this regime.
According to the formula (6) with C em = Tr (Q) = 1, the coefficient α in (11) should be related to the pion decay constant f π via α = N c /(2πf π ) 2 . Using our estimate α ≈ 14 a 2 , from this relation we roughly estimate a f π ≈ 0.06, which has a reasonable order of magnitude when compared to the mass of the ρ meson, a m ρ ≈ a m π /0.4 ≈ 0.4 in our calculations. Thus fitting our data with the formula (11) implies the ratio f π /m ρ ≈ 0.15, as compared to f π /m ρ ≈ 0.12 with f π ≈ 93 MeV, m ρ ≈ 770 MeV in QCD.
In order to present further evidence for the scaling of σ CSE with ρ V , in Fig. 6 we show the dependence of the ratios σ CSE (k min ) /µ (on the left) and σ CSE (k min ) / a 2 ρ V (µ) (on the right) on the chemical potential µ at different temperatures with L t ≥ 16. It is quite obvious from these figures that the ratio σ CSE (k min ) / a 2 ρ V (µ) shows smaller relative deviations from a constant value α a −2 = 14 than the ratio σ CSE (k min ) /µ does from N c /(2π 2 ). In particular, for L s = 30 all data points contain the value α a −2 = 14 within their error bars. The data points for other ensembles deviate from this value by not more than 50%, FIG. 6. The ratio of the low-momentum CSE transport coefficient σCSE (kmin) and the chemical potential µ (left), compared to that of σCSE (kmin) and the charge density ρV (µ) (right). The solid black line in the left plot corresponds to the free continuum quark result σ 0 CSE (k → 0) /µ = Nc/(2π 2 ). In the plot on the right, the solid black line marks the optimal value α = 14 a 2 in our phenomenologically motivated formula (11).
whereas for the ratio σ CSE (k min ) /µ the relative deviation between different data points is much larger, up to a factor of 5. Note that the ratio σ CSE (k min ) / a 2 ρ V (µ) has larger error bars than σ CSE (k min ) /µ because ρ V (µ) also has statistical errors, in contrast to µ.
While these observations give some qualitative support to the formula (11), because of a conceptually different status of axial-singlet and axial-non-singlet currents in low-energy chiral effective theory, at the quantitative level one can expect further corrections to this formula.
A general conclusion that we can make based on our results is that the CSE becomes more and more suppressed for low temperatures in the confinement regime and for values of the chemical potential roughly smaller than half of the pion mass. At higher temperatures and densities, the CSE transport coefficient approaches its value for free quarks.

V. CONCLUSIONS
To summarize, we have found that in a gauge theory with dynamical quarks the Chiral Separation Effect is very close to the free quark result in the hightemperature quark-gluon plasma regime, and is gradually suppressed towards lower temperatures and densities in the low-temperature hadron resonance-gas regime with broken chiral symmetry at T T c and µ m π /2. Exactly this regime of SU (2) gauge theory is similar to the low-temperature, low-density phase of real finite-density QCD, thus our findings should be also relevant for real QCD at least qualitatively.
Note that our conclusions on the CSE suppression do not contradict the results of a previous study [38] in quenched SU (3) lattice gauge theory, where exactly chiral valence quarks were used, thus the pion mass was formally zero and the region with µ m π /2 was absent.
The suppression of the CSE at low densities and temperatures can be approximately described if one assumes that the CSE current is proportional to the charge density rather than chemical potential, as in equation (6). While the formula (6) was derived in [35] for the axial non-singlet current, we see that at the qualitative level it also applies to the axial singlet current, which has a different status within the chiral effective theory.
Contributions from disconnected fermionic diagrams to the Chiral Separation Effect appear to be consistent with zero within our statistical errors. The latter are relatively large at low temperatures and densities. For this reason we cannot rule out that the disconnected contribution might become important when the connected one is strongly suppressed. This scenario would certainly be interesting from a theoretical viewpoint. However, as the sum of both contributions still appears to be small when compared to the free massless fermion result, it is probably not very relevant for the experimental detection of anomalous transport phenomena.
One potential source of systematic errors in our study is the multiplicative renormalization of the axial current. However, comparison of the data obtained with Domain Wall and Wilson-Dirac fermions, as well as previous lattice studies of axial current renormalization, suggest that the renormalization effects are small.
Another potential improvement of our study would be the use of twisted boundary conditions in one of the spatial directions in order to access σ CSE (k) for nonzero momenta k that are smaller than k min = 2π Ls and perform a proper extrapolation of the numerical data for σ CSE (k) to the limit k → 0. Fig. 3 suggests that extrapolation from k min towards k = 0 can potentially make the estimates of the limit σ CSE (k → 0) considerably larger. To avoid this uncertainty of extrapolation to k = 0, in the present work we have compared our data with the free quark result for σ CSE (k) at the same nonzero momentum k min , but comparison at smaller values of k would be more reliable. The simulations were also performed on the GPU cluster at the Institute for Theoretical Physics at Giessen University. Many thanks to Dominik Schweitzer for keeping this cluster alive during #JLUoffline.