Complex Langevin simulation of QCD at finite density and low temperature using the deformation technique

We study QCD at finite density and low temperature by using the complex Langevin method. We employ the gauge cooling to control the unitarity norm and introduce a deformation parameter in the Dirac operator to avoid the singular-drift problem. The reliability of the obtained results are judged by the probability distribution of the magnitude of the drift term. By making extrapolations with respect to the deformation parameter using only the reliable results, we obtain results for the original system. We perform simulations on a $4^3\times 3$ lattice and show that our method works well even in the region where the reweighting method fails due to the severe sign problem. As a result we observe a delayed onset of the baryon number density as compared with the phase-quenched model, which is a clear sign of the Silver Blaze phenomenon.


Introduction
Investigating the QCD phase diagram at finite density and temperature is one of the central issues in modern theoretical physics. In particular, QCD at high density and low temperature is interesting due to its relevance to the state of dense matter that is considered to be realized in astronomical objects like neutron stars. It is known, however, that QCD in this parameter region suffers from a severe sign problem, which makes it practically inaccessible by conventional lattice QCD simulations based on the importance sampling.
One promising approach to the sign problem that has attracted much attention in recent years is the complex Langevin method (CLM) [1,2]. The CLM is based on the stochastic time evolution for the complexified dynamical variables using the complex Langevin equation. Since this method does not rely on the probabilistic interpretation of the Boltzmann weight in the path integral formulation of the original theory, one has a chance to solve the sign problem. It is known, however, that the equivalence between the CLM and the original path integral does not always hold [3,4]. In order to prove this equivalence, we actually need to satisfy the condition that the probability distribution of the magnitude of drift term decays exponentially or faster [5]. If this condition is met, one can show that the CLM gives the correct result. If the drift distribution falls off more slowly than exponential, the proof of the equivalence does not go through, and the results obtained by the CLM can be simply wrong. Thus, in applying the CLM, it is extremely important to judge the reliability of the obtained results by monitoring the asymptotic behavior of the drift distribution.
In fact, it is known that there are two causes for the slow decay of the drift distribution. One is the excursion problem, where the complexified dynamical variables run deeply into the imaginary direction. The other is the singular-drift problem, where the probability distribution of the configurations does not vanish at the singularity of the drift term. In order to overcome these problems and to enlarge the applicability region of the CLM, various techniques have been developed [6][7][8][9]. For instance, the gauge cooling [6] is a technique for avoiding the excursion problem by making use of the complexified gauge transformation. (It can be used also to avoid the singular-drift problem [10].) In order to avoid the singular-drift problem, one can use the deformation technique [8]. In the case of QCD, for instance, one can deform the Dirac operator by adding a deformation parameter in such a way that the eigenvalue distribution does not cover the origin. The result of the original theory can then be recovered by extrapolating the deformation parameter to zero using only the data points that pass the reliability test.
Here we show for the first time that the CLM can be used to investigate QCD at high density and low temperature. We perform simulations on a 4 3 × 8 lattice and calculate the baryon number density as a function of the quark chemical potential. (See also ref. [11] for related work.) We use the gauge cooling for the excursion problem and the deformation technique for the singular-drift problem. To test the reliability of the obtained results, we monitor the probability distribution of the drift term. We find, in particular, that the onset of the baryon number density occurs at smaller chemical potential than in the phase-quenched model, which is a clear sign of the so-called Silver Blaze phenomenon [12].
The rest of this article is organized as follows. In Sec. 2, we first review the application of the CLM to lattice QCD at finite density. We then explain the condition required for the validity of the CLM as well as the techniques we use to meet this condition such as the gauge cooling and the deformation technique. In Sec. 3, we present our results obtained on a 4 3 ×8 lattice. Sec. 4 is devoted to a summary and discussions.

Application of the CLM to lattice QCD
We consider lattice QCD defined on a four-dimensional lattice of the size N t and N s in the temporal and spatial directions, respectively. (The lattice spacing is set to unity throughout this article.) The partition function is given by where x = (x 1 , x 2 , x 3 , x 4 ) labels the sites on the lattice, and U xµ ∈ SU(3) (µ = 1, 2, 3, 4) are the link variables. The gauge action S g (U) is given by where U x,µν = U xµ U x+μ,ν U −1 x+ν,µ U −1 xν withμ being the unit vector in the µ direction. The fermion determinant detM(U, µ) in (1) is defined with the fermion matrix for the staggered fermion with four flavors given by where m is the degenerate quark mass and µ is the quark chemical potential. As usual, we introduce the sign factors η µ = (−1) x 1 +···+x µ−1 , which play the role of the gamma matrices. Note that ǫ x M(U, µ) xy ǫ y = M(U, −µ * ) † holds, where ǫ x = (−1) x 1 +x 2 +x 3 +x 4 plays the role of γ 5 . This identity suggests for µ 0 that detM(U, µ) can be complex and the sign problem occurs.
Here we investigate the lattice QCD (1) at µ 0 using the CLM, in which the link variables U xµ ∈ SU(3) are complexified to U xµ ∈ SL(3, C). The discretized version of the complex Langevin equation is given by where t is the discretized time with the step size ǫ, λ a (a = 1, · · · , 8) are the SU(3) generators normalized by tr(λ a λ b ) = δ ab and η axµ (t) is the real Gaussian noise with the distribution exp{− 1 The CLM enables us to calculate the expectation values of observables that admit the holomorphic extension by taking the ensemble average with thermalized configurations obtained by solving eq. (4). As a typical observable, we consider the baryon number density where N V = N t N 3 s and the symbol · · · represents the expectation value with respect to (1).

Condition for the validity of the CLM
It is known that the results obtained by the CLM are not always correct [3,4], and therefore one has to judge whether they are reliable or not. For that purpose, we need to investigate how the probability distribution of the drift term behaves asymptotically at large magnitude [5]. It is shown under certain assumptions that the CLM gives the correct result if the distribution decays exponentially or faster. Roughly speaking, frequent appearance of a large drift invalidates the equivalence between the CLM and the path integral of the original theory. Note that this problem cannot be solved by making the step size smaller or by using the adaptive step size. In this work, we use the magnitude of the drift term defined by and monitor its probability distribution for thermalized configurations. There are actually two causes for the frequent appearance of a large drift, which makes the CLM fail. One is the excursion problem, which occurs when the complexified link variables U xµ ∈ SL(3, C) go far away from SU(3) [3,4]. The other is the singular-drift problem, which occurs when the fermion matrix has near-zero eigenvalues [13]. In order to avoid these problems and to make the CLM work, one needs to modify the complex Langevin dynamics in some way or another.

Gauge cooling for the excursion problem
For the excursion problem, we use the gauge cooling [6], which amounts to making a complexified gauge transformation after each Langevin update so that the complexified link variables become as close to SU(3) as possible. The gauge transformation is determined by minimizing the unitarity norm which measures the distance of the complexified link variables U xµ ∈ SL(3, C) from SU (3). While the gauge cooling changes the Langevin dynamics nontrivially, the changes do not affect the proof of the equivalence between the CLM and the path integral of the original theory due to the covariance of the drift term and the invariance of the observable under the complexified gauge transformation [14]. Therefore, the proof goes through if the same condition for reliability is satisfied. Indeed this technique played a crucial role in applying the CLM to finite density QCD at high temperature [15].

Deformation for the singular-drift problem
In the high density and low temperature region, we have to deal also with the singular-drift problem.
For that purpose, we use the deformation technique [8], which amounts to deforming the fermion matrix by introducing a parameter so that the singular-drift problem is avoided. In this work, we consider the deformation M(U, µ) xy → M(U, µ) xy + iαη 4 (x)δ xy , where M(U, µ) xy is the fermion matrix defined in (3) and α is the deformation parameter. Note that this deformation corresponds formally to adding an imaginary chemical potential in the continuum theory, which implies that the chiral symmetry and the Lorentz symmetry are kept intact. For a sufficiently large α, the eigenvalue spectrum of the fermion matrix develops a gap in the imaginary direction, and thus the singular-drift problem can be avoided. Using the symmetry under α ↔ −α, we can fit the data points to a polynomial in α 2 and obtain results for the original theory by making α → 0 extrapolations using only the data points that pass the reliability test.

Results
In this section, we show our results for finite density QCD with β = 5.7, m = 0.05 and various values of the quark chemical potential up to µ = 0.7 on a 4 3 × 8 lattice. In applying the CLM, we use a fixed step size ǫ = 10 −4 and the total Langevin time 50 ∼ 150. The gauge cooling and the deformation technique are used to avoid the excursion problem and the singular-drift problem, respectively. We calculate the baryon number density (6) for the deformed theory, and make extrapolations to the zero deformation parameter using only the reliable data points. Let us first explain how we choose the data points to be used for the extrapolation by taking the µ = 0.7 case as an example. In figure 1 (Top-Left), we plot the baryon number density against α 2 .
First we have to choose the data points that pass the reliability test. For that, we look at the probability distribution of the drift term shown in figure 1 (Top-Right), where each line corresponds to each data point in figure 1 (Top-Left). We find that the results for α = 0.2, 0.3 are not reliable because the probability distribution of the drift term falls off with a power law.
The behavior of the drift distribution in figure 1 (Top-Right) is consistent with the eigenvalue distribution of the fermion matrix shown in figure 1 (Bottom-Left) and (Bottom-Right) for α = 0.2 and α = 0.4, respectively. For α = 0.2, we find that the eigenvalue distribution covers the origin, which causes the frequent appearance of a large drift term. For α = 0.4, on the other hand, we observe no eigenvalues close to the origin due to the gap made by the relatively large α, which avoids the appearance of a large drift term.
While the results for α 0.4 are reliable, we cannot use the data points at too large α, which should include higher order corrections with respect to α 2 . Note, in particular, that the baryon number density seems to vanish identically for α ≥ 0.6, which suggests that there is a phase transition at α ∼ 0.6. Hence, for µ = 0.7, we use the data points obtained with α = 0.4, 0.45, 0.5 and make a linear Let us add that the deformation technique is useful not only in avoiding the singular-drift problem, but also in stabilizing the simulations. In fact, at µ = 0.4, the singular-drift problem may not occur even without the deformation [16]. However, the history of observables have large spikes, which makes it difficult to reduce the statistical errors within reasonable amount of statistics. This problem is solved by the deformation even with α as small as 0.1, which clearly suggests another virtue of using the deformation in the CLM.
In figure 3, the circles represent the baryon number density in the original theory obtained by the α → 0 extrapolations shown in figure 2. We find that the baryon number density is zero for µ 0.5 and starts to increase at 0.5 < µ 0.6. Obviously there should be certain systematic errors associated with the α → 0 extrapolations, in particular for µ = 0.7, where we cannot use the data points at small α. Note, however, that the α-dependence of the baryon number density shown in figure 2 changes qualitatively between µ = 0.6 and µ = 0.7. Since the extrapolations take this α-dependence into account, we consider that the rapid increase of the baryon number density at µ ∼ 0.7 is robust.
In figure 3, the squares represent the results obtained by the RHMC simulation of the phasequenched model. In this case, we find that the baryon number density becomes nonzero at µ ∼ 0.4 and grows gradually with increasing µ. Thus we find that the onset of the baryon number density in the CLM occurs at larger µ than in the phase-quenched model. This is consistent with the theoretically expected behavior known as the Silver Blaze phenomenon [12]: at zero temperature, the onset of the baryon number density in full QCD is expected to occur at µ = m N /3, while that in the phase- quenched model occurs at µ = m π /2(< m N /3), where m N is the nucleon mass and m π is the pion mass.
Obviously, the complex phase of the fermion determinant should play a crucial role in realizing the zero baryon number density within m π /2 < µ < m N /3 in full QCD. The delayed onset of the baryon number density seen in our results of the CLM is a clear sign of the Silver Blaze phenomenon, which strongly suggests that the complex phase of the fermion determinant is treated correctly in the CLM. Let us also emphasize that we are able to obtain results with the present setup up to µ/T = 0.7 × 8 = 5.6, where the sign problem is extremely severe and conventional lattice QCD simulations have never been successful.

Summary and discussions
In this work we showed for the first time that it is possible to study QCD at high density and low temperature by the CLM. We have performed lattice QCD simulations with four-flavor staggered fermion on a 4 3 × 8 lattice and calculated the baryon number density as a function of the quark chemical potential within 0.4 ≤ µ ≤ 0.7. We made use of the gauge cooling and the deformation technique for the excursion problem and the singular-drift problem, respectively. The deformation technique was shown to be useful also in stabilizing the simulations. We checked the reliability of the obtained results by the probability distribution of the drift term. The results for the original theory was obtained by extrapolating the reliable results to the zero deformation parameter. The obtained results showed a sharp onset of the baryon number density at µ ∼ 0.6. We also find that the onset of the baryon number density in the CLM occurs at larger µ than in the phase-quenched model, which is a clear sign of the Silver Blaze phenomenon. In order to confirm this phenomenon, we clearly need to increase the lattice size. We have already started simulations with a 8 3 × 16 lattice [17]. In this case, it turned out that correct results can be obtained without deformation and extrapolations even in the region where the quark number density is not very small. In that region, we obtained preliminary results showing that the quark number density increases rapidly twice, which may be interpreted as the transition to the nuclear matter phase and that to the quark matter phase. We hope to report on these results as well as the confirmation of the Silver Blaze phenomenon in the near future.