TwoColor QCD with Nonzero Chiral Chemical Potential
Abstract
The phase diagram of twocolor QCD with nonzero 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.
pacs:
12.38.Aw, 12.38.Gc, 11.15.HaITEPLAT201501, HUEP15/12
I Introduction
For the vacuum state of QCD as well as for the properties of finitetemperature QCD the existence of nontrivial topological excitations is important. Well known are instantons Belavin et al. (1975) as classical solutions in Euclidean space as tunneling events between topologically different zerotemperature vacua. The role of topology for the solution of the famous problem has been recognized very early Witten (1979); Veneziano (1979). Presently, the anomalous breaking of the symmetry above the deconfinement transition/crossover is under intense investigation in the lattice community (see, e.g., Ref. Dick et al. (2015)).
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 Ilgenfritz et al. (2008), while localized instantonlike structures, actually prevailing at an infrared scale, are believed to explain chiral symmetry breaking Schäfer and Shuryak (1998); Shuryak and Schäfer (1997).
The gluon fields contributing to the path integral at finite temperature correspondingly may contain calorons Kraan and van Baal (1998); Lee and Lu (1998). Adapted to the nontrivial 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 Ilgenfritz et al. (2014a); Bornyakov et al. (2014).
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 and violating domains (violating also translational invariance) in heavy ion collisions Kharzeev (2006); Fukushima et al. (2008). It has been demonstrated by detailed numerical calculations Kharzeev (2006); Kharzeev et al. (2002) that macroscopic domains of (anti)parallel colorelectric and colormagnetic 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 electromagnetic field) Kharzeev et al. (2008). The resulting charge asymmetry of quarks would become observable in terms of recombined hadrons (chiralmagnetic effect) Abelev et al. (2009, 2010). 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 Kharzeev (2006); Fukushima et al. (2008); Adamczyk et al. (2014).
In recent years the dependence of the chiral and deconfinement transitions on the magnetic field has been investigated both in models and abinitio lattice simulations, see e.g. Shovkovy (2013); D’Elia (2013). 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 violating, topologically nontrivial gluonic background in a heavy ion collision event) by inducing a chiral imbalance, which is provided (i.e. frozen) by a nonzero chiral chemical potential. In this setting, the modification of the phase diagram by the chiral chemical potential has been studied mainly by means of effective models Fukushima et al. (2010); Chernodub and Nedelin (2011); Gatto and Ruggieri (2012); Andrianov et al. (2014); Planells et al. (2013); Chao et al. (2013); Yu et al. (2014), the predictions of which will later be compared with our results.
On the lattice, contrary to the case of nonzero baryon chemical potential,
simulations with nonvanishing are not hampered by a sign
problem
In Refs. Braguta et al. (2015, 2014) we have carried out the first lattice study of the phase diagram with nonzero chiral chemical potential. It was performed in 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 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 twocolour QCD computations with an external magnetic field Ilgenfritz et al. (2012); Ilgenfritz et al. (2014b). 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 gauge group. We employ the standard Wilson plaquette action
(1) 
For the fermionic part of the action we use staggered fermions
(2) 
where the are the standard staggered phase factors: for . The lattice spacing is denoted by , the bare fermion mass by , and is the value of the chiral chemical potential. In the chirality breaking term , represents a shift to the diagonally opposite site in a spatial elementary cube. The combination of three links connecting sites and ,
(3) 
is symmetrized over the 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 spinorflavor basis KlubergStern et al. (1983); Montvay and Münster (1994) as follows
(4) 
We would like to emphasize that the chiral chemical potential, introduced in Eq. (2), corresponds to the tastesinglet operator in the continuum limit.
It should be also noted here that the usual baryonic chemical potential after the discussion in Ref. Hasenfratz and Karsch (1983), and also the chiral chemical potential as used in Ref. Yamamoto (2011b), are introduced to the action as a modification of the temporal links by corresponding exponential factors in order to eliminate chemicalpotential 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 Yamamoto (2011b). Therefore, we decided to introduce 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 . The measured observables are

the Polyakov loop
(5) 
the chiral condensate
(6) 
the Polyakov loop susceptibility
(7) 
the disconnected part of the chiral susceptibility
(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 and the bare chiral chemical potential (both in dimensionless units) as well as the lattice size. The physical temperature and the volume are given by
(9) 
Polyakov loop  chiral condensate 
To perform the scale setting (calibration) of the lattice, we use the results of Ref. Ilgenfritz et al. (2012). There, the static potential in the same theory has been measured at zero temperature for several values of . The Sommer parameter was calculated in lattice units and compared to its physical value, which was supposed to be the same as in QCD, i.e. fm. It was found that the scaling function 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 in units of MeV. For more details, see Ref. Ilgenfritz et al. (2012).
Iii Results of the calculation
We first performed simulations on a lattice of size for five fixed values of MeV and different values of . The fermion mass was kept fixed in physical units MeV( MeV). The expectation values of the Polyakov loop and the chiral condensate are shown in Fig. 1. The sharp change of the observables as functions of 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.
Polyakov loop susceptibility  chiral susceptibility 
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 MeV in Fig. 2 and for values 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 and extracted the critical temperature . The resulting dependence of the critical temperature on the value of the chiral chemical potential is shown on the Fig. 4 and Table 1. The values of and and their uncertainties are calculated from a fit of 56 points in the vicinity of the peak by the Gaussian function given above. It should be noted that for small values of the dependence of on is well described by the function , but the larger is the larger is the deviation from this simple formula.
Polyakov loop susceptibility  chiral susceptibility 
In addition to the calculations on the lattice we carried out simulations on a larger lattice of size . Note that this allows us to investigate larger values of . 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 MeV ( MeV).
Polyakov loop  chiral condensate 
On this lattice we also calculated the observables and as functions of for different values of . The results for the Polyakov loop and the chiral condensate are presented in Fig. 5. They underpin the fact that the critical temperature grows when one increases . The value of the chiral chemical potential used in our simulations was rather large (up to MeV), but we have not seen signals of a firstorder phase transition, although the transition becomes sharper when the value of the chiral chemical potential grows.
It is interesting to study how the observables and depend on for fixed temperature. To do this we have measured the Polyakov loop and the chiral condensate as functions of for three different values of , (158 MeV), (186 MeV) and (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 , it remains confined at . Moreover, we observe the Polyakov loop to drop down with increasing both in the deconfinement phase and in the transition region. Thus, the system goes into the confinement phase for sufficiently large . 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 twocolor QCD in an external magnetic fieldIlgenfritz et al. (2012); Ilgenfritz et al. (2014a).
Polyakov loop  chiral condensate 
In order to study the dependence of the observables on the value of the quark mass for different values of chiral chemical potential, we chose four values of : one in the confinement phase, (corresponding to a temperature MeV), two in the transition region, ( MeV) and ( MeV), and one in the deconfinement phase ( MeV). The lattice size used here is . We chose four values of the fermion mass and three values of the chemical potential MeV, MeV and MeV. The expectation values of the Polyakov loop and the chiral condensate are presented as functions of for different in Fig. 7, Fig. 8, Fig. 9, and Fig. 10.
In Fig. 7 we show the results in the confinement region at ( 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 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 Pisarski and Wilczek (1984). See also Wallace and Zia (1975); Hasenfratz and Leutwyler (1990); Smilga and Stern (1993). In that case we may use the ansatz , where
(10) 
with the nonanalytic term coming from the Goldstone bosons. Such a parametrization has been used in Bazavov et al. (2012) in the context of the finite temperature transition in QCD. As an alternative we also use the chiral extrapolation relevant for zero temperature, namely , where
(11) 
The fits performed with the Eq. (10) are shown as lines in Fig. 7. The corresponding parameters of both are given in Table 2. As can be seen from table II and the corresponding figure, the chiral condensate extrapolates to a nonzero value in the chiral limit.
At larger values ( MeV), ( MeV), and ( MeV) and at zero chiral chemical potential the system is in the chirally restored phase, when goes to zero. However, at ( MeV) and MeV, the chiral condensate has nonzero expectation values. We have performed the same extrapolation technique as for and present also these results in Table 2, showing that for these values of the data are consistent with a nonzero values in the chiral limit. It implies that also in the chiral limit a nonzero 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 the values of Polyakov loop for is several times larger than for nonzero values of , so that the system for nonzero is deeper in the confinement region than for zero chiral chemical potential.
At ( MeV) and ( 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 ( MeV) at the largest chiral chemical potential under our investigation, 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 Ilgenfritz et al. (2012). 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 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 twocolor 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 and the chiral chemical potential on lattices of size and . The main result is that at nonzero 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 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 Fukushima et al. (2010); Chernodub and Nedelin (2011); Gatto and Ruggieri (2012); Chao et al. (2013); Yu et al. (2014) where the critical temperature is said to decrease as increases.
We concede that we study SU(2) QCD with quarks which is different from what is considered in the effective models, where SU(3) QCD with 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 obtained in papers Fukushima et al. (2010); Chernodub and Nedelin (2011) is opposite to the result obtained in paper Planells et al. (2013). 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 ( MeV in Ref. Fukushima et al. (2010), MeV in Ref. Chernodub and Nedelin (2011)) – 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 , but we don’t see a first order phase transition up to 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 –plane was studied in papers Wang et al. (2015); Xu et al. (2015) in a framework of DysonSchwinger equations. The authors of these papers found that the critical temperature rises with 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 Hanada and Yamamoto (2011). In this paper the authors address the question of universality of phase diagrams in QCD and QCDlike theories through the large equivalence. Using the results of this paper one can show that at large and the chiral limit there is equivalence between QCD phase diagram at finite and QCD phase diagram at finite isospin chemical potential . The chiral condensate of QCD with can be mapped to the pion condensate of QCD with . In latter theory the pion condensate and critical temperature of pion condensation increase with . Despite the fact that we considered one can expect that chiral condensate and in our theory also increase with with . We believe this is one more fact in favour of our results.
Acknowledgements
The authors are grateful to V. I. Zakharov and V. G. Bornyakov for interesting and stimulating discussions. The authors are grateful to N. Yamamoto who paid our attention to universality in large 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 130201387a, 140201185a, 150207596a, by a grant of the president of the RF, MD3215.2014.2, and by a grant of the FAIRRussia 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
(12)  
Here and are mass and chiral chemical potential in lattice units, are color indices, the sum is taken over all possible values of and .
In the limit the condensate per one fermion flavour can be written as
(13) 
To calculate the integral in formula (13) we use the standard approach, which separates the main divergence from the rest of the integral
(14)  
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 can be found in Capitani (2003)
(15) 
The second term in (14) is proportional to and contains a logarithmic divergence in the continuum limit. This divergence can be calculated using saddlepoint method. To calculate the second term in (14) we subtract the leading divergence and calculate the rest of the integral in the limit . The result of the calculation can be written as follows
(16) 
We have checked the last formula numerically. Introducing a physical mass and a physical chiral chemical potential one can write the expression for the condensate in physical units as
(17) 
From the last formula one sees that at oneloop level the inclusion of a nonzero chiral chemical potential leads to an additional logarithmic divergence. We believe that this conclusion persists if one takes into account radiative corrections to formula (17). To see this, consider some Feynman graph of radiative corrections to formula (17). The expansion of this graph in powers of can be reduced to the inclusion of a dimension3 operator to the graph that diminishes the power of divergence by one unit per one power of . Expansion of the chiral condensate in powers of contains only even powers. The main divergence of the chiral condensate is . So, the next to leading term in expansion is two powers smaller than , i.e. it is at most logarithmically divergent.
Further it is important to understand how the additional divergence connected to nonzero effects the results of this paper. To estimate this we fixed the physical values of temperature MeV, quark mass MeV and measured the condensate in the deconfinement phase for and MeV for different values of lattice spacing . In particular, we took the following lattice parameters: , , and . The difference between the measured chiral condensates at zero and nonzero 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 (17) to the difference .
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 contribution to the chiral condensate obtained from the oneloop expression for the logarithmic divergence is by a factor 23 smaller than the lattice results. The difference can be attributed to a nonzero temperature effect and radiative corrections. This fact allows one to state that the renormalization effect of the oneloop expression for the logarithmic divergence (17) 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 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 estimated according to formula (17). For instance, from Fig. 1 one can see that at temperature MeV (close to the phase transition at ) the difference . At the same time the additional ultraviolet divergence according to formula (17) gives . 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 divergence in the chiral condensate.
Appendix B Ultraviolet divergences in the Polyakov loop
At the leading order approximation in the strong coupling constant the Polyakov loop can be written in the following form
(18) 
where and 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 is assumed. From eq. (18) it is clear that additional divergences connected to nonzero at oneloop level can appear from the fermion selfenergy part of the gluon propagator. The study of the divergences with ”naive” fermions in the selfenergy 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 oneloop selfenergy of the gluon propagator can be written in the following form
(19) 
where is the external momentum, .
To proceed we expand equation (19) in . Evidently each new term in this expansion diminishes the power of divergence by two units. That means that beginning from term the series is ultraviolet convergent. The first term is just the expression for the selfenergy without chiral chemical potential. So, to calculate the divergence which results from nonzero one needs to study only the divergence of the second term. Taking the derivative of equation (19) with respect to we get
(20)  
One can easily check that this expression is free from ultraviolet divergence. Moreover, a calculation shows that at large expression (20) behaves as , where are some constants. So, at oneloop level there is no ultraviolet divergence in the Polyakov loop which results from nonzero chiral chemical potential. Similarly to the previous section, one can argue that there is no ultraviolet divergence in the contribution to the Polyakov loop at higher loops.
In order to numerically study the role of divergences in the Polyakov loop due to nonzero chiral chemical potential, we used the following definition of renormalized Polyakov loop
(21) 
where is the Polyakov loop and is the Polyakov loop at the same lattice but with . Evidently 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 MeV, quark mass MeV and measured the in the deconfinement phase for MeV and different values of lattice spacing . The result of this measurement is shown in Fig. 12
From this figure one sees that the renormalized Polyakov loop is consistent with a constant in the region of lattice spacings investigated. In the same region the unrenormalized Polyakov loop changes by a factor of two. From these facts we conclude that there is no additional ultraviolet divergence in Polyakov loop due to nonzero chiral chemical potential.
Footnotes
 In there is no sign problem even in presence of a baryon chemical potential.
References
 A. A. Belavin, A. M. Polyakov, A. S. Shvarts, and Y. S. Tyupkin, Phys. Lett. B59, 85 (1975).
 E. Witten, Nucl. Phys. B156, 269 (1979).
 G. Veneziano, Nucl. Phys. B159, 213 (1979).
 V. Dick, F. Karsch, E. Laermann, S. Mukherjee, and S. Sharma (2015), eprint 1502.06190.
 E.M. Ilgenfritz, D. Leinweber, P. Moran, K. Koller, G. Schierholz, and V. Weinberg, Phys.Rev. D77, 074502 (2008), eprint 0801.1725.
 T. Schäfer and E. V. Shuryak, Rev.Mod.Phys. 70, 323 (1998), eprint hepph/9610451.
 V. Shuryak, E and T. Schäfer, Nucl.Phys.Proc.Suppl. 53, 472 (1997).
 T. C. Kraan and P. van Baal, Nucl.Phys. B533, 627 (1998), eprint hepth/9805168.
 K.M. Lee and C.H. Lu, Phys.Rev. D58, 025011 (1998), eprint hepth/9802108.
 E.M. Ilgenfritz, B. Martemyanov, and M. MüllerPreussker, Phys.Rev. D89, 054503 (2014a), eprint 1309.7850.
 V. G. Bornyakov, E.M. Ilgenfritz, B. V. Martemyanov, and M. MüllerPreussker (2014), eprint 1410.4632.
 D. E. Kharzeev, Phys.Lett. B633, 260 (2006), eprint hepph/0406125.
 K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys.Rev. D78, 074033 (2008), eprint 0808.3382.
 D. E. Kharzeev, A. Krasnitz, and R. Venugopalan, Phys.Lett. B545, 298 (2002), eprint hepph/0109253.
 D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, Nucl.Phys. A803, 227 (2008), eprint 0711.0950.
 B. I. Abelev et al. (STAR Collaboration), Phys.Rev.Lett. 103, 251601 (2009), eprint 0909.1739.
 B. I. Abelev et al. (STAR Collaboration), Phys.Rev. C81, 054908 (2010), eprint 0909.1717.
 L. Adamczyk et al. (STAR Collaboration), Phys.Rev.Lett. 113, 052302 (2014), eprint 1404.1433.
 I. A. Shovkovy, Lect.Notes Phys. 871, 13 (2013), eprint 1207.5081.
 M. D’Elia, Lect.Notes Phys. 871, 181 (2013), eprint 1209.0374.
 K. Fukushima, M. Ruggieri, and R. Gatto, Phys.Rev. D81, 114031 (2010), eprint 1003.0047.
 M. N. Chernodub and A. S. Nedelin, Phys.Rev. D83, 105008 (2011), eprint 1102.0188.
 R. Gatto and M. Ruggieri, Phys.Rev. D85, 054013 (2012), eprint 1110.4904.
 A. A. Andrianov, D. Espriu, and X. Planells, Eur.Phys.J. C74, 2776 (2014), eprint 1310.4416.
 X. Planells, A. A. Andrianov, V. A. Andrianov, and D. Espriu, PoS QFTHEP2013, 049 (2013), eprint 1310.4434.
 J. Chao, P. Chu, and M. Huang, Phys.Rev. D88, 054009 (2013), eprint 1305.1100.
 L. Yu, H. Liu, and M. Huang, Phys.Rev. D90, 074009 (2014), eprint 1404.6969.
 A. Yamamoto, Phys.Rev.Lett. 107, 031601 (2011a), eprint 1105.0385.
 A. Yamamoto, Phys.Rev. D84, 114504 (2011b), eprint 1111.4681.
 V. Braguta, V. Goi, M. Ilgenfritz, A. Y. Kotov, A. Molochkov, et al., JETP Lett. 100, 547 (2015).
 V. Braguta, V. Goy, E. M. Ilgenfritz, A. Y. Kotov, A. Molochkov, et al. (2014), eprint 1411.5174.
 E.M. Ilgenfritz, M. Kalinowski, M. MüllerPreussker, B. Petersson, and A. Schreiber, Phys.Rev. D85, 114504 (2012), eprint 1203.3360.
 E.M. Ilgenfritz, M. MüllerPreussker, B. Petersson, and A. Schreiber, Phys.Rev. D89, 054512 (2014b), eprint 1310.7876.
 H. KlubergStern, A. Morel, O. Napoly, and B. Petersson, Nucl.Phys. B220, 447 (1983).
 I. Montvay and G. Münster, Quantum fields on a lattice (Cambridge University Press, 1994).
 P. Hasenfratz and F. Karsch, Phys.Lett. B125, 308 (1983).
 R. D. Pisarski and F. Wilczek, Phys.Rev. D29, 338 (1984).
 D. Wallace and R. Zia, Phys.Rev. B12, 5340 (1975).
 P. Hasenfratz and H. Leutwyler, Nucl.Phys. B343, 241 (1990).
 A. V. Smilga and J. Stern, Phys.Lett. B318, 531 (1993).
 A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. Ding, et al., Phys.Rev. D85, 054503 (2012), eprint 1111.1710.
 B. Wang, Y.L. Wang, Z.F. Cui, and H.S. Zong, Phys.Rev. D91, 034017 (2015).
 S.S. Xu, Z.F. Cui, B. Wang, Y.M. Shi, Y.C. Yang, et al., Phys.Rev. D91, 056003 (2015).
 M. Hanada and N. Yamamoto, PoS LATTICE2011, 221 (2011), eprint 1111.3391.
 S. Capitani, Phys.Rept. 382, 113 (2003), eprint heplat/0211036.