Lattice study of QCD at finite chiral density: topology and confinement

Abstract

In this paper we study the properties of QCD at nonzero chiral density $\rho _5$, which is introduced through chiral chemical potential $\mu _5$. The study is performed within lattice simulation of QCD with dynamical rooted staggered fermions. We first check that $\rho _5$ is generated at nonzero $\mu _5$ and in the chiral limit observe $\rho _5 \sim \varLambda _{QCD}^2 \mu _5$. We also test the possible connection between confinement and topological fluctuations. To this end, we measured the topological susceptibility $\chi _{\mathrm{top}}$ and string tension $\sigma$ for various values of $\mu _5$. We observed that string tension grows with $\mu _5$. It seems that topological susceptibility also rises with $\mu _5$, but to state this more reliably the uncertainties should be reduced. We believe that our results indicate possible connection between topological fluctuations and the strength of confinement.

1 Introduction

Quantum Chromodynamics (QCD) is believed to be the theory of strong interaction. While microscopic QCD Lagrangian is well known, the theory itself is extremely complicated and possesses a plenty of not fully understood nontrivial properties and phenomena. The most well-known examples include color confinement and chiral symmetry breaking.

One of the possible ways to shed light on these phenomena and their mechanism is to investigate QCD or QCD-like theories under extreme conditions. These extreme conditions include finite temperature studies [1,2,3,4], the influence of large magnetic field on QCD properties [5,6,7,8,9,10], QCD and QCD-like theories at finite baryon density [10,11,12,13,14,15,16,17,18] and QCD at finite isospin density [19,20,21,22].

Among others are the properties of QCD at nonzero chiral density. Systems with nonzero chiral density attract considerable attention because of unusual phenomena which take place in such systems. The renowned example of such phenomena is the chiral magnetic effect (CME) [23, 24], the appearance of electric current in chiral medium along applied magnetic field.

Nonzero chiral density can be generated in heavy ion collisions either due to sphaleron transitions in quark-gluon plasma [25,26,27] or due to the axial anomaly in parallel electric and magnetic fields [28]. Although its possible magnitude is under active debate, there are indications that it might reach relatively high values$\rho _5=8.85\times 10^{-1} \text {fm}^-3$ [29]. Nonzero chiral density is not conserved and depletes with time. Nonetheless, its relaxation time is expected to be relatively large, comparable with the life-time of Quark–Gluon Plasma [30, 31], or much larger [32]. It allows to conclude that nonzero chiral density can be treated as external parameter and makes it interesting to study the properties of Quantum Chromodynamics with nonzero chiral density. There are a lot of studies of QCD properties with chiral density which is introduced through nonzero chiral chemical potential [33,34,35,36,37,38,39,40,41,42,43,44].

One of the interesting questions which can be addressed is how the confinement and the chiral symmetry breaking in QCD are affected by nonzero chiral density. The influence of nonzero chiral chemical potential on the chiral symmetry breaking was considered in a number of theoretical papers [33,34,35,36,37, 41, 44] as well as in the lattice studies [42, 43]. Today it is clear that in any system the chiral chemical potential either creates or enhances the dynamical chiral symmetry breaking depending on the strength of interactions between constituents in the media. This phenomenon was called the chiral catalysis and the mechanism responsible for this phenomenon was first explained in [44]. The essence of this phenomenon is that nonzero chiral density generates additional fermionic states which take part in the formation of the chiral condensate.

In this paper we mainly address three questions. First, we show that introduction of nonzero $\mu _5$ to the system Hamiltonian leads to generation of nonzero chiral density $\rho _5$. We study its dependence within lattice simulation of QCD and compare the observed behavior with the existing models such as ChPT and NJL. Second, we study the influence of chiral density on topological structure of QCD. Finally, we study the confinement in QCD with nonzero $\mu _5$ and its connection to the topology of QCD. The possible link between these phenomena was introduced in [45, 46]. Namely, the authors suggested to modify the gluon propagator to have the form $G(p) = (p^2 + \chi _{\mathrm{top}} / p^2)^{-1}$ due to Veneziano ghosts tunnelling between different topological sectors of QCD. This form of gluon propagator implies maximum propagation range of order $\chi _{\mathrm{top}}^{-1/4}$ and suggests enhancement of confinement with the growth of topological susceptibility. To check this connection, we study the string tension $\sigma$ between heavy quark and antiquark at nonzero $\mu _5$ and its correlation with topological susceptibility $\chi _{\mathrm{top}}$.

It is well known that introduction of baryon chemical potential leads to the sign problem in SU(3) theory and spoils the LQCD simulations. On contrary, introduction of the chiral chemical potential does not lead to the sign problem [23], which allows us to carry out this study within lattice simulation of QCD.

This paper is organized as follows. In the next section we discuss the chiral density generated by nonzero chiral chemical potential in QCD. In the Sect. 3 we describe the details of our lattice simulation. Our results are presented in the Sect. 4. In the last section we discuss our results and draw the conclusions. In Appendix A we derive the chiral density for free naïve fermions and study divergences in the chiral density.

2 Nonzero chiral chemical potential in QCD

In this paper we are going to study the properties of QCD with nonzero chiral density $\rho _5=\bar{\psi }\gamma _4 \gamma _5 \psi$. It is well known that nonzero baryon density can be introduced to statistical system through modification of the Hamiltonian in the partition function $\hat{H} \rightarrow \hat{H} - \mu \int d^3 x \bar{\psi }\gamma _4 \psi$.¹ Similarly, one can modify the Hamiltonian by the term with chiral chemical potential $\mu _5$

$$\begin{aligned} \hat{H} \rightarrow \hat{H} - \mu _5 \int d^3 x \bar{\psi }\gamma _4 \gamma _5 \psi . \end{aligned}$$

(1)

We would like to stress that the chiral chemical potential is different to the baryon chemical potential since chiral density is not conserved. There are two operators resulting in the non-conservation of the chiral density

$$\begin{aligned} \frac{d}{dt} \int d^3 x \rho _5&= \frac{\alpha _s N_c}{4 \pi } \int d^3 x F^a_{\mu \nu } {\tilde{F}}^a_{\mu \nu } \nonumber \\&\quad +2 m \int d^3 x \bar{\psi }\gamma _5 \psi . \end{aligned}$$

(2)

The first operator $\sim F^a_{\mu \nu } {\tilde{F}}^a_{\mu \nu }$ is the anomalous contribution due to quantum corrections. The second operator $\sim m \bar{\psi }\gamma _5 \psi$ results from the equation of motion for massive fermions. Note that chirality is not well-defined for massive fermions due to the possible spin flipping process. The dynamical fermion mass generation $\sim \varLambda _{QCD}$ due to the chiral symmetry breaking can significantly increase the effect of spin flipping. Thus, the physical meaning of modification (1) should be discussed more carefully. It is clear that the $\rho _5$ operator becomes the true chiral density only in the massless limit $\rho _5\bigl |_{m\rightarrow 0}=(Q_R-Q_L)/V$. For massive quarks the meaning of the $\rho _5$ operator should be considered in more detail.

Chemical potential is usually introduced with respect to conserved charge. In our study we consider $\mu _5 \rho _5$ as the new term in the Hamiltonian and the conservation of $\rho _5$ is not required. We expect that the modification (1) leads to nonzero averaged value of the chiral density operator $\langle \rho _5 \rangle \ne 0$ even for nonzero quark mass. The situation with $\mu _5$ and $\rho _5$ is similar to the one with the fermion mass term $m \bar{\psi }\psi$. The conservation of the $\bar{\psi }\psi$ operator is not required and once this operator is introduced to the Hamiltonian it leads to the generation of nonzero condensate $\langle \bar{\psi }\psi \rangle \ne 0$. To show that it is very likely that non-zero $\mu _5$ will result in non-zero $\rho _5$ generation even at finite quark mass, let us consider various models of QCD.

First, in terms of fermionic spectrum, the modification of the Hamiltonian (1) modifies the dispersion relation $E^2(p)=(|\mathbf {p}| - s \mu _5)^2 + m^2$ [33], where $s=\pm 1$ is the fermion helicity. For $\mu _5>0$ this implies that at fixed momentum $|\mathbf {p}|$ the fermion with helicity $s=+1$ has smaller energy than the one with $s=-1$. In thermodynamic equilibrium there will be a larger number of fermions with helicity $s=+1$ than that with $s=-1$. So one can expect that the modification (1) leads to nonzero helicity even at nonzero quark mass.

We proceed with the consideration of $SU(N_c)$ QCD with finite chemical potential in the large $N_c$ limit $N_c \rightarrow \infty$. At low temperature T the chiral perturbation theory (ChPT) [47, 48] can be applied. In the leading order in $1/N_c$ there is no contribution of the anomalous term. Modification (1) only adds the flavour singlet axial current $A_{\mu } = \mu _5 \delta _{\mu 4} \hat{1}$ and the modification of the partition function due to the introduction of this axial current within ChPT reads

$$\begin{aligned} Z(\mu _5)=Z_{QCD}\times \exp {( \beta V N_f f_{\pi }^2 \mu _5^2)}, \quad \beta = \frac{1}{T}. \end{aligned}$$

(3)

From Eq. (3) it is seen that nonzero $\mu _5$ leads to the additional constant factor in the QCD partition function, i.e. it is not related to dynamical degrees of freedom of the ChPT. This is because the ChPT accounts the chiral symmetry breaking in QCD but it does not provide its mechanism. In more complicated models [49, 50] which consider the chiral symmetry breaking mechanism, the $\mu _5$ couples to the scalar $\sigma$ and $\eta '$ fields, which leads to enhancement of the chiral symmetry breaking with $\mu _5$ and the chiral catalysis phenomenon in QCD [44].

From (3) for two flavor QCD one has

$$\begin{aligned} \langle \rho _5 \rangle = \frac{1}{\beta V} \frac{\partial \log Z(\mu _5)}{\partial \mu _5} = 4 f_{\pi }^2 \mu _5. \end{aligned}$$

(4)

So, one can see that the modification of the Hamiltonian (1) indeed leads to nonzero $\langle \rho _5 \rangle$ in the limit $N_c \rightarrow \infty$ even at nonzero quark mass.

Similar study can be carried out in the Nambu–Jona–Lasinio (NJL) model [51], which successfully describes low energy phenomenology of QCD. Since the NJL model is usually studied within the saddle point approximation, a lot of results are obtained within the $N_c \rightarrow \infty$ assumption. Within the NJL model, the chiral symmetry breaking leads to the generation of the dynamical quark mass $m\sim \varLambda _{QCD}$. The calculation of the chiral density with the quark mass $m\sim \varLambda _{QCD}$ gives $\langle \rho _5 \rangle \sim \varLambda _{QCD}^2 \mu _5$ (see Appendix A). This is another argument in favor of the hypothesis that non-zero chiral density can be generated at finite mass $\langle \rho _5 \rangle \sim \varLambda _{QCD}^2 \mu _5 \ne 0$.

On the one hand the approximation $N_c \rightarrow \infty$ works quite well for real QCD. So, one might expect that $\langle \rho _5 \rangle \sim \varLambda _{QCD}^2 \mu _5 \ne 0$ for $N_c=3$. However, the anomaly contribution which appears in higher orders $\sim 1/N_c$–corrections can modify the $N_c \rightarrow \infty$ result for the chiral density. To clarify the $N_c = 3$ behavior in Sect. 4 we conduct lattice study of chiral density $\rho _5$ at nonzero chiral chemical potential.

3 Lattice setup

In this paper we are going to study QCD with two flavours and nonzero chiral chemical potential. To this end we perform lattice simulations with the SU(3) gauge group and employed the tree level improved Symanzik gauge action [52, 53]. For the fermionic part of the action we used staggered fermions with the action [43]

$$\begin{aligned} \begin{aligned} S_f&=ma\sum _x {\bar{\psi }_x}\psi _x \\&\quad +\frac{1}{2}\sum _{x\mu } \eta _{\mu }(x)({\bar{\psi }_{x+\mu }U_{\mu }(x)}\psi _x -{\bar{\psi }_x}U^{\dag }_{\mu }(x) \psi _{x+\mu }) \\&\quad + \frac{1}{2}\mu _5a\sum _x s(x)({\bar{\psi }}_{x+\delta }{\bar{U}}_{x+\delta , x}\psi _x -{\bar{\psi }}_{x}{\bar{U}}_{x+\delta , x}^{\dag }\psi _{x+\delta }), \end{aligned} \end{aligned}$$

(5)

where the $\eta _{\mu }(x)$ are the standard staggered phase factors: $\eta _1(x)=1,\eta _{\mu }(x)=(-1)^{x_1+\ldots +x_{\mu -1}}$ for $\mu =2,3,4$². The lattice spacing is denoted by a, the bare fermion mass by m, and $\mu _5$ is the chiral chemical potential. In the chirality breaking term $s(x)=(-1)^{x_2}$, $\delta =(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+\delta$,

$$\begin{aligned} \begin{aligned} {\bar{U}}_{x+\delta ,x}=\frac{1}{6}\sum \limits _{i,j,k= \text {perm}(1,2,3)}U_i(x+e_j+e_k)U_j(x+e_k)U_k(x) \end{aligned} \end{aligned}$$

(6)

is symmetrized over the 6 shortest paths between these sites. In the partition function, after integrating out fermions, one obtains the corresponding fermionic determinant. In order to obtain two flavours in the continuum limit we apply the rooting procedure. Simulations were performed using Rational Hybrid Monte Carlo algorithm.

In the continuum limit and after rooting procedure our lattice action can be rewritten in the Dirac spinor-flavor basis [56, 57] as follows

$$\begin{aligned} \begin{aligned} S_f \rightarrow S^{(cont)}_f&= \int d^4x \sum _{i=1}^2 \bar{q_i} (\partial _{\mu }\gamma _{\mu } +igA_{\mu }\gamma _{\mu }\\&\qquad +m+\mu _5\gamma _5\gamma _4)q_i. \end{aligned} \end{aligned}$$

(7)

We would like to emphasize that the chiral chemical potential introduced in Eq. (5) corresponds to the taste-singlet operator $\gamma _5\gamma _4\otimes \varvec{1}$ in the continuum limit.

Table 1 Lattice parameters and typical number of configurations used in the simulations

It should be also noted here that the baryonic chemical potential [58] and the chiral chemical potential as in [59], are introduced to the action as the modification of the temporal links by the corresponding exponential factors in order to eliminate chemical-potential dependent quadratic divergences. For staggered fermions with the baryonic chemical potential this modification can be performed. However, in the case of $\mu _5$ this method would lead to a highly non-local action [59]. Therefore, we introduce $\mu _5$ in Eq. (5) in the additive way similarly to the mass term. It is known that the additive introduction of the chemical potential might lead to additional divergences in observables. In this paper we perform lattice measurement of chiral density and gluonic observables: the topological charge, the topological susceptibility and the string tension. In what follows we account ultraviolet divergences in the chiral density. We also believe that there are no additional divergences due to chiral chemical potential in gluon observables, because the chiral chemical potential term can be considered as some vertex with coupling constant of dimension of energy. It is known that the inclusion of such vertex to Feynman diagrams reduces the power of ultraviolet divergences. Since the fermion loops in QCD diverge as powers of $\log a$, the chiral chemical potential does not give rise to additional divergences. The ultraviolet divergences in QCD with chiral chemical potential are also discussed in [42, 43].

The physical lattice spacing a was determined from setting Sommer parameter $r_0$ [60] to its physical values $r_0 = 0.468(4)$ fm [61]. Simulation for scale setting were performed with the lattice size $24^4$, $\mu _5 = 0$ and fixed $ma = 0.01$. Since the Sommer scale very mildly depends on the quark mass [62], the physical units are almost independent from the quark mass. Notice also that as was shown in papers [42, 43] nonzero $\mu _5$ does not affect to the scale setting procedure.

In the calculation we employed three different lattices with different lattice spacings to keep the physical volume fixed at approximately $1.7\,\hbox {fm}^3$: $14^4$ with $a = 0.128(3)\,\hbox {fm}$ ($\beta = 3.9$), $16^4$ with $a = 0.1054(11)\,\hbox {fm}$ ($\beta = 4.0$) and $20^4$ with $a = 0.0856(14)\,\hbox {fm}$ ($\beta = 4.1$). To investigate chiral properties for each of the listed lattices three values of pion mass were considered: $m_\pi = 563,\,762,\,910\,\hbox {MeV}$. We summarize lattice parameters of the simulations and number of configurations in Table 1, during the generation of ensembles only configurations separated by 4 MD-trajectories were saved, and the typical length of Molecular Dynamics trajectory was close to 1. We note that the simulations performed in this paper indicate that the required simulation time grows with the chiral chemical potential. Lattice simulations at the largest values of chiral chemical potential are numerically very expensive. Statistical error analysis was performed with binned jack-knife method and typical bin size O(100) configurations. We checked, that this bin size is big enough to incorporate autocorrelation time for all observables considered in this study.

4 Results of the calculation

4.1 The chiral density

In this section we perform lattice measurement of the chiral density $\rho _5$ for all lattice spacings and pion masses under study. The chiral densities as a function of the chiral chemical potential for different pion masses and $a=0.105\,\hbox {fm}$ are shown in Fig. 1. The chiral densities for other lattice spacings look similar. For this reason we do not show them. From the upper panel of Fig. 1 one sees that the data are well described by the linear dependence. It turns out that the coefficient of this linear dependence can be mostly attributed to the ultraviolet divergence in $\rho _5$. However, our data are rather accurate. Typical uncertainty of the calculation is $\sim 0.1 \%$, for this reason we can extract the sub-leading terms on the background of leading ultraviolet divergence.

Fig. 1

The chiral densities as a function of the chiral chemical potential for different pion masses and $a=0.105\,\hbox {fm}$, unrenormalized (upper panel) and after the subtraction of divergent parts (lower panel)

To proceed we need to know the structure of the divergences in $\rho _5$. In Appendix A the study of the ultraviolet divergences in $\rho _5$ for free naïve fermions is presented. In particular, it is shown that there are two ultraviolet divergences in the term linear in $\mu _5$. The leading divergence is quadratic and the next-to-leading divergence is logarithmic. Additionally, the linear in $\mu _5$ term contains finite contribution. Finally higher terms in $\mu _5$ expansion do not contain ultraviolet divergences.

In this paper we are going to use the following anzats for $\rho _5$ which accounts for the results obtained in Appendix A:

$$\begin{aligned} a^3 \rho _5&= E (a \mu _5)^3 + (A + a^2 B + C_1 (m a)^2 + \nonumber \\&\quad +D (ma)^2 \log (ma)^2 + F a^2 (ma)^2 + X a^4) \times (a \mu _5). \end{aligned}$$

(8)

This fit gives decent description of the data $\chi ^2 / \hbox {ndof} \sim 3$. Since the measurements of $\rho _5$ are quite accurate (at some points the error is only 0.05 %), we are able to fix all the parameters with the error of not worse than 15 %. Removing of any of the terms in (8) leads to significant growth of $\chi ^2 / \hbox {ndof}.$ However, adding higher powers of ma and a to the fit does not improve the quality.

Among the coefficients in Eq. (8), A and D correspond to the divergencies. In the lower panel of Fig. 1 we present the results for the chiral density $\rho _5$ after the subtraction of these divergent parts.

It is important to notice that the coefficient B from (8) is non-zero and $B = (340(10)\,\hbox {MeV})^2$. This coefficient parameterizes the chiral density in the continuum and in the chiral limits. For this reason we can state that $B \sim \varLambda _{QCD}^2$ or $\rho _5 \sim \varLambda _{QCD}^2 \mu _5$. Notice, however, that it is not possible to write exactly $\rho _5 = B\mu _5$ since the multiplicative renormalization of $\rho _5$ might be important but goes beyond the scope of this paper. To summarize, the results of this section allow us to state that finite $\mu _5$ generates nonzero chiral density $\rho _5 \sim \varLambda _{QCD}^2 \mu _5 + O(\mu _5^3)$. This finding up to the proportionality coefficient is in striking agreement with the recent calculations performed within effective models [63, 64], although in these papers authors get somewhat smaller values of the coefficient B. We believe, that this may be attributed to possible multiplicative renormalization of the results. Also, we would like to point our that the fit parameters obtained within Eq. 8 are in quantitative agreement with the free naive fermion computation performed in Appendix A. For instance, the value obtained within the fitting $A = -0.3006(3)$ might be just a multiplicative renormalization of $-0.4648$ from Eq. 18. Similarly, the leading $\mu _5^3$ coefficient $E = 0.0967(2)$ is suggestively the renormalized $J_2$ from Eq. 18. We thus conclude that the results obtained within fully interacting theory qualitatively agree with naive computation in free theory.

4.2 The topological charge and topological susceptibility

Our next task is to study how nonzero chiral chemical potential influences the topological properties of QCD. To this end we measure the topological charge $\left\langle Q \right\rangle$ and the topological susceptibility $\left\langle Q^2 \right\rangle$ for different values of the chiral chemical potential under study.

Our measurement of the topological charge and the topological susceptibility mainly follows [65]. We smoothen each configuration using the Gradient Flow [66, 67]. Topological charge is measured on the smoothened configurations

$$\begin{aligned} Q_L=-\frac{1}{512\pi ^2}\sum _x\sum _{\mu \nu \rho \sigma =\pm 1}^{\pm 4}\tilde{\epsilon }_{\mu \nu \rho \sigma }{\text {Tr}}U_{\mu \nu }(x)U_{\rho \sigma }(x)\,, \end{aligned}$$

(9)

where $U_{\mu \nu }(x)$ is the plaquette at the point x in directions $\mu$ and $\nu$. In order to reduce the lattice artifacts we used the following estimators of the topological charge Q:

$$\begin{aligned} Q={\text {round}}\left( \alpha Q_L\right) , \end{aligned}$$

(10)

where ${\text {round}}$ gives the closest integer to its argument and the factor $\alpha$ is chosen in such a way that it minimizes

$$\begin{aligned} \langle \left( \alpha Q_L-{\text {round}}\left( \alpha Q_L\right) \right) ^2 \rangle \,. \end{aligned}$$

(11)

Fig. 2

The topological susceptibility as the function of the chiral chemical potential for different pion masses in the continuum limit performed using the non-singlet pion mass correction of $\chi ^{1/4}(a)$ procedure described in Appendix B

In other words, we rescale our definition of the topological charge $Q_L$ so that its peaks become closer to integer values and then round the result to this integer value. The topological susceptibility is then defined as

$$\begin{aligned} \chi _{\mathrm{top}}=\frac{\langle Q^2 \rangle }{V_4}\,, \end{aligned}$$

(12)

where $V_4$ is the four-dimensional volume of the lattice. We have found that for Gradient Flow times $t/a^2>3.0$ the dependence of the topological susceptibility $\chi _{\mathrm{top}}$ on the value of Gradient Flow time exhibits a plateau with almost no dependence on the value of t. The value at this plateau was taken as a final estimation for the topological susceptibility $\chi _{\mathrm{top}}$. In the Appendix B we study systematic effects, such as topological freezing and discretization errors of the topological susceptibility.

Our results for the topological properties of QCD are the following. The topological charge is zero within the uncertainty of the calculations for all pion masses, lattice spacings and chiral chemical potentials under study.

In Fig. 2 we show the topological susceptibility as the function of chiral chemical potential for different pion masses in the continuum limit performed according to the non-Goldstone correction procedure described in Appendix B. Note that errorbars contain both statistical and systematic uncertainty, the detailed description of error estimation can be found in Appendix B. The chiral chemical potential seems to enhance the topological fluctuations in QCD for all pion masses, but our data need to be improved to make more strict statements.

A possible explanation of this behaviour is the following. As we know nonzero chiral chemical potential leads to generation of nonzero chiral charge in the system with some average value $Q_5$. Due to the anomaly this chiral charge can annihilate to gluon configurations with nonzero Chern–Simons number, which is compensated by an inverse process: creation of the chiral charge from gluon background with Chern-Simons number. In the thermodynamic equilibrium these processes compensate each other leading to some fixed average value of the chiral density. Notice also that both processes result from the chiral anomaly. Further let us consider the process of annihilation of the chiral charge as a number of elementary processes in which one quark and one antiquark annihilate to gluon configuration with nonzero Chern–Simons number. It is reasonable to assume that the larger the chiral charge the larger the number of elementary annihilations per time unit in the system. In other words the larger the chiral charge $Q_5$ the larger average $\langle d Q_5 / dt \rangle _{annihilation}$ for the annihilation processes. Notice that is completely compensated by the inverse process leading to the total $\langle d Q_5 / dt \rangle =0$ From this picture one can expect that the larger the chiral charge the larger the topological fluctuations in the system under investigation. Our results imply that $\mu _5$ is the parameter which allows influencing the topological sector of QCD through the anomaly equation.

4.3 The string tension

In order to study how nonzero chiral density influences the confinement properties of QCD we calculated the interaction potential of static charges through the measurement of Wilson loops. To obtain reasonable signal-to-noise ratio for Wilson loops the smearing techniques were employed. One step of the hypercubic blocking [68] with parameters $\alpha = (1.0,\,1.0,\,0.5)$ [69] was performed for the temporal links only, followed by 24 steps of the APE smearing [70] with $\alpha _{APE} = 0.165$.

The quark-antiquark interaction potential is related to Wilson loops as

$$\begin{aligned} V(R) = \lim _{t \rightarrow \infty } \text {log} \Bigl [ \frac{\left\langle W(R,t) \right\rangle }{\left\langle W(R,t+1) \right\rangle } \Bigr ]\,. \end{aligned}$$

(13)

This logarithm exhibits a clear plateau at large times $t/a\in [5; 9]$, its height was extracted as V(R).

String tension $\sigma$ was obtained from fitting of the potential in the range $R \in [3.5 a ; L_s / 2]$ by the Cornell fit

$$\begin{aligned} V(R) = A - \alpha / R + \sigma R\,. \end{aligned}$$

(14)

This fit provides $\chi ^2 / \hbox {ndof} \lesssim 1$ for all values of chiral chemical potentials. To estimate systematic uncertainty the left fitting range was varied in the interval [3a; 4a] and the produced small change of $\sim 0.5\%$ in the string tension was added to the statistical error. Change of the right boundary of R in the fit does not alter the results in a noticeable way. Statistical errors for the fit parameters were estimated with the jackknife method.

Fig. 3

The ratio of the string tension $\sigma$ to the string tension at zero chiral chemical potential $\sigma _0$. Points for different pion masses are slightly shifted in horizontal axis for better visibility

It is worth to note, that the Wilson loop corresponds to an operator, that creates the static color sources and a string between them, and this operator has a small overlap with the state, corresponding to a broken string [71], thus the string breaking phenomenon can not be observed from the Wilson loops. On the other hand, for finite lattice due to the p.b.c. in spatial directions the maximal achievable separation between q and $\bar{q}$ is $L_s/2$, which in our case corresponds to 0.85 fm ($L_s \approx 1.7$ fm, see Table 1). The string breaking for the physical pion mass appears near 1 fm, and in our study due to the heavier pions this should occur ever at larger $q \bar{q}$ separation. Thus extraction of quark-antiquark interaction potentials from Wilson loops does not lead to any problems as far as we investigate V(R) at distances, which are smaller than the string breaking distance. The same argument applies to the choice of Cornell potential for V(R) fitting.

We do not observe any significant dependence of the string tension on the lattice spacing a. Thus, we perform the constant fit of the string tension versus a to average over different lattice spacings. The quality of such fit is good, $\chi ^2 / \hbox {ndof} < 1$. The ratio of the string tension $\sigma$ (extrapolated to continuum in the above described way) to the string tension at zero chiral chemical potential $\sigma _0$ is presented in Fig. 3. It is seen from Fig. 3 that the string tension rises with the chiral chemical potential i.e. with the chiral density.

4.4 Topological fluctuations and confinement

The phenomenon of the QCD confinement is not well understood on the present day. However, papers [45, 46] have established a possible link between confinement properties and QCD topology. In their setup, gluon propagator is modified by the interaction with Veneziano ghosts tunneling between different topological sectors. The gluon propagator then reads $G(p) = (p^2 + \chi _{\mathrm{top}} / p^2)^{-1}$, where $\chi _{\mathrm{top}}$ is the topological susceptibility. The propagator has only complex poles $p^2 = \pm i \chi _{\mathrm{top}}^{1/2}$, thus gluons cannot propagate as free particles. The typical range of gluon propagation decreases as $\chi _{\mathrm{top}}^{-1/4}$ with the growth of topological susceptibility.

Our results might be in agreement with this picture of confinement. As one can observe from Sect. 4.2, it seems that the topological susceptibility is enhanced by $\mu _5$. We also observe ( see Sect. 4.3) that the confining properties, namely the string tension, are also enhanced by $\mu _5$. Notice, however, that to make our statement more reliable one should reduce the uncertainties in the topological susceptibility.

5 Conclusion and discussion

In this paper we studied the properties of QCD at nonzero chiral density $\rho _5$, which is introduced through the chiral chemical potential $\mu _5$. Contrary to the baryon chemical potential introduction of the chiral chemical potential does not lead to the sign problem. For this reason our study of QCD with nonzero chemical potential can be performed within lattice simulation. In the simulations we employed the tree level improved Symanzik gauge action and rooted staggered fermions which in the continuum limit correspond to $N_f=2$ dynamical quarks.

In the calculation we employed three different lattices with different lattice spacings to keep the physical volume fixed at approximately $1.7\,\hbox {fm}^3$: $14^4$ with $a = 0.128(3)\,\hbox {fm}$ ($\beta = 3.9$), $16^4$ with $a = 0.1054(11)\,\hbox {fm}$ ($\beta = 4.0$) and $20^4$ with $a = 0.0856(14)\,\hbox {fm}$ ($\beta = 4.1$). To investigate the chiral properties for each of the listed lattices three values of pion mass were considered, $m_\pi = 563,\,762,\,910\,\hbox {MeV}$.

The first observable considered in this paper is the chiral density. We found that nonzero chiral chemical potential leads to generation of nonzero chiral density in QCD. Our lattice results support ChPT formula for the chiral density $\rho _5 \sim \varLambda _{QCD}^2 \mu _5$.

The next question is the influence of nonzero chiral chemical potential on the topological properties of QCD. To address this question we measured the topological charge and the topological susceptibility for various values of $\mu _5$. We found that the topological charge is zero for all values of the chiral chemical potential under investigation. It seems that topological susceptibility rises with $\mu _5$.

The last observable studied in this paper is the string tension. We calculated the static potential from Wilson loops and determined the string tension for all values of chiral chemical potentials at lattice parameters studied. We found that the string tension rises with rising chiral chemical potential.

It would be interesting to understand the mechanism how confinement in QCD is enhanced by nonzero chiral chemical potential. One possible explanation can be based on the results of [45, 46], where the authors considered the gluon propagator, modified due to Veneziano ghosts tunneling between different topological sectors, making gluons confined at typical distances $\sim \chi _{\mathrm{top}}^{-1/4}$ where $\chi _{\mathrm{top}}$ is the topological susceptibility. Our results might confirm this picture of confinement. The topological susceptibility seems to be enhanced by $\mu _5$. The string tension, describing the confining properties, also grows with $\mu _5$. However, to make our statement more reliable, the uncertainties in the topological susceptibility should be reduced.

Another possible explanation is that the gluon fields generated in the system due to fluctuations of $\rho _5$ might have nontrivial properties which give rise to the confinement. In particular, if the gluon fields are self-dual due to the $\rho _5$ fluctuations they might enhance the confinement [72,73,74,75]. Unfortunately, quite large uncertainties of the calculation do not allow us to draw any strong conclusion about the origin of confinement enhancement with $\mu _5$. This question including the mechanism of self-dual gluon fields is the subject for further research.

Appendices

Ultraviolet divergences in the chiral density for free naïve fermions

To get an idea about the ultraviolet divergences in the chiral density at nonzero chiral chemical potential in this section we are going to derive the chiral density for free naïve fermions. The fermion propagator including the chiral chemical potential for naïve lattice fermions can be written in the following form:

$$\begin{aligned} S^{\alpha \beta }(x,y)= & {} \frac{\delta ^{\alpha \beta }}{L_t L_s^3} \sum _{ \{p\} } \sum _s e^{i p (x-y)}\nonumber \\&\frac{-i \sum _{\mu } \gamma _{\mu } \sin (p_{\mu }) + ma + (\mu _5a) \gamma _4 \gamma _5 }{\sin ^2 (p_4) + (|p|-s (\mu _5a))^2 +(ma)^2} \times P(s) ,\nonumber \\ P(s)= & {} \frac{1}{2} \biggl ( 1 - i s \sum _i \frac{\gamma _{i} \sin (p_{i}) }{|p|} \gamma _0 \gamma _5 \biggr ),~~~i=1,2,3, \nonumber \\ |p|^2= & {} \sin ^2 (p_1) + \sin ^2 (p_2) + \sin ^2 (p_3), \nonumber \\ p_i= & {} \frac{2 \pi }{L_s} n_i,~~ i=1,2,3,~~ n_i=0,...,L_s-1, \nonumber \\ p_4= & {} \frac{2 \pi }{L_t} n_4 + \frac{\pi }{L_t},~~~ n_4=0,...,L_t-1. \end{aligned}$$

(15)

Here m and $\mu _5$ are mass and chiral chemical potential in physical units, $\alpha , \beta$ are color indices, the sum is taken over all possible values of $(n_1, n_2, n_3, n_4)$, $s= \pm 1$ and a is a lattice spacing.

In the limit $L_s, L_t \rightarrow \infty$ the chiral density in lattice units for two fermion flavours can be written as

$$\begin{aligned} \langle \bar{\psi }\gamma _4 \gamma _5 \psi \rangle _{lat}= & {} - \frac{3 N_f}{16} Sp\bigl [ \gamma _4 \gamma _5 S(x,x) \bigr ] \nonumber \\= & {} -\frac{3}{4} \sum _{s=\pm 1}\nonumber \\&\int \frac{d^4 p}{(2\pi )^4} \frac{s |p| - \mu _5a }{\sin ^2 p_4 {+} (|p|-s(\mu _5a))^2 {+} (ma)^2 }\nonumber \\ \end{aligned}$$

(16)

It should be noted here that the factor 3 in the first equality is due to the sum over the fermion colors. Now let us expand the chiral density (16) in powers of the chiral chemical potential. It turns out that it is sufficient to keep only two terms: $\sim \mu _5$ and $\sim \mu _5^3$. Higher order terms in this expansion do not contain ultraviolet divergences. The calculation of the integrals which appear in this expansion is rather cumbersome but straightforward. For this reason we don’t show the details of the calculation. We would like only to mention that the integrals which appear in the expansion of (16) in $\mu _5$ can be found in [76]. The resulting expression for the chiral density $\rho _5$ in physical units can be written in the following form

$$\begin{aligned} \rho _5= & {} \frac{1}{a^3} \langle \bar{\psi }\gamma _4 \gamma _5 \psi \rangle _{lat} = {\mu _5} J_1 + {\mu _5^3} J_2 + O(\mu _5^5) \nonumber \\ J_1= & {} -0.464800 \frac{1}{a^2} - \frac{3}{\pi ^2} m^2 \log (ma)^2 + 0.807241 m^2 \nonumber \\ J_2= & {} 0.242419 \end{aligned}$$

(17)

Now few comments are in order.

• From Eq. (17) we notice that there are two divergences in the linear in the $\mu _5$ term. The leading divergence is quadratic and the next-to-leading divergence is logarithmic.

• In addition to the divergences the linear in the $\mu _5$ term contains finite contribution which is proportional to the fermion mass in the second power $\rho _5 \sim m^2 \mu _5$. Now recall that in Nambu–Jona–Lasinio model [51], which successfully describes low energy phenomenology of QCD, the chiral symmetry breaking leads to generation of the dynamical fermion mass $m \sim \varLambda _{QCD}$. For this reason one can expect that due to chiral symmetry breaking in QCD the renormalized $\rho _5 \sim \varLambda _{QCD}^2 \mu _5$. ChPT confirms this statement (see Sect. 2).

• Notice also that the logarithmic ultraviolet divergence is also possible in the $\mu _5^3$ term. However, final result does not contain the logarithmic ultraviolet divergence.

• If one takes the chiral limit in formula (5), it is possible to get rid of logarithmic divergence as well as final term proportional to $m^2$. Unfortunately it is not possible to get rid of the $1/a^2$ divergence which results from the additive way of introducing the chiral chemical potential. In addition to the divergence the additive chemical potential modifies the coefficient in front of the $\mu _5^3$ contribution. For free chiral fermions this coefficient is determined by Fermi distribution and for three colors and two flavours it is $\rho _5(\mu _5)=2/\pi ^2 \cdot \mu _5^3\simeq 0.202642 \cdot \mu _5^3$ ([23]). Comparing this value with the $J_2$ in formula (17) it is seen that lattice artificial contribution is rather small but it is present. In this paper we concentrate on the linear in $\mu _5$ term, thus the fact that the coefficient of the $\mu _5^3$ term is modified by lattice artifacts does not affect the results of this paper.

Topological susceptibility in the continuum limit

Fig. 4

Monte Carlo history of topological charge(upper panel) and histogram of topological charge distribution(lower panel). Lattice size is $20^4$, lattice step $a=0.0856$fm, quark mass is $ma=0.00658$ and chiral chemical potential $\mu 5=0$

It is known that measuring the topological susceptibility on the lattice is a challenging task. First, lattice simulations suffer from topological freezing: transitions between sectors with different topological charge values can be suppressed and system can be stuck in one topological sector. We have checked, that in our simulations this was not the case. In Fig. 4 we present the Monte Carlo history of topological charge versus the configuration number and also the histogram of the topological charge per configuration for the lattice size $20^4$, the lattice spacing $a=0.0856$ fm, bare quark mass is $ma=0.00658$, chiral chemical potential $\mu _5=0$. MC histories for other lattice spacings (see Table 1) look the same, thus we do not present them here. One may clearly see, that in our simulations transitions between various topological sectors occur often enough, also the histogram for topological charge is broad.

Another problem is that topological susceptibility $\chi$ suffers from large discretization errors. To check the dependency of our data on a finite lattice step, we performed continuum extrapolation for $\chi$, using two various procedures, described in [77, 78]. In the case of zero chiral chemical potential $\mu _5=0$ we compare our results with the predictions of the Chiral Perturbation Theory (ChPT).

First of all, in [78] it is noted that the dominant source of lattice artefacts in $\chi ^{1/4}(a)$ is the chiral symmetry breaking present at finite lattice spacing in the staggered discretization. The dependence of $\chi ^{1/4}(a)$ on the lattice spacing can be significantly reduced if instead of $\chi ^{1/4}(a)$ the quantity

$$\begin{aligned} \chi ^{1/4}_{tc}(a) = \frac{m_{\pi }}{m_{\mathrm{ngb}}(a)} \chi ^{1/4}(a) \end{aligned}$$

(18)

is considered. Here $m_{\mathrm{ngb}}(a)$ is the mass of one of the non-Goldstone pions, i.e. of a state that becomes massless in the chiral limit only if the continuum limit is taken. Clearly, $m_{\mathrm{ngb}} \rightarrow m_{\pi }$ as $a \rightarrow 0$, so in the continuum limit $\chi ^{1/4}_{tc}(a) \rightarrow \chi ^{1/4}(a).$ Following [78], as the $m_{\mathrm{ngb}}(a)$ we used the state with the taste structure $\gamma _i \gamma _{\mu }$, which mass is close to the root mean square of all other taste masses.

Fig. 5

The topological susceptibility $\chi ^{1/4}(a)$ naïve fit (linear hatching) and the topological susceptibility NGB–corrected fit in accordance with 18$\chi ^{1/4}_{tc}(a)$ (cross hatching). The continuum extrapolation is done using the fit $\sim A + a^2 B$. Chiral chemical potential is zero $\mu _5=0$

In Fig. 5 we show a simple naïve fit $\sim A + a^2 B$ of the data for the topological susceptibility $\chi ^{1/4}(a)$ (linear hatching) and the NGB–corrected topological susceptibility $\chi ^{1/4}_{tc}(a)$ (cross hatching) for all pion masses

$563\,\hbox {MeV},\,762\,\hbox {MeV},\,910\,\hbox {MeV}$ and zero chiral chemical potential $\mu _5=0$. In both cases the data are described by this simple quadratic function $\sim A + a^2 B$ with $\chi ^2 / \hbox {ndof} < 1$. Note also that the results of the naïve fit of $\chi ^{1/4}(a)$ agree with the NGB–corrected fit $\chi ^{1/4}_{tc}(a)$ in the continuum limit. The same agreement was observed in [78]. Since both naïve and NGB-corrected extrapolations give similar results, we expect that our continuum extrapolation is under control, although data at smaller lattice step would be desirable to make more precise estimations of systematic uncertainties of the extrapolation.

Another two-step way to check the discretization errors and also compare result with the ChPT is discussed in [77]. To do so, we fit the inverse topological susceptibility versus squared pion mass $m_{\pi }^2$ for every fixed lattice spacing a and fixed chiral chemical potential $\mu _5$ with the ChPT–motivated ansatz

$$\begin{aligned} \frac{1}{\chi (a)} = \frac{A(a)}{m_{\pi }^2} + B(a). \end{aligned}$$

(19)

When the fitting parameters A(a) and B(a) are obtained, we perform the continuum extrapolation using the simple square ansatz:

$$\begin{aligned} \begin{aligned} A(a)&= A_0 + a^2 A_1,\\ B(a)&= B_0 + a^2 B_1. \end{aligned} \end{aligned}$$

(20)

Fig. 6

Topological susceptibility at zero chiral chemical potential for all lattice spacings plotted as the function of squared pion mass $m_{\pi }^2$. The orange cross-hatched region corresponds to continuum extrapolation result in accordance with the Eq. 19. Violet left-pointing triangles and pink triangles show the result of the naïve continuum extrapolation and NGB–corrected extrapolation from Fig. 5 respectively

Both fitting stages provide a good description of the data with $\chi ^2 / \hbox {ndof} \approx 1$ for all lattice spacings and both A(a) and B(a). We then interpret the function $\chi ^{-1}_0 = A_0 / m_{\pi }^2 + B_0$ as the continuum extrapolation of the inverse topological susceptibility. In Fig. 6 we show our data for all lattice spacings and pion masses together with the continuum extrapolation result (yellow cross-hatched region) obtained within the two-step procedure discussed in [77]. We also plot violet left-pointing triangles and pink triangles that the result of the naïve continuum extrapolation and NGB–corrected extrapolation from Fig. 5 respectively. The pictures for $\mu _5 > 0$ look the same and thus we do not show them here.

Note that the results of all three extrapolation procedures (naïve, NGB-corrected and two-step) are in agreement with each other. Finally, the ChPT states that $A_0 = 4 / f_{\pi }^2$. The extracted from two-step extrapolation value of the pion decay constant is $f_{\pi }=92(9)$ MeV, which agrees with the physical $f_{\pi }^{ph}\approx 93$ MeV. These results suggest that continuum extrapolation employed in this paper is under control, although simulations at smaller lattice step a are required to make more precise estimations of systematic uncertainty.

As the final result for topological susceptibility we take the values obtained by NGB–corrected procedure. We consider the difference between the NGB–corrected and two-step methods as an estimation of systematic uncertainty of the extrapolation. These results are presented in Fig. 2, where both statistical and systemic uncertainties are taken into account in the errorbars.