Two-Color QCD with Non-zero Chiral Chemical Potential

The phase diagram of two-color QCD with non-zero chiral chemical potential is studied by means of lattice simulation. We focus on the influence of a chiral chemical potential on the confinement/deconfinement phase transition and the breaking/restoration of chiral symmetry. The simulation is carried out with dynamical staggered fermions without rooting. The dependences of the Polyakov loop, the chiral condensate and the corresponding susceptibilities on the chiral chemical potential and the temperature are presented. The critical temperature is observed to increase with increasing chiral chemical potential.


I. INTRODUCTION
For the vacuum state of QCD as well as for the properties of finite-temperature QCD the existence of nontrivial topological excitations is important. Well known are instantons [1] as classical solutions in Euclidean space as tunneling events between topologically different zero-temperature vacua. The role of topology for the solution of the famous U A (1) problem has been recognized very early [2,3]. Presently, the anomalous breaking of the U A (1) symmetry above the deconfinement transition/crossover is under intense investigation in the lattice community (see, e.g., Ref. [4]).
It is known by now from lattice QCD at zero temperature that a (fractal) low-dimensional (laminar) topological vacuum structure is discernible at very fine resolution scale [5], while localized instanton-like structures, actually prevailing at an infrared scale, are believed to explain chiral symmetry breaking [6,7].
The gluon fields contributing to the path integral at finite temperature correspondingly may contain calorons [8,9]. Adapted to the non-trivial holonomy they have a richer structure (in terms of "dyonic" constituents) than instantons. The changes of this structure at the QCD phase transition are presently under study [10,11].
Some time ago the gluonic topological structure and the famous axial anomaly have been proposed to be immediately observable (and controllable) through the generation of P and CP violating domains (violating also translational invariance) in heavy ion collisions [12,13]. It has been demonstrated by detailed numerical calculations [12,14] that macroscopic domains of (anti)parallel color-electric and color-magnetic field can emerge in a heavy ion collision creating an increasing chiral imbalance among the quarks which are deconfined due to the high energy density. In this situation, the magnetic field created by the spectator nucleons may initiate a charge separation relative to the reaction plane (parallel to the electro-magnetic field) [15]. The resulting charge asymmetry of quarks would become observable in terms of recombined hadrons (chiral-magnetic effect) [16,17]. The strength (and particularly the dependence on the collision energy) of this effect has been theoretically studied and proposed to be a signal of the transient existence of liberated quarks [12,13,18].
In recent years the dependence of the chiral and deconfinement transitions on the magnetic field has been investigated both in models and ab-initio lattice simulations, see e.g. [19,20]. It remains an open question whether the phase transition from quarks to hadrons, i.e. the onset of confinement and chiral symmetry breaking (and vice versa), depends on the chiral imbalance.
In this paper we study the change of the phase structure induced by a chiral chemical potential in an equilibrium lattice simulation. We mimic the topological content (of a CP violating, topologically nontrivial gluonic background in a heavy ion collision event) by inducing a chiral imbalance, which is provided (i.e. frozen) by a non-zero chiral chemical potential. In this setting, the modification of the phase diagram by the chiral chemical potential µ 5 has been studied mainly by means of effective models [21][22][23][24][25][26][27], the predictions of which will later be compared with our results.
On the lattice, contrary to the case of non-zero baryon chemical potential, simulations with non-vanishing µ 5 are not hampered by a sign problem 1 . Thus, one can employ the standard hybrid Monte Carlo algorithms. Such lattice simulations with µ 5 = 0 were already performed in Ref. [28,29]. The main goal of these papers, however, was the study of the chiral magnetic effect. Therefore, the phase diagram was not systematically studied.
In Refs. [30,31] we have carried out the first lattice study of the phase diagram with non-zero chiral chemical potential. It was performed in SU (2) QCD with four flavors, which we have considered as a simplified model of QCD. In this article we extend this investigation generating considerably more data and performing a more detailed analysis. In this paper we do not attempt to analyse the topological structure that would reflect the presence of the chiral chemical potential.
One reason for choosing the SU (2) gauge group is that less computational resources are required for this pilot study than for full QCD. The second reason is that we have already carried out two-colour QCD computations with an external magnetic field [32,33]. Furthermore, the four flavor case results from our choice of staggered fermions as lattice regularization while we avoid to take the root of the fermion determinant, which would allow to reduce the number of flavors. The "rooting procedure" might introduce further systematic errors at finite lattice spacing.
In section II we introduce the model and its lattice implementation and define the quantities we measure. In section III we present our results, and section IV is devoted to their discussion and to the formulation of conclusions. In the Appendix we discuss the question of renormalizations refering to explicit analytical calculations in perturbation theory.

II. DETAILS OF THE SIMULATIONS
We have performed simulations with the SU (2) gauge group. We employ the standard Wilson plaquette action For the fermionic part of the action we use staggered fermions where the η µ (x) are the standard staggered phase factors: The lattice spacing is denoted by a, the bare fermion mass by m, and µ 5 is the value of the chiral chemical potential. In the chirality breaking term s(x) = (−1) x2 , δ = (1, 1, 1, 0) represents a shift to the diagonally opposite site in a spatial 2 3 elementary cube. The combination of three links connecting sites x and x + δ, is symmetrized over the 6 shortest paths between these sites. In the partition function, after formally integrating over the fermions, one obtains the corresponding determinant. As mentioned above, we do not take the fourth root of this determinant in order to represent each flavor ("taste") independently of the others. Thus, the continuum limit of our model corresponds to a theory of four (degenerate) flavors.
In the continuum limit Eq. (2) can be rewritten in the Dirac spinor-flavor basis [34,35] as follows We would like to emphasize that the chiral chemical potential, introduced in Eq. (2), corresponds to the taste-singlet operator γ 5 γ 4 ⊗ 1 in the continuum limit.
It should be also noted here that the usual baryonic chemical potential after the discussion in Ref. [36], and also the chiral chemical potential as used in Ref. [29], are introduced to the action as a modification of the temporal links by corresponding exponential factors in order to eliminate chemical-potential dependent quadratic divergencies. For staggered fermions this modification can be performed as well for the baryonic chemical potential. However, for the chiral chemical potential such a modification would lead to a highly nonlocal action [29]. Therefore, we decided to introduce µ 5 in Eq. (2) in an additive way similar to the mass term.
It is known that the additive way of introducing the chemical potential leads to additional divergencies in observables. In the Appendix we present analytical and numerical investigations of additional divergencies in the Polyakov loop and the chiral condensate. Our study shows that there is no additional divergency in the Polyakov loop, whereas there is an additional logarithmic divergency in the chiral condensate. The latter is numerically small and does not effect the results of this paper.
We have performed simulations with two lattice sizes N τ × N 3 σ = 6 × 20 3 , 10 × 28 3 . The measured observables are Tr Nτ n4=1 U 4 (n 1 , n 2 , n 3 , n 4 ) , • the chiral condensate • the Polyakov loop susceptibility • the disconnected part of the chiral susceptibility The Polyakov loop and the corresponding susceptibility are sensitive to the confinement/deconfinement phase transition, whereas the chiral condensate in principle responds to chiral symmetry breaking/restoration. The simulations have been carried out with a CUDA code in order to perform the simulations using the Hybrid Monte Carlo algorithm on GPU's.
The parameters of our lattice calculation are the inverse of the bare coupling constant β, the bare mass ma and the bare chiral chemical potential aµ 5 (both in dimensionless units) as well as the lattice size. The physical temperature and the volume are given by To perform the scale setting (calibration) of the lattice, we use the results of Ref. [32]. There, the static potential in the same theory has been measured at zero temperature for several values of β. The Sommer parameter r 0 was calculated in lattice units and compared to its physical value, which was supposed to be the same as in QCD, i.e. r 0 = 0.468(4) fm. It was found that the scaling function a(β) in the region were we perform our measurements is well described by the two loop β-function. Thus, for given β we can obtain e.g. the temperature T in units of MeV. For more details, see Ref. [32].  Fig. 1. The sharp change of the observables as functions of T indicates the onset of the deconfinement and the chiral restoration phase transition. It is seen that the temperature of both phase transitions increases with the chiral chemical potential. One also sees that the phase trasition becomes sharper for increasing chiral chemical potential. To study the change of the critical temperature more quantitatively, we also calculated the chiral and the Polyakov loop susceptibilities. In order to make the figures readable we plot separately the susceptibilities for values µ 5 = 0, 475, 950 MeV in Fig. 2 and for values µ 5 = 0, 150, 300 MeV in Fig. 3. We see that increasing the value of the chiral chemical potential moves the position of the peaks of the chiral and Polyakov loop susceptibilities to larger values of β. This means that the transition temperature increases. Our results do not show any splitting between the chiral and the deconfinement transition. We have fitted the data for the chiral susceptibility near the peak with a gaussian function χ = a 0 + a 1 exp − (β−βc) 2 2σ 2 and extracted the critical temperature T c (β c ). The resulting dependence of the critical temperature on the value of the chiral chemical potential µ 5 is shown on the Fig. 4 and Table I. The values of β c and T c and their uncertainties are calculated from a fit of 5-6 points in the vicinity of the peak by the Gaussian function given above. It should be noted that for small values of µ 5 the dependence of T c on µ 5 is well described by the function T c = a + bµ 2 5 , but the larger µ 5 is the larger is the deviation from this simple formula.   In addition to the calculations on the 6 × 20 3 lattice we carried out simulations on a larger lattice of size 10 × 28 3 . Note that this allows us to investigate larger values of µ 5 . The susceptibilities require large statistics, and therefore our current computational resources do not allow us to measure them on the larger lattice. In the simulations we kept the physical fermion mass fixed at m = 18.5 MeV (m π ≈ 550 MeV).
On this lattice we also calculated the observables L and ψψ as functions of T for different values of that the critical temperature grows when one increases µ 5 . The value of the chiral chemical potential used in our simulations was rather large (up to µ 5 = 3345 MeV), but we have not seen signals of a first-order phase transition, although the transition becomes sharper when the value of the chiral chemical potential grows.
It is interesting to study how the observables L and ψψ depend on µ 5 for fixed temperature. To do this we have measured the Polyakov loop and the chiral condensate as functions of µ 5 a for three different values of β, β = 1.87 (158 MeV), 1.91 (186 MeV) and 1.95 (219 MeV), which for a vanishing chiral chemical potential corresponds to a temperature below the transition, in the transition region, and in the high temperature phase, respectively. The results of these measurements are shown in Fig. 6. As can be seen from the figure, in the confinement phase the Polyakov loop remains almost constant with increasing chiral chemical potential. It means that if the system was in the confinement phase at µ 5 = 0, it remains confined at µ 5 > 0. Moreover, we observe the Polyakov loop to drop down with increasing µ 5 both in the deconfinement phase and in the transition region. Thus, the system goes into the confinement phase for sufficiently large µ 5 . With other words, we conclude that the critical temperature increases with an increasing chiral chemical potential in agreement with our results obtained on the smaller lattice. It is worth mentioning that the behavior described above looks quite similar to the behavior obtained for two-color QCD in an external magnetic field [10,32]. In Fig. 7 we show the results in the confinement region at β = 1.80 (T = 119 MeV). The Polyakov loop is small and does not show any nontrivial behaviour. The chiral condensate remains almost constant when the value of the fermion mass changes.
To extrapolate to the chiral limit (ma = 0) we performed chiral extrapolations in two ways. Since the value of the temperature is not very far from the transition temperature, we suppose that we can use the formula for the reduced three dimensional model [37]. See also [38][39][40]. In that case we may use the ansatz a 3 <ψψ >= f 1 (ma),   where with the non-analytic term coming from the Goldstone bosons. Such a parametrization has been used in [41] in the context of the finite temperature transition in QCD. As an alternative we also use the chiral extrapolation relevant for zero temperature, namely a 3 <ψψ >= f 2 (ma), where The fits performed with the Eq. (10) are shown as lines in Fig. 7. The corresponding parameters of both are given in Table II. As can be seen from table II and the corresponding figure, the chiral condensate extrapolates to a non-zero value in the chiral limit.
At larger values β = 1.90 (T = 178 MeV), β = 2.00 (T = 268 MeV), and β = 2.10 (T = 404 MeV) and at zero chiral chemical potential the system is in the chirally restored phase, when m goes to zero. However, at β = 1.90 (T = 178 MeV) and µ 5 = 1115, 2230 MeV, the chiral condensate has non-zero expectation values. We have performed the same extrapolation technique as for β = 1.80 and present also these results in Table II, showing that for these values of µ 5 the data are consistent with a non-zero values in the chiral limit. It implies that also in the chiral limit a non-zero chiral chemical potential shifts the position of the phase transition to larger temperatures. The plot for the Polyakov loop (right panel of Fig. 8) confirms this observation: for the two smaller values of the fermion mass ma the values of Polyakov loop for µ 5 a = 0 is several times larger than for non-zero values of µ 5 , so that the system for non-zero µ 5 is deeper in the confinement region than for zero chiral chemical potential.
At β = 2.00 (T = 268 MeV) and β = 2.10 (T = 404 MeV) the chiral condensate goes to zero in the chiral limit (Fig. 9, 10), thus the system is in the chirally restored phase. The dependence of the chiral condensate on the value of the mass is almost a linear function, apart from the behavior for β = 2.00 (T = 268 MeV) at the largest chiral chemical potential under our investigation, µ 5 = 2230 MeV. This behaviour can be explained if one assumes that this point is near the phase transition. It is known that increasing the mass shifts the position of the transition to larger temperatures [32]. Thus, at smaller masses the system is in the chirally restored phase, while increasing the mass lets the system undergo the chiral phase transition. Again, the plot of the Polyakov loop confirms this suggestion. The Polyakov loop is small at larger masses implying that the system is in the confinement phase and increases when the mass is decreased, i.e. the system becomes deconfining. This behavior suggests that at the largest value of µ 5 = 2230 MeV the transition happens at even larger temperature than we have investigated. These results are in total agreement with the results above.

IV. DISCUSSION AND CONCLUSION
In this paper we have presented an investigation of the phase diagram of two-color QCD with a chiral chemical potential using lattice simulations with dynamical staggered fermions without rooting, i.e. four flavors in the continuum limit. We have calculated the chiral condensate, the Polyakov loop, the Polyakov loop susceptibility and the chiral susceptibility for different values of the temperature T and the chiral chemical potential µ 5 on lattices of size 6 × 20 3 and 10 × 28 3 . The main result is that at non-zero values of the chiral chemical potential the critical temperatures of the confinement/deconfinement phase transition and of the chiralsymmetry breaking/restoration phase transition still coincide and that the common transition shifts to larger temperatures as µ 5 increases. It is shown that this conclusion remains true in the limit of zero quark mass. This result is in contradiction to the results within different effective models of QCD [21-23, 26, 27] where the critical temperature is said to decrease as µ 5 increases.
We concede that we study SU (2) QCD with N f = 4 quarks which is different from what is considered in the effective models, where SU(3) QCD with N f = 2 quarks is considered. It would therefore be very interesting for a closer comparison to examine the results of these models by a lattice simulation.
It is likely, however, that the contradictions are rather a consequence of the fact that the critical temperature is not a universal parameter but crucially depends on the structure and parameters of the effective models. It is unclear to what extent they describe the actual behavior of finite temperature QCD. It should be also noted that some predictions obtained in different effective models contradict among each others. For instance, the dependence of the chiral condensate on µ 5 obtained in papers [21,22] is opposite to the result obtained in paper [25]. The latter paper claims that the chiral condensate grows with the chiral chemical potential, which is in agreement with our results.
In addition to the dependence of the critical temperature on the chiral chemical potential, some effective models predict that -beginning from some critical value of µ c 5 (µ c 5 ≈ 400 MeV in Ref. [21], µ c 5 = 50 MeV in Ref. [22]) -the transition from the hadronic to the plasma phase turns into a first order transition. In our simulations we see that the transition becomes sharper as we increase µ 5 , but we don't see a first order phase transition up to µ 5 = 3345 MeV, which would manifest itself as discontinuity in Polyakov loop and chiral condensate at a sufficiently large spatial extension of the lattice. In this connection a finite size scaling analysis would be valuable but is outside the scope of this paper.
Besides within effective models, the phase diagram of QCD in the (µ 5 , T )-plane was studied in papers [42,43] in a framework of Dyson-Schwinger equations. The authors of these papers found that the critical temperature rises with µ 5 and the "phase transition" actually is always a crossover. These results are corroborating the results of our paper.
We would like also to mention the paper [44]. In this paper the authors address the question of universality of phase diagrams in QCD and QCD-like theories through the large-N c equivalence. Using the results of this paper one can show that at large-N c and the chiral limit there is equivalence between QCD phase diagram at finite µ 5 and QCD phase diagram at finite isospin chemical potential µ I . The chiral condensate of QCD with µ 5 = 0 can be mapped to the pion condensate of QCD with µ I = 0. In latter theory the pion condensate and critical temperature T c of pion condensation increase with µ I . Despite the fact that we considered N c = 2 one can expect that chiral condensate and T c in our theory also increase with with µ 5 . We believe this is one more fact in favour of our results.
The simulations were performed at GPUs of the K100 supercomputer of the Institute of Applied Mathematics of the Russian Academy of Sciences in Moscow and at GPUs of the Particle Phenomenology group at the Institute of Physics of the Humboldt University Berlin. The work was supported by Far Eastern Federal University, by RFBR grants 13-02-01387-a, 14-02-01185-a, 15-02-07596-a, by a grant of the president of the RF, MD-3215.2014.2, and by a grant of the FAIR-Russia Research Center.
Appendix A: Ultraviolet divergences in the chiral condensate The fermion propagator including the chiral chemical potential for "naive" lattice fermions can be written in the following form Here m and µ 5 are mass and chiral chemical potential in lattice units, α, β are color indices, the sum is taken over all possible values of (n 1 , n 2 , n 3 , n 4 ) and s = ±1.
In the limit L s , L t → ∞ the condensate per one fermion flavour can be written as To calculate the integral in formula (A2) we use the standard approach, which separates the main divergence from the rest of the integral The first term in this expression is just the loop integral for ψ ψ without chiral chemical potential. The expression for this integral in the continuum limit a → 0 can be found in [45] π −π The second term in (A3) is proportional to µ 2 5 and contains a logarithmic divergence in the continuum limit. This divergence can be calculated using saddle-point method. To calculate the second term in (A3) we subtract the leading divergence and calculate the rest of the integral in the limit a → 0. The result of the calculation can be written as follows We have checked the last formula numerically. Introducing a physical mass m p a = m and a physical chiral chemical potential µ p 5 a = µ 5 one can write the expression for the condensate in physical units as From the last formula one sees that at one-loop level the inclusion of a non-zero chiral chemical potential leads to an additional logarithmic divergence. We believe that this conclusion persists if one takes into account radiative corrections to formula (A6). To see this, consider some Feynman graph of radiative corrections to formula (A6). The expansion of this graph in powers of µ 5 can be reduced to the inclusion of a dimension-3 operator to the graph that diminishes the power of divergence by one unit per one power of µ 5 . Expansion of the chiral condensate in powers of µ 5 contains only even powers. The main divergence of the chiral condensate is ∼ 1/a 2 . So, the next to leading term in µ 5 expansion is two powers smaller than ∼ 1/a 2 , i.e. it is at most logarithmically divergent. Further it is important to understand how the additional divergence connected to non-zero µ 5 effects the results of this paper. To estimate this we fixed the physical values of temperature T = 270 MeV, quark mass m p = 33 MeV and measured the condensate ψ ψ in the deconfinement phase for µ 5 = 0 and µ 5 = 330 MeV for different values of lattice spacing a. In particular, we took the following lattice parameters: 16 3 × 8 β = 1.9500, 20 3 × 10 β = 2.0047, 24 3 × 12 β = 2.0493 and 28 3 × 14 β = 2.0870. The difference between the measured chiral condensates at zero and non-zero µ 5 for different lattice spacings is shown in Fig.11. In addition to the lattice measurement, in Fig.11 we plot the contribution of logarithmic divergence (A6) to the difference ψ ψ − ψ ψ µ5=0 . From Fig.11 one sees that uncertainty of the calculation doesn't allow to confirm the logarithmic type of the divergence. However, the variation of the lattice results with variation of the lattice spacing is very slow what allows us to assume that there are no other divergences different from the logarithmic one. Note also that the value of the µ 5 = 0 contribution to the chiral condensate obtained from the one-loop expression for the logarithmic divergence is by a factor 2-3 smaller than the lattice results. The difference can be attributed to a non-zero temperature effect and radiative corrections. This fact allows one to state that the renormalization effect of the one-loop expression for the logarithmic divergence (A6) is not very large. So one can use it to estimate the contribution of ultraviolet logarithmic divergence in the confinement phase.
The calculation shows that the characteristic values of the difference ψ ψ − ψ ψ µ5=0 in the confinement phase for different values of lattice parameters used in this paper are approximately by one order of magnitude larger than the additional ultraviolet logarithmic divergence due to µ 5 = 0 estimated according to formula (A6). For instance, from Fig. 1 one can see that at temperature T = 200 MeV (close to the phase transition at µ 5 = 0) the difference ( ψ ψ µ5=300MeV − ψ ψ µ5=0 )/T 3 ≈ 2. At the same time the additional ultraviolet divergence according to formula (A6) gives ( ψ ψ add )/T 3 ≈ 0.1. Note also that the position of the phase transition manifests itself as a peak of the susceptibility. Evidently, there is no peak due to the additional ultraviolet divergence. All this allows us to state that the conclusions obtained in this paper are not affected by the additional ultraviolet µ 5 = 0 divergence in the chiral condensate.
where D ab 00 (t, x) andD aa 00 (q 0 , q) represent the propagator of the temporal components of the gluon field in coordinate and momentum space, correspondingly, β is inverse temperature, and the summation over color index a is assumed. From eq. (B1) it is clear that additional divergences connected to non-zero µ 5 at one-loop level can appear from the fermion self-energy part of the gluon propagator. The study of the divergences with "naive" fermions in the self-energy part of the propagator is cumbersome. For this reason we use momentum cut of regularization procedure to study the divergences in the Polyakov loop. The fermion contribution to the one-loop self-energy of the gluon propagatorD ab 00 can be written in the following form Π 00 ( q) = 1 4 s1,s2 Λ d 4 k (2π) 4 −k 2 0 + (µ 5 − s 1 | k|)(µ 5 − s 2 | k + q|) + m 2 (k 2 0 + (| k| − s 1 µ 5 ) 2 + m 2 )(k 2 0 + (| k + q| − s 2 µ 5 ) 2 + m 2 ) 1 + s 1 s 2 k · ( k + q) where q is the external momentum, s 1 , s 2 = ±1.
In order to numerically study the role of divergences in the Polyakov loop due to non-zero chiral chemical potential, we used the following definition of renormalized Polyakov loop where L is the Polyakov loop and L µ5=0 is the Polyakov loop at the same lattice but with µ 5 = 0. Evidently L ren doesn't contain divergences that are usual for simulations with zero chiral chemical potential. Similarly to the previous section we fixed the physical values of temperature T = 270 MeV, quark mass m p = 33 MeV and measured the L ren in the deconfinement phase for µ 5 = 330 MeV and different values of lattice spacing a.
The result of this measurement is shown in Fig. 12 From this figure one sees that the renormalized Polyakov loop L ren is consistent with a constant in the region of lattice spacings investigated. In the same region the unrenormalized Polyakov loop L changes by a factor of two. From these facts we conclude that there is no additional ultraviolet divergence in Polyakov loop due to non-zero chiral chemical potential.