Complex Langevin calculations in QCD at finite density

We demonstrate that the complex Langevin method (CLM) enables calculations in QCD at finite density in a parameter regime in which conventional methods, such as the density of states method and the Taylor expansion method, are not applicable due to the severe sign problem. Here we use the plaquette gauge action with $\beta = 5.7$ and four-flavor staggered fermions with degenerate quark mass $m a = 0.01$ and nonzero quark chemical potential $\mu$. We confirm that a sufficient condition for correct convergence is satisfied for $\mu /T = 5.2 - 7.2$ on a $8^3 \times 16$ lattice and $\mu /T = 1.6 - 9.6$ on a $16^3 \times 32$ lattice. In particular, the expectation value of the quark number is found to have a plateau with respect to $\mu$ with the height of 24 for both lattices. This plateau can be understood from the Fermi distribution of quarks, and its height coincides with the degrees of freedom of a single quark with zero momentum, which is 3 (color) $\times$ 4 (flavor) $\times$ 2 (spin) $=24$. Our results may be viewed as the first step towards the formation of the Fermi sphere, which plays a crucial role in color superconductivity conjectured from effective theories.


Introduction
QCD at finite temperature and density attracts a lot of interest due to its rich phase diagram predicted by effective theories. Heavy-ion collision experiments are being performed to elucidate the phase structure, whereas the observation of gravitational waves is expected to provide significant information on the equation of state of neutron stars reflecting the phase structure. In parallel, many efforts have been made toward theoretical understanding of QCD at finite temperature and density. The difficulty in theoretical analyses, however, is that nonperturbative calculations based on lattice QCD suffer from the sign problem at finite density. This problem occurs because of the complex fermion determinant, which prevents us from applying standard Monte Carlo methods based on important sampling by identifying the Boltzmann weight as the probability distribution.
Here we focus on the complex Langevin method (CLM) [1,2], which has recently proven a promising method for solving the sign problem. It is a complex extension of the stochastic quantization based on the Langevin equation, where dynamical variables are complexified and physical quantities as well as the drift term are extended holomorphically. An expectation value can be obtained as an average over the fictitious time evolution by the complex Langevin equation after thermalization. See Ref. [3] for a summary of the recent progress concerning this method and other methods for solving the sign problem. In particular, the CLM has been tested extensively in lattice QCD at finite density [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19].
An important issue in the CLM is that physical observables converge to wrong results depending on the parameters, the model, or even on how the method is implemented. It was not until recently that the causes of this incorrect convergence were clarified [20,21,22,23,24,25,26,27]. There are actually two kinds of causes; one is the excursion problem [20] and the other is the singular drift problem [22]. In the case of finite density QCD, the excursion problem occurs when the link variables have long excursions away from the SU(3) group manifold. This problem can be circumvented by adding the gauge cooling procedure after each Langevin update as proposed in Refs. [28,29] and justified in Refs. [23,24]. The singular drift problem occurs, on the other hand, when the fermion matrix has eigenvalues close to zero since the drift term involves the inverse of the fermion matrix.
Both these problems can be detected by just probing the magnitude of the drift term [24], which is calculated anyway for the Langevin evolution. If the histogram of the drift falls off exponentially or faster, one can trust the results, as has been shown from a refined argument for justifying the CLM [24] based on the discrete Langevin-time formulation. The validity of this criterion has been confirmed explicitly in simple one-variable models [24] as well as in exactly solvable models involving infinite degrees of freedom such as chiral Random Matrix Theory [27].
An alternative criterion for correct convergence has been discussed from the viewpoint of the boundary terms [30,31], which appear in the original argument [20,21] for justifying the CLM based on the continuous Langevin-time formulation. Note, however, that the limit of taking the Langevin-time stepsize to zero and the notion of time-evolved observables, which are crucial in the original argument, can be subtle when the drift histogram does not fall off fast enough [24]. These subtleties are taken into account in the refined argument, which led to the above criterion for the validity of the CLM.
In applications to QCD at finite density, it is therefore of primary importance to determine the parameter region in which the CLM gives correct results. We address this issue using staggered fermions corresponding to QCD with four-flavor quarks, which is known to have a first order deconfining phase transition at finite temperature T for zero quark chemical potential µ = 0 [32]. The phase structure of four-flavor QCD on the T − µ plane has been investigated by various methods including the CLM using either staggered fermions [4,6,11,12,13,14,18,33,34,35,36,37,38,39,40] or Wilson fermions [19,41,42,43,44,45]. The strong coupling expansion was also applied as reviewed in Refs. [46,47].
In our calculations, we use Wilson's plaquette action with β = 5.7 on 8 3 ×16 and 16 3 ×32 lattices. The quark mass for the four-flavor staggered fermions is set to ma = 0.01. The criterion for correct convergence of the CLM is found to be satisfied for µ/T = 5.2 − 7.2 on a 8 3 × 16 lattice and µ/T = 1.6 − 9.6 on a 16 3 × 32 lattice, where conventional methods, such as the density of states method and the Taylor expansion method, are not applicable due to the severe sign problem.
In particular, our calculations reveal, for the first time, a plateau behavior of the quark number with respect to the quark chemical potential with the height of 24, which can be understood from the Fermi distribution of quarks. It actually coincides with the number of degrees of freedom for a single quark with zero momentum, which is 3 (the number of colors) × 4 (the number of flavors) × 2 (the number of spin degrees of freedom). This can be regarded as the first step towards the formation of the Fermi surface, which plays a crucial role in color superconductivity. We also investigate other observables such as the Polyakov loop and the chiral condensate. Part of our results have been presented in proceedings articles [12,13,14,18]. This paper is organized as follows. In Section 2 we explain how we apply the CLM to lattice QCD at finite density. In Section 3 we discuss the validity of the CLM in our simulation setup and present our results with emphasis on their physical interpretation. Section 4 is devoted to a summary and discussions.

CLM for QCD at finite density
The partition function of QCD on a four-dimensional lattice is given by where U x,µ ∈ SU(3) is a link variable in the µ(= 1, 2, 3, 4) direction with x = (x 1 , x 2 , x 3 , x 4 ) representing a lattice site. For the gauge action S g [U], we use the plaquette action where β = 6/g 2 with the gauge coupling g andμ represents the unit vector in the µ direction. In this paper, we consider the four-flavor staggered fermions with degenerate quark mass m and quark chemical potential µ, which corresponds to the fermion matrix Periodic boundary conditions are assumed except that we impose anti-periodic boundary conditions on the fermion fields in the temporal direction so that the temporal extent of the lattice represents the inverse temperature T −1 . We apply the CLM to this system. In this method, the link variables U x,µ ∈ SU(3) are complexified to U x,µ ∈ SL(3,C), and we consider a fictitious time evolution governed by which is the discrete version of the complex Langevin equation with the stepsize ε.
Here η x,µ (t) is the noise term, which is a traceless 3 × 3 Hermitian matrix obeying the Gaussian where the symbol · η represents the expectation value with respect to the Gaussian noise η.
The v x,µ in eq. (2.4) is the drift term, which is defined by holomorphically extending using the notation of continuum Langevin time t for simplicity. Here t 0 represents the time needed for thermalization and τ represents the total Langevin time for taking the average, which should be large enough to achieve good statistics. The CLM is justified if the expectation valueŌ obtained by the CLM agrees with the expectation value The criterion for justification [24] is based on the magnitude of the drift term (2 If the histogram of this quantity falls off exponentially or faster, we can trust the results. This can be violated either by the excursion problem or by the singular drift problem as we mentioned in the Introduction. In the former case, it is the drift coming from the gauge action that shows slow fall-off in the histogram, while in the latter case, it is the drift coming from the fermion determinant.
The excursion problem can also be probed by the unitarity norm which represents the distance of a configuration from the SU(3) manifold with L t and L s being the lattice size in the temporal and spatial directions, respectively. The unitarity norm (2.10) is positive semi-definite and it becomes zero if and only if all the link variables are unitary. Rapid growth of the unitarity norm typically signals the occurrence of the excursion problem. We note, however, that recent work [48] on 2D U(1) theory with a θ term suggests that the CLM works even if the unitarity norm becomes large as far as the drift histogram falls off fast. Therefore, one cannot tell the validity of the CLM by looking at the unitarity norm alone.
In order to avoid the excursion problem, we use the gauge cooling [28,29], which amounts to making a complexified gauge transformation in such a way that the unitarity norm is minimized. Adding this procedure after each step of the Langevin-time evolution does not spoil the argument for justification as is shown in Refs. [23,24]. The gauge cooling keeps all the link variables as close to unitary matrices as possible during the Langevin-time evolution.
As for physical observables, we calculate the Polyakov loop, the quark number and the chiral condensate. The Polyakov loop is given by with x = (x 1 , x 2 , x 3 ) being a spatial coordinate on the lattice. The quark number N q and the chiral condensate Σ are defined by where the latter trace Tr is taken not only for the color index but also for the spacetime index. These two quantities (2.13) and (2.14) are calculated by the so-called noisy estimator using 20 noise vectors.

Results
We

Validity of the CLM
Let us first discuss the validity of the CLM in our simulation setup. In Fig. 1 we plot the distribution p(v) for the magnitude of the drift term (2.9) coming from the gauge action (Left) and from the fermion determinant (Right), respectively, obtained by simulations on a 8 3 × 16 lattice with β = 5.7 andm = 0.01. For the sake of visibility, we separate the data points into two regionsμ ≤ 0.3 andμ ≥ 0.325 and show them in the upper and lower panels, respectively. We find that the drift term coming from the gauge action shows slow fall-off forμ = 0.2, which implies that the excursion problem occurs only in this case.    Figure 2 shows similar plots for a 16 3 × 32 lattice with the same β = 5.7 andm = 0.01. Here we find that the excursion problem does not occur for any values ofμ, whereas the singular drift problem occurs forμ = 0.325. Thus we find that the CLM is expected to give correct results for 0.05 ≤μ ≤ 0.3 on the 16 3 × 32 lattice.
In Fig. 3, we plot the Langevin-time histories of the unitarity norm (2.10). For the 8 3 × 16 lattice (Left), we find that the history forμ = 0.2 looks quite violent, which is consistent with what we observe in Fig. 1 (Top-Left). For the 16 3 × 32 lattice (Right), on the other hand, the unitarity norm seems to be well under control; Note that the scale of the vertical axis here is an order of magnitude smaller than that in the Left panel. This is also consistent with what we observe in Fig. 2 (Left). Let us emphasize, however, that from the histories of the unitarity norm alone, we cannot judge the validity of the CLM unambiguously. Note also that the unitarity norm has a long autocorrelation time as one can see from Fig. 3. Fortunately, we find that the physical observables we investigate are not correlated with the unitarity norm, and they have a much shorter autocorrelation time. This is important because it allows us to calculate their expectation values reliably within a reasonable length of the total Langevin time.

Physical observables
In what follows, we present only the data points in the parameter region in which the criterion for justification is satisfied. In Fig. 4 (Top) we plot the real part of the Polyakov loop (2.12) againstμ for β = 5.7 andm = 0.01 on the 8 3 × 16 and the 16 3 × 32 lattices. The results for the 8 3 × 16 lattice are slightly nonzero and the results for the 16 3 × 32 lattice are consistent with zero. Let us recall here that the Polyakov loop is an order parameter for the deconfining transition. The interpretation of our results requires some care, though. Note that the spatial size of our lattice is aL s = 0.36 fm and 0.68 fm for L s = 8 and 16, respectively, which are smaller than the typical length scale of QCD, namely Λ −1 LQCD ∼ 1 fm. Thus the situation we are simulating should be regarded as QCD in a small box, where the notion of quark confinement does not make sense. In fact, the temperature is T ∼ 290 MeV and 145 MeV for L t = 16 and 32, respectively, which are higher than or close to the critical temperature T c ∼ 170 MeV for the deconfining transition in QCD with four-flavor staggered quarks [51] with m/T = 0.2, where large physical volume is implicitly assumed. The Polyakov loop being either small or zero in our setup simply confirms that we are probing the "low temperature" behavior of such a finite size system due to the chosen aspect ratio of our lattice.
In Fig. 4 (Middle) we plot the quark number (2.13) againstμ for β = 5.7 andm = 0.01  On both lattices, we observe a plateau at the height of N q = 24. In order to understand this behavior, let us recall that the physical extent of our lattice is too small to create a baryon in it. The effective gauge coupling is small due to the asymptotic freedom, which makes the Fermi distribution of quarks qualitatively valid. At sufficiently largeμ, the path integral is therefore dominated by a state obtained from the vacuum by creating quarks with momentum p satisfying p 2 + m 2 eff ≤μ, where m eff is the effective mass including quantum corrections 1 . It should be noted here that the momentum is discretized in a finite box as p = (2π/L s ) n with n being a 3D integer vector. In particular, for m eff ≤μ ≤ µ 1 , where µ 1 = (2π/L s ) 2 + m 2 eff , only the quarks with p = 0 are created. The height of the plateau is therefore given by the internal degrees of freedom N f × N c × N spin = 4 × 3 × 2 = 24, where N f , N c and N spin are the number of flavors, the number of colors and the number of spin degrees of freedom, respectively 2 . In Fig. 4 (Middle), we also observe that our data start to leave the plateau for largerμ, which can be understood as the effects of quarks with the first non-zero momenta being created. The value ofμ at which this growth of N q occurs is smaller than µ 1 defined above, which can be understood as a result of finite temperature effects. Note that µ 1 becomes smaller as the lattice becomes larger, which is clearly reflected in our results.
The appearance of such plateaus in the quark number for QCD in a finite box was discussed in the case of free theory using the naive lattice action for fermions [53]. See also Ref. [54] for such behaviors based on one-loop perturbative calculations in the continuum QCD on a small S 3 . In a separate paper [55], we will report on our results of the CLM for larger β, which are compared with perturbative results obtained with staggered fermions. This plateau behavior should not be confused with the quark number saturation that occurs at much largerμ. In that case, the path integral is dominated by a state with all the sites being occupied by fermions. Taking the internal degrees of freedom into account, the height of the plateau becomes N f × N c × N spin × L 3 s = 24 × L 3 s , which is much higher than 24. In Fig. 4 (Bottom) we show our results for the chiral condensate (2.14). The plateau behaviors appear here as well because of the "low temperature", where changingμ a little cannot create quarks at higher energy levels. The plateau corresponding to the state with the zero-momentum quarks appears with the height only slightly lower than that corresponding to the state that dominates atμ = 0. This suggests that the chiral symmetry for m = 0, which is considered to be spontaneously broken atμ = 0, does not get restored by increasingμ in the present parameter regime.
The plateau behaviors observed above in the quark number and the chiral condensate are analogous to those in two-color QCD 3 using two-flavor Wilson fermions on a 3 3 × 64 lattice at β = 24 with finite µ [56]. In that case, however, the height of the plateau in the quark number does not agree with the free fermion results, which is in contrast to our results for the SU(3) gauge group on 8 3 × 16 and 16 3 × 32 lattices.

Summary and Discussions
In this paper we have applied the CLM to QCD at finite density with the plaquette gauge action and the four-flavor staggered fermions on 8 3 × 16 and 16 3 × 32 lattices. While the spatial size of our lattice is still as small as aL s = 0.36 fm and 0.68 fm for L s = 8 and 16, we find that the criterion for correct convergence is satisfied for µ = 1.5 − 2.1 GeV on the 8 3 × 16 lattice and for µ = 0.23 − 1.4 GeV on the 16 3 × 32 lattice with T ∼ 290 MeV and 145 MeV, respectively. These parameter regimes cannot be reached by conventional methods, such as the density of states method and the Taylor expansion method. Thus our results clearly demonstrate a big advantage of the CLM in overcoming the sign problem in finite density QCD.
Let us also mention that the previous work [11] shows that the CLM works on a 4 3 × 8 lattice with the same β but only with the aid of the deformation technique [57], which is actually not needed for the lattice size in the present work. Thus we find that the situation becomes better for a larger lattice, which is also seen by comparing our results for 8 3 × 16 and 16 3 × 32 lattices in section 3.1.
One of our main physical results is that the quark number exhibits a plateau behavior as a function of the quark chemical potential with the height of 24 at sufficiently large µ. This has been interpreted as the creation of quarks with zero momentum, which has the internal degrees of freedom N f × N c × N spin = 4 × 3 × 2 = 24. We may regard it as the first step towards the formation of the Fermi surface, which plays a crucial role in color superconductivity.
It is of particular importance to perform similar calculations with larger lattices. That will enable us to observe the growth of the Fermi sphere with moderate values of µ thanks to better momentum resolution. Note that color superconductivity is expected to occur due to Cooper pairing of quarks near the Fermi surface, which is possible even at weak coupling or in a small physical volume [58]. For this reason, we are currently exploring the larger β regime, where we can compare our results against perturbative calculations [55].
In particular, we are trying to observe a departure from perturbative behaviors as β gets smaller than some critical value, which would signal the onset of color superconductivity.