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.


Introduction
For the vacuum state of QCD as well as for the properties of finite temperature QCD the existence of non-trivial 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) lowdimensional (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 still 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 JHEP06(2015)094 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. But 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 2 we introduce the model and its lattice implementation and define the quantities we measure. In section 3 we present our results, and section 4 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.

Details of the simulations
We have performed simulations with the SU(2) gauge group. We employ the standard Wilson plaquette action

JHEP06(2015)094
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) x 2 , δ = (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.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.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 divergences. 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.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 divergences in observables. In the appendix we present analytical and numerical investigations of additional divergences in the Polyakov loop and the chiral condensate. Our study shows that there is no additional divergence in the Polyakov loop, whereas there is an additional logarithmic divergence in the chiral condensate. The latter is numerically small and does not effect the results of this paper.

JHEP06(2015)094
We have performed simulations with two lattice sizes N τ × N 3 σ = 6 × 20 3 , 10 × 28 3 . The measured observables are Tr Nτ n 4 =1 U 4 (n 1 , n 2 , n 3 , n 4 ) , (2.5) • the chiral condensate • the Polyakov loop susceptibility • the disconnected part of the chiral susceptibility (2.8) 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].   It is seen that the temperature of both phase transitions increases with the chiral chemical potential. One also sees that the phase transition 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  We have fitted the data for both the susceptibilities near the peak with a Gaussian func- and extracted the critical temperatures T χ c (β χ c ) using the chiral susceptibility and T L c (β L c ) using the Polyakov loop susceptibility. The resulting dependences of both critical temperatures on the value of the chiral chemical potential µ 5 are shown on the figure 4 and table 1. The values of β χ c , β L c and T χ c , T L 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. One sees that for all points except for µ 5 = 950 MeV the critical temperatures T χ c and T L c coincide within errors. There is a slight discrepancy for µ 5 = 950 MeV, but because of possible systematic uncertainties in the fit as well as finite size effects, we cannot claim that the transition temperatures are different.
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 µ 5 . The results for the Polyakov loop and the chiral condensate are presented in figure 5. They underpin the fact 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.
Polyakov loop Chiral condensate  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 figure 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].  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  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

JHEP06(2015)094
The fits performed with the eq. (3.1) are shown as lines in figure 7. The corresponding parameters of both are given in table 2. As can be seen from this table 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 2, showing that for these values of µ 5 the data are consistent with 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 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 (figure 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 behavior 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 described above.

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. with four flavor degrees of freedom in the continuum limit.
We have calculated the chiral condensate and the Polyakov loop for different values of the temperature T and chiral chemical potential µ 5 on lattices of size 6 × 20 3 and 10 × 28 3 as well as their respective susceptibilities (for the smaller lattice size). Our main result is the observation that the finite temperature phase transition becomes clearly shifted to larger critical temperature with rising chiral chemical potential. It was seen that this conclusion remains true when one is extrapolating to vanishing quark mass. Additionally by comparing the peak positions of the Polyakov loop susceptibility with those of the chiral susceptibility we saw the respective critical temperatures T L c and T χ c to agree within errors at least up to µ 5 ≈ 0.5 GeV. For larger µ 5 values we were not able to draw a final conclusion about the (dis)agreement.
Our result is in contradiction to results obtained with various effective models of QCD [21-23, 26, 27], where the critical temperature was said to decrease as µ 5 increases.
We concede that we study two-color QCD with N f = 4 quarks which is different from what is mostly considered in the effective models. For a closer comparison, it would be very interesting to examine the corresponding model case or to carry out the corresponding lattice simulation in the SU(3) case.
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 are in conflict with each other. For instance, the dependence of the chiral condensate on µ 5 obtained in papers [21,22] is opposite to the result obtained JHEP06(2015)094 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 ≈ 3.3 GeV, which would manifest itself as discontinuity in the 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 in the chiral limit there is an equivalence between the QCD phase diagram at finite µ 5 and the 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 the latter theory the pion condensate and the critical temperature T c of pion condensation increase with µ I . Despite the fact that we considered N c = 2 one can expect that the chiral condensate and T c in our theory also increase with µ 5 . We believe this is one more fact in favor of our results. to universality in large-N c QCD-like theories. 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 and 15-32-21117, by a grant of the president of the RF, MD-3215.2014.2, and by a grant of the FAIR-Russia Research Center. The work of A. Yu. Kotov was also partially supported by the Dynasty foundation.

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 |p| 2 = sin 2 (p 1 ) + sin 2 (p 2 ) + sin 2 (p 3 ), 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 flavor can be written as To calculate the integral in formula (A.2) 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]  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 (A.6). To see this, consider some Feynman graph of radiative corrections to formula (A.6). 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 nonzero µ 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 figure 11. In addition to the lattice measurement, in figure 11 we plot the contribution of logarithmic divergence (A.6) to the difference ψ ψ − ψ ψ µ 5 =0 .
From figure 11 one sees that the uncertainty of the calculation does not 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 (A.6) is not very large. So one can use it to estimate the contribution of the 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 (A.6). For instance, from figure 1 one can see that at temperature T = 200 MeV (close to the phase transition at µ 5 = 0) the difference ( ψ ψ µ 5 =300 MeV − ψ ψ µ 5 =0 )/T 3 ≈ 2. At the same time the additional ultraviolet divergence according to formula (A.6) 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.

B Ultraviolet divergences in the Polyakov loop
At the leading order approximation in the strong coupling constant g the Polyakov loop can be written in the following form 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. (B.1) it is clear that additional divergences connected to non-zero µ 5 at one-loop level can appear from the fermion selfenergy 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.
One can easily check that this expression is free from ultraviolet divergence. Moreover, a calculation shows that at large q 2 expression (B.3) behaves as Π 00 ( q) − Π 00 ( q)| µ 5 =0 ∼ µ 2 5 (c 1 log q 2 +c 2 ), where c 1 , c 2 are some constants. So, at one-loop level there is no ultraviolet divergence in the Polyakov loop which results from non-zero chiral chemical potential. Similarly to the previous section, one can argue that there is no ultraviolet divergence in the µ 5 = 0 contribution to the Polyakov loop at higher loops.
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 figure 12

JHEP06(2015)094
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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.