Complex Langevin calculations in QCD at finite density
Yuta Ito, Hideo Matsufuru, Yusuke Namekawa, Jun Nishimura, Shinji Shimasaki, Asato Tsuchiya, Shoichiro Tsutsui
aa r X i v : . [ h e p - l a t ] J u l KEK-TH-2230, RIKEN-QHP-479
Complex Langevin calculations in QCD at finite density
Yuta I to ab , Hideo M atsufuru c , Yusuke N amekawa a , Jun N ishimura ad ,Shinji S himasaki e , Asato T suchiya f , Shoichiro T sutsui g a Theory Center, Institute of Particle and Nuclear Studies,High Energy Accelerator Research Organization (KEK),1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan b National Institute of Technology, Tokuyama College,Gakuendai, Shunan, Yamaguchi 745-8585, Japan c Computing Research Center, High Energy Accelerator Research Organization (KEK),1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan d Department of Particle and Nuclear Physics,School of High Energy Accelerator Science,Graduate University for Advanced Studies (SOKENDAI),1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan e Research and Education Center for Natural Sciences, Keio University,Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan f Department of Physics, Shizuoka University,836 Ohya, Suruga-ku, Shizuoka 422-8529, Japan g Theoretical Research Division, Nishina Center, RIKEN,Wako, Saitama 351-0198, Japan E-mail address : [email protected] E-mail address : [email protected] E-mail address : [email protected]; Present address is Yukawa Institute for Theoretical Physics,Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502 Japan E-mail address : [email protected] E-mail address : [email protected] E-mail address : [email protected] E-mail address : [email protected] bstract
We demonstrate that the complex Langevin method (CLM) enables calculations inQCD at finite density in a parameter regime in which conventional methods, such asthe density of states method and the Taylor expansion method, are not applicabledue to the severe sign problem. Here we use the plaquette gauge action with β =5 . ma = 0 .
01 andnonzero quark chemical potential µ . We confirm that a sufficient condition for correctconvergence is satisfied for µ/T = 5 . − . ×
16 lattice and µ/T = 1 . − . ×
32 lattice. In particular, the expectation value of the quark number is found tohave a plateau with respect to µ with the height of 24 for both lattices. This plateaucan be understood from the Fermi distribution of quarks, and its height coincideswith the degrees of freedom of a single quark with zero momentum, which is 3 (color) × × Introduction
QCD at finite temperature and density attracts a lot of interest due to its rich phase diagrampredicted by effective theories. Heavy-ion collision experiments are being performed toelucidate the phase structure, whereas the observation of gravitational waves is expectedto provide significant information on the equation of state of neutron stars reflecting thephase structure. In parallel, many efforts have been made toward theoretical understandingof 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 atfinite density. This problem occurs because of the complex fermion determinant, whichprevents us from applying standard Monte Carlo methods based on important sampling byidentifying the Boltzmann weight as the probability distribution.Here we focus on the complex Langevin method (CLM) [1, 2], which has recentlyproven a promising method for solving the sign problem. It is a complex extension ofthe stochastic quantization based on the Langevin equation, where dynamical variables arecomplexified and physical quantities as well as the drift term are extended holomorphi-cally. An expectation value can be obtained as an average over the fictitious time evolutionby the complex Langevin equation after thermalization. See Ref. [3] for a summary ofthe recent progress concerning this method and other methods for solving the sign prob-lem. 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 resultsdepending on the parameters, the model, or even on how the method is implemented. Itwas 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, theexcursion 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 procedureafter each Langevin update as proposed in Refs. [28, 29] and justified in Refs. [23, 24]. Thesingular drift problem occurs, on the other hand, when the fermion matrix has eigenvaluesclose 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 driftfalls off exponentially or faster, one can trust the results, as has been shown from a refinedargument for justifying the CLM [24] based on the discrete Langevin-time formulation. Thevalidity of this criterion has been confirmed explicitly in simple one-variable models [24]1s well as in exactly solvable models involving infinite degrees of freedom such as chiralRandom Matrix Theory [27].An alternative criterion for correct convergence has been discussed from the viewpointof the boundary terms [30, 31], which appear in the original argument [20, 21] for justifyingthe CLM based on the continuous Langevin-time formulation. Note, however, that the limitof 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 notfall 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 deter-mine the parameter region in which the CLM gives correct results. We address this issueusing staggered fermions corresponding to QCD with four-flavor quarks, which is known tohave a first order deconfining phase transition at finite temperature T for zero quark chem-ical potential µ = 0 [32]. The phase structure of four-flavor QCD on the T − µ plane hasbeen 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 . ×
16 and 16 × ma = 0 .
01. Thecriterion for correct convergence of the CLM is found to be satisfied for µ/T = 5 . − . ×
16 lattice and µ/T = 1 . − . ×
32 lattice, where conventional methods,such as the density of states method and the Taylor expansion method, are not applicabledue to the severe sign problem.In particular, our calculations reveal, for the first time, a plateau behavior of the quarknumber with respect to the quark chemical potential with the height of 24, which can beunderstood from the Fermi distribution of quarks. It actually coincides with the numberof degrees of freedom for a single quark with zero momentum, which is 3 (the number ofcolors) × × CLM for QCD at finite density
The partition function of QCD on a four-dimensional lattice is given by Z = Z Y x,µ dU x,µ det M [ U ] e − S g [ U ] , (2.1)where U x,µ ∈ SU(3) is a link variable in the µ (= 1 , , ,
4) direction with x = ( x , x , x , x )representing a lattice site. For the gauge action S g [ U ], we use the plaquette action S g [ U ] = − β X x X µ = ν tr (cid:0) U x,µ U x +ˆ µ,ν U − x +ˆ ν,µ U − x,ν (cid:1) , (2.2)where β = 6 /g with the gauge coupling g and ˆ µ represents the unit vector in the µ direction. In this paper, we consider the four-flavor staggered fermions with degeneratequark mass m and quark chemical potential µ , which corresponds to the fermion matrix M xy [ U ] = ˜ mδ xy + X µ =1 η µ ( x )2 (cid:0) e ˜ µδ µ U x,µ δ x +ˆ µ,y − e − ˜ µδ µ U − x − ˆ µ,µ δ x − ˆ µ,y (cid:1) , (2.3)with η µ ( x ) ≡ ( − x + ··· + x µ − being a site-dependent sign factor. In (2.3) we have defineddimensionless parameters ˜ m = ma and ˜ µ = µa , where a represents the lattice spacing. Thedeterminant det M [ U ] in (2.1) becomes complex for ˜ µ = 0, which causes the sign problem.Periodic boundary conditions are assumed except that we impose anti-periodic boundaryconditions on the fermion fields in the temporal direction so that the temporal extent ofthe lattice represents the inverse temperature T − .We apply the CLM to this system. In this method, the link variables U x,µ ∈ SU(3) arecomplexified to U x,µ ∈ SL(3,C), and we consider a fictitious time evolution governed by U x,µ ( t + ε ) = exp (cid:8) i (cid:0) − εv x,µ [ U ( t )] + √ εη x,µ ( t ) (cid:1)(cid:9) U x,µ ( t ) , (2.4)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 × − tr η x,µ ( t )] so that h η ijx,µ ( s ) η kly,ν ( t ) i η = 2 δ xy δ µν δ st (cid:18) δ il δ jk − δ ij δ kl (cid:19) , (2.5)where the symbol h · i η 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 v x,µ [ U ] = X a λ a ∂∂α S [ e iαλ a U x,µ ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α =0 (2.6)3efined for unitary configurations U , where S [ U ] = S g [ U ] − ln det M [ U ] and λ a ( a = 1 , . . . , λ a λ b ) = δ ab . Note that the drift term is notHermitian for ˜ µ = 0 even for unitary link variables, which makes the time-evolved linkvariables inevitably non-unitary.The expectation value of an observable O ( U ) is calculated as¯ O = 1 τ Z t + τt dt O ( U ( t )) (2.7)using the notation of continuum Langevin time t for simplicity. Here t represents thetime needed for thermalization and τ represents the total Langevin time for taking theaverage, which should be large enough to achieve good statistics. The CLM is justified ifthe expectation value ¯ O obtained by the CLM agrees with the expectation value hO ( U ) i = 1 Z Z dU O ( U ) det M [ U ] e − S g [ U ] (2.8)defined in the path integral formalism.The criterion for justification [24] is based on the magnitude of the drift term (2.6) v = max x,µ r
13 tr (cid:16) v † x,µ v x,µ (cid:17) . (2.9)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 aswe mentioned in the Introduction. In the former case, it is the drift coming from the gaugeaction that shows slow fall-off in the histogram, while in the latter case, it is the driftcoming from the fermion determinant.The excursion problem can also be probed by the unitarity norm N = 112 L t L X x,µ tr ( U † x,µ U x,µ − ) , (2.10)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 unitaritynorm (2.10) is positive semi-definite and it becomes zero if and only if all the link variablesare unitary. Rapid growth of the unitarity norm typically signals the occurrence of theexcursion 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 thedrift histogram falls off fast. Therefore, one cannot tell the validity of the CLM by lookingat the unitarity norm alone. 4n order to avoid the excursion problem, we use the gauge cooling [28, 29], which amountsto making a complexified gauge transformation δ g U x,µ = g x U x,µ g − x +ˆ µ , g ∈ SL(3,C) (2.11)in such a way that the unitarity norm is minimized. Adding this procedure after each stepof the Langevin-time evolution does not spoil the argument for justification as is shown inRefs. [23, 24]. The gauge cooling keeps all the link variables as close to unitary matrices aspossible during the Langevin-time evolution.As for physical observables, we calculate the Polyakov loop, the quark number and thechiral condensate. The Polyakov loop is given by P = 13 L X ~x tr L t Y x =1 U ( ~x,x ) , ! , (2.12)with ~x = ( x , x , x ) being a spatial coordinate on the lattice. The quark number N q andthe chiral condensate Σ are defined by N q = 1 L t ∂∂ ˜ µ log Z = 1 L t *X x η ( x )2 tr (cid:16) e ˜ µ M − x +ˆ4 ,x U x, + e − ˜ µ M − x − ˆ4 ,x U − x − ˆ4 , (cid:17)+ , (2.13)Σ = 1 L L t ∂∂ ˜ m log Z = (cid:10) Tr M − (cid:11) , (2.14)where the latter trace Tr is taken not only for the color index but also for the spacetimeindex. These two quantities (2.13) and (2.14) are calculated by the so-called noisy estimatorusing 20 noise vectors. We have performed complex Langevin simulations at β = 5 . m = 0 .
01 on 8 ×
16 and 16 ×
32 lattices. We determine the lattice spacing as a − =4 . ×
48 lattice with ˜ µ = 0. The quark chemical potential is varied withinthe range ˜ µ = 0 . − . ×
16 lattice and ˜ µ = 0 . − .
325 on the 16 ×
32 lattice.The Langevin-time stepsize is chosen initially as ε = 1 . × − and reduced adaptively [50]when the magnitude (2.9) of the drift becomes larger than the threshold 0.01. The totalLangevin time after thermalization is τ = 70 −
140 for the 8 ×
16 lattice and τ = 10 − ×
32 lattice. The error bars of the physical observables are evaluated by thejackknife method with the bin sizes of 2 Langevin time on 8 ×
16 and 0 . ×
32. 5 -10 -8 -6 -4 -2 p ( v g ) v g β =5.7, m~=0.01 on 8 x 16 µ ~=0.100 µ ~=0.150 µ ~=0.200 µ ~=0.300 10 -10 -8 -6 -4 -2 p ( v f ) v f β =5.7, m~=0.01 on 8 x 16 µ ~=0.100 µ ~=0.150 µ ~=0.200 µ ~=0.30010 -10 -8 -6 -4 -2 p ( v g ) v g β =5.7, m~=0.01 on 8 x 16 µ ~=0.325 µ ~=0.400 µ ~=0.450 µ ~=0.475 µ ~=0.500 10 -10 -8 -6 -4 -2 p ( v f ) v f β =5.7, m~=0.01 on 8 x 16 µ ~=0.325 µ ~=0.400 µ ~=0.450 µ ~=0.475 µ ~=0.500 Figure 1: The probability distributions p ( v ) of the drift term obtained from simulations ona 8 ×
16 lattice with β = 5 . m = 0 .
01 are plotted for the drift term v g coming fromthe gauge action (Left) and v f from the fermion determinant (Right). The upper and lowerpanels show the data points for ˜ µ ≤ . µ ≥ . Let us first discuss the validity of the CLM in our simulation setup. In Fig. 1 we plot thedistribution 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 ona 8 ×
16 lattice with β = 5 . m = 0 .
01. For the sake of visibility, we separate thedata points into two regions ˜ µ ≤ . µ ≥ .
325 and show them in the upper and lowerpanels, respectively. We find that the drift term coming from the gauge action shows slowfall-off for ˜ µ = 0 .
2, which implies that the excursion problem occurs only in this case. Thedrift term coming from the fermion determinant, on the other hand, shows slow fall-off for0 . ≤ ˜ µ ≤ . . ≤ ˜ µ indicating the occurrence of the singular drift problem inthese cases. Thus we conclude that the CLM gives correct results in the regions ˜ µ ≃ . . ≤ ˜ µ ≤ .
45 on the 8 ×
16 lattice. 6 -10 -8 -6 -4 -2 p ( v g ) v g β =5.7, m~=0.01 on 16 x 32 µ ~=0.050 µ ~=0.100 µ ~=0.200 µ ~=0.300 µ ~=0.325 10 -10 -8 -6 -4 -2 p ( v f ) v f β =5.7, m~=0.01 on 16 x 32 µ ~=0.050 µ ~=0.100 µ ~=0.200 µ ~=0.300 µ ~=0.325 Figure 2: The probability distributions p ( v ) of the drift term obtained from simulations ona 16 ×
32 lattice with β = 5 . m = 0 .
01 are plotted for the drift term v g coming fromthe gauge action (Left) and v f from the fermion determinant (Right). N Langevin time β =5.7, m~=0.01 on 8 x16 µ ~=0.100 µ ~=0.200 µ ~=0.300 µ ~=0.500 0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 0.0014 0.0016 0.0018 0 10 20 30 40 50 60 70 80 N Langevin time β =5.7, m~=0.01 on 16 x32 µ ~=0.050 µ ~=0.100 µ ~=0.200 µ ~=0.300 µ ~=0.325 Figure 3: The Langevin-time histories of the unitarity norm (starting from the initialconfiguration) are plotted for β = 5 . m = 0 .
01 on the 8 ×
16 lattice (Left) and the16 ×
32 lattice (Right). 7igure 2 shows similar plots for a 16 ×
32 lattice with the same β = 5 . m = 0 . µ , whereas thesingular drift problem occurs for ˜ µ = 0 . . ≤ ˜ µ ≤ . ×
32 lattice.In Fig. 3, we plot the Langevin-time histories of the unitarity norm (2.10). For the8 ×
16 lattice (Left), we find that the history for ˜ µ = 0 . ×
32 lattice (Right), onthe other hand, the unitarity norm seems to be well under control; Note that the scale ofthe vertical axis here is an order of magnitude smaller than that in the Left panel. Thisis also consistent with what we observe in Fig. 2 (Left). Let us emphasize, however, thatfrom the histories of the unitarity norm alone, we cannot judge the validity of the CLMunambiguously. Note also that the unitarity norm has a long autocorrelation time as onecan see from Fig. 3. Fortunately, we find that the physical observables we investigate arenot 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 withina reasonable length of the total Langevin time.
In what follows, we present only the data points in the parameter region in which thecriterion for justification is satisfied. In Fig. 4 (Top) we plot the real part of the Polyakovloop (2.12) against ˜ µ for β = 5 . m = 0 .
01 on the 8 ×
16 and the 16 ×
32 lattices.The results for the 8 ×
16 lattice are slightly nonzero and the results for the 16 × 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Λ − ∼ T ∼
290 MeV and 145 MeV for L t = 16 and 32, respectively, which are higher than orclose to the critical temperature T c ∼
170 MeV for the deconfining transition in QCD withfour-flavor staggered quarks [51] with m/T = 0 .
2, where large physical volume is implicitlyassumed. The Polyakov loop being either small or zero in our setup simply confirms that weare probing the “low temperature” behavior of such a finite size system due to the chosenaspect ratio of our lattice.In Fig. 4 (Middle) we plot the quark number (2.13) against ˜ µ for β = 5 . m = 0 . 〈 R e ( P ) 〉 µ ~ β =5.7, m~ =0.018 x1616 x32 0 8 16 24 32 40 0 0.1 0.2 0.3 0.4 0.5 0.6 N q µ ~ β =5.7, m~ =0.018 x1616 x32 0 0.02 0.04 0.06 0.08 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 Σ µ ~ β =5.7, m~ =0.018 x1616 x32 Figure 4: The real part of the Polyakov loop P (Top), the quark number N q (Middle) andthe chiral condensate Σ (Bottom) are plotted against the quark chemical potential ˜ µ for β = 5 . m = 0 .
01 on the 8 ×
16 and 16 ×
32 lattices.9n the 8 ×
16 and 16 ×
32 lattices. On both lattices, we observe a plateau at the heightof N q = 24. In order to understand this behavior, let us recall that the physical extentof our lattice is too small to create a baryon in it. The effective gauge coupling is smalldue to the asymptotic freedom, which makes the Fermi distribution of quarks qualitativelyvalid. At sufficiently large ˜ µ , the path integral is therefore dominated by a state obtainedfrom the vacuum by creating quarks with momentum ~p satisfying p ~p + m ≤ ˜ µ , where m eff is the effective mass including quantum corrections . It should be noted here thatthe momentum is discretized in a finite box as ~p = (2 π/L s ) ~n with ~n being a 3D integervector. In particular, for m eff ≤ ˜ µ ≤ µ , where µ = p (2 π/L s ) + m , only the quarkswith ~p = 0 are created. The height of the plateau is therefore given by the internal degreesof freedom N f × N c × N spin = 4 × × N f , N c and N spin are the number offlavors, the number of colors and the number of spin degrees of freedom, respectively . InFig. 4 (Middle), we also observe that our data start to leave the plateau for larger ˜ µ , whichcan 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 µ defined above, whichcan be understood as a result of finite temperature effects. Note that µ becomes smalleras 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 wasdiscussed in the case of free theory using the naive lattice action for fermions [53]. See alsoRef. [54] for such behaviors based on one-loop perturbative calculations in the continuumQCD on a small S . In a separate paper [55], we will report on our results of the CLM forlarger β , which are compared with perturbative results obtained with staggered fermions.This plateau behavior should not be confused with the quark number saturation that occursat much larger ˜ µ . In that case, the path integral is dominated by a state with all the sitesbeing occupied by fermions. Taking the internal degrees of freedom into account, the heightof the plateau becomes N f × N c × N spin × L = 24 × L , which is much higher than 24.In Fig. 4 (Bottom) we show our results for the chiral condensate (2.14). The plateaubehaviors appear here as well because of the “low temperature”, where changing ˜ µ a littlecannot create quarks at higher energy levels. The plateau corresponding to the state withthe zero-momentum quarks appears with the height only slightly lower than that corre-sponding 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 byincreasing ˜ µ in the present parameter regime. Here we neglect finite lattice spacing effects assuming that we are close to the continuum limit. In the case of two-flavor Wilson fermions, the height of the plateau is found to be 12, which agreeswith N f × N c × N spin = 2 × × using two-flavor Wilson fermions on a 3 × β = 24 with finite µ [56]. In that case, however, the height of the plateau inthe quark number does not agree with the free fermion results, which is in contrast to ourresults for the SU(3) gauge group on 8 ×
16 and 16 ×
32 lattices.
In this paper we have applied the CLM to QCD at finite density with the plaquette gaugeaction and the four-flavor staggered fermions on 8 ×
16 and 16 ×
32 lattices. While thespatial size of our lattice is still as small as aL s = 0 .
36 fm and 0 .
68 fm for L s = 8 and16, we find that the criterion for correct convergence is satisfied for µ = 1 . − . ×
16 lattice and for µ = 0 . − . ×
32 lattice with T ∼
290 MeVand 145 MeV, respectively. These parameter regimes cannot be reached by conventionalmethods, such as the density of states method and the Taylor expansion method. Thus ourresults clearly demonstrate a big advantage of the CLM in overcoming the sign problem infinite density QCD.Let us also mention that the previous work [11] shows that the CLM works on a 4 × β but only with the aid of the deformation technique [57], which isactually not needed for the lattice size in the present work. Thus we find that the situationbecomes better for a larger lattice, which is also seen by comparing our results for 8 × ×
32 lattices in section 3.1.One of our main physical results is that the quark number exhibits a plateau behavioras 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 theinternal degrees of freedom N f × N c × N spin = 4 × × µ thanksto better momentum resolution. Note that color superconductivity is expected to occurdue to Cooper pairing of quarks near the Fermi surface, which is possible even at weakcoupling or in a small physical volume [58]. For this reason, we are currently exploring thelarger β regime, where we can compare our results against perturbative calculations [55]. This model has no sign problem even at finite density and hence the standard hybrid Monte Carloalgorithm is applicable.
11n particular, we are trying to observe a departure from perturbative behaviors as β getssmaller than some critical value, which would signal the onset of color superconductivity. Acknowledgements
We thank E. Itou, A. Ohnishi and M. Scherzer for important discussions, which enabledus to establish the interpretation of the plateau behavior observed in our results for thequark number. The authors are also grateful to M.P. Lombardo for bringing our atten-tion to Ref. [53] and to Y. Asano, T. Kaneko and T. Yokota for collaborations in theongoing projects related to the present work. This research was supported by MEXT as“Priority Issue on Post-K computer” (Elucidation of the Fundamental Laws and Evolu-tion of the Universe) and as “Program for Promoting Researches on the SupercomputerFugaku” (Simulation for basic science: from fundamental laws of particles to creation ofnuclei). It is also supported by Joint Institute for Computational Fundamental Science(JICFuS). Computations were carried out using computational resources of the K com-puter provided by the RIKEN Center for Computational Science through the HPCI Sys-tem Research project (Project ID:hp180178, hp190159) and the Oakbridge-CX providedby the Information Technology Center at the University of Tokyo through the HPCI Sys-tem Research project (Project ID:hp200079). Y. N. and J. N. were supported in part byJSPS KAKENHI Grant Numbers JP16H03988, JP18K03638. S. S. was supported by theMEXT-Supported Program for the Strategic Research Foundation at Private Universities“Topological Science” (Grant No. S1511006). S. T. was supported by the RIKEN SpecialPostdoctoral Researchers Program.
References [1] G. Parisi,
On complex probabilities , Phys. Lett. (1983) 393–395.[2] J. R. Klauder,
Coherent state Langevin equations for canonical quantum systems withapplications to the quantized Hall effect , Phys. Rev.
A29 (1984) 2036–2047.[3] C. E. Berger, L. Rammelm¨uller, A. C. Loheac, F. Ehmann, J. Braun, and J. E. Drut,
Complex Langevin and other approaches to the sign problem in quantum many-bodyphysics , arXiv:1907.10183 .[4] D. Sexty, Simulating full QCD at nonzero density using the complex Langevinequation , Phys. Lett.
B729 (2014) 108–111, [ arXiv:1307.7748 ].125] G. Aarts, E. Seiler, D. Sexty, and I.-O. Stamatescu,
Simulating QCD at nonzerobaryon density to all orders in the hopping parameter expansion , Phys. Rev.
D90 (2014), no. 11 114505, [ arXiv:1408.3770 ].[6] Z. Fodor, S. Katz, D. Sexty, and C. Torok,
Complex Langevin dynamics fordynamical QCD at nonzero chemical potential: A comparison with multiparameterreweighting , Phys. Rev. D (2015), no. 9 094516, [ arXiv:1508.05260 ].[7] D. K. Sinclair and J. Kogut, Exploring complex-Langevin methods for finite-densityQCD , PoS
LATTICE2015 (2016) 153, [ arXiv:1510.06367 ].[8] D. Sinclair and J. Kogut,
Complex Langevin for lattice QCD at T = 0 and µ ≥ PoS
LATTICE2016 (2016) 026, [ arXiv:1611.02312 ].[9] D. Sinclair and J. Kogut,
Complex Langevin simulations of QCD at finite density –Progress report , EPJ Web Conf. (2018) 07031, [ arXiv:1710.08465 ].[10] D. Sinclair and J. Kogut,
Complex Langevin for lattice QCD , PoS
LATTICE2018 (2018) 143, [ arXiv:1810.11880 ].[11] K. Nagata, J. Nishimura, and S. Shimasaki,
Complex Langevin calculations in finitedensity QCD at large µ/T with the deformation technique , Phys. Rev.
D98 (2018),no. 11 114513, [ arXiv:1805.03964 ].[12] Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki, A. Tsuchiya, and S. Tsutsui,
Exploring the phase diagram of finite density QCD at low temperature by the complexLangevin method , PoS
LATTICE2018 (2018) 146, [ arXiv:1811.12688 ].[13] S. Tsutsui, Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki, and A. Tsuchiya,
Canthe complex Langevin method see the deconfinement phase transition in QCD at finitedensity? , PoS
LATTICE2018 (2018) 144, [ arXiv:1811.07647 ].[14] S. Tsutsui, Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki, and A. Tsuchiya,
Exploring the finite density QCD based on the complex Langevin method , JPS Conf.Proc. (2019) 024012.[15] J. B. Kogut and D. K. Sinclair, Applying complex Langevin simulations to latticeQCD at finite density , Phys. Rev.
D100 (2019), no. 5 054512, [ arXiv:1903.02622 ].[16] D. Sexty,
Calculating the equation of state of dense quark-gluon plasma using thecomplex Langevin equation , Phys. Rev. D (2019), no. 7 074503,[ arXiv:1907.08712 ]. 1317] D. K. Sinclair and J. B. Kogut,
Applying complex Langevin to lattice QCD at finite µ , in , vol. LATTICE2019, p. 245, 2019. arXiv:1910.11412 .[18] S. Tsutsui, Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki, and A. Tsuchiya, Exploring the QCD phase diagram at finite density by the complex Langevin methodon a × lattice , in , 2019. arXiv:1912.00361 .[19] M. Scherzer, D. Sexty, and I.-O. Stamatescu, Deconfinement transition line with thecomplex Langevin equation up to µ/T ∼ arXiv:2004.05372 .[20] G. Aarts, E. Seiler, and I.-O. Stamatescu, The complex Langevin method: When canit be trusted? , Phys. Rev.
D81 (2010) 054508, [ arXiv:0912.3360 ].[21] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu,
Complex Langevin: Etiologyand diagnostics of its main problem , Eur. Phys. J.
C71 (2011) 1756,[ arXiv:1101.3270 ].[22] J. Nishimura and S. Shimasaki,
New insights into the problem with a singular driftterm in the complex Langevin method , Phys. Rev.
D92 (2015), no. 1 011501,[ arXiv:1504.08359 ].[23] K. Nagata, J. Nishimura, and S. Shimasaki,
Justification of the complex Langevinmethod with the gauge cooling procedure , PTEP (2016), no. 1 013B01,[ arXiv:1508.02377 ].[24] K. Nagata, J. Nishimura, and S. Shimasaki,
Argument for justification of the complexLangevin method and the condition for correct convergence , Phys. Rev.
D94 (2016),no. 11 114515, [ arXiv:1606.07627 ].[25] L. L. Salcedo,
Does the complex Langevin method give unbiased results? , Phys. Rev.
D94 (2016), no. 11 114505, [ arXiv:1611.06390 ].[26] G. Aarts, E. Seiler, D. Sexty, and I.-O. Stamatescu,
Complex Langevin dynamics andzeroes of the fermion determinant , JHEP (2017) 044, [ arXiv:1701.02322 ].[Erratum: JHEP01,128(2018)]. 1427] K. Nagata, J. Nishimura, and S. Shimasaki, Testing the criterion for correctconvergence in the complex Langevin method , JHEP (2018) 004,[ arXiv:1802.01876 ].[28] E. Seiler, D. Sexty, and I.-O. Stamatescu, Gauge cooling in complex Langevin forQCD with heavy quarks , Phys. Lett.
B723 (2013) 213–216, [ arXiv:1211.3709 ].[29] G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty, and I.-O. Stamatescu,
Controllingcomplex Langevin dynamics at finite density , Eur. Phys. J.
A49 (2013) 89,[ arXiv:1303.6425 ].[30] M. Scherzer, E. Seiler, D. Sexty, and I.-O. Stamatescu,
Complex Langevin andboundary terms , Phys. Rev.
D99 (2019), no. 1 014512, [ arXiv:1808.05187 ].[31] M. Scherzer, E. Seiler, D. Sexty, and I. O. Stamatescu,
Controlling complex Langevinsimulations of lattice models by boundary term analysis , Phys. Rev.
D101 (2020),no. 1 014501, [ arXiv:1910.09427 ].[32] M. Fukugita, H. Mino, M. Okawa, and A. Ukawa,
Finite size test for the finitetemperature chiral phase transition in lattice QCD , Phys. Rev. Lett. (1990)816–819.[33] Z. Fodor and S. D. Katz, A new method to study lattice QCD at finite temperatureand chemical potential , Phys. Lett.
B534 (2002) 87–92, [ hep-lat/0104001 ].[34] M. D’Elia and M.-P. Lombardo,
Finite density QCD via imaginary chemicalpotential , Phys. Rev.
D67 (2003) 014505, [ hep-lat/0209146 ].[35] M. D’Elia and M. P. Lombardo,
QCD thermodynamics from an imaginary mu(B):Results on the four flavor lattice model , Phys. Rev. D (2004) 074509,[ hep-lat/0406012 ].[36] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, Finite density QCD: A newapproach , JHEP (2004) 010, [ hep-lat/0409157 ].[37] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, Phase diagram of QCD withfour quark flavors at finite temperature and baryon density , Nucl. Phys.
B723 (2005)77–90, [ hep-lat/0503010 ].[38] Z. Fodor, S. D. Katz, and C. Schmidt,
The density of states method at non-zerochemical potential , JHEP (2007) 121, [ hep-lat/0701022 ].1539] M. D’Elia, F. Di Renzo, and M. P. Lombardo, The strongly interacting quark gluonplasma, and the critical behaviour of QCD at imaginary mu , Phys. Rev.
D76 (2007)114509, [ arXiv:0705.3814 ].[40] G. Endrodi, Z. Fodor, S. D. Katz, D. Sexty, K. K. Szabo, and C. Torok,
Applyingconstrained simulations for low temperature lattice QCD at finite baryon chemicalpotential , Phys. Rev.
D98 (2018), no. 7 074508, [ arXiv:1807.08326 ].[41] H.-S. Chen and X.-Q. Luo,
Phase diagram of QCD at finite temperature and chemicalpotential from lattice simulations with dynamical Wilson quarks , Phys. Rev.
D72 (2005) 034504, [ hep-lat/0411023 ].[42] P. de Forcrand and S. Kratochvila,
Finite density QCD with a canonical approach , Nucl. Phys. Proc. Suppl. (2006) 62–67, [ hep-lat/0602024 ].[43] A. Li, A. Alexandru, K.-F. Liu, and X. Meng,
Finite density phase transition of QCDwith N f = 4 and N f = 2 using canonical ensemble method , Phys. Rev.
D82 (2010)054502, [ arXiv:1005.4158 ].[44] S. Takeda, Y. Kuramashi, and A. Ukawa,
On the phase of quark determinant inlattice QCD with finite chemical potential , Phys. Rev.
D85 (2012) 096008,[ arXiv:1111.6363 ].[45] X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, and A. Ukawa,
Finite size scalingstudy of N f = 4 finite density QCD on the lattice , Phys. Rev.
D88 (2013), no. 9094508, [ arXiv:1307.7205 ].[46] O. Philipsen,
Constraining the QCD phase diagram at finite temperature and density ,in , 2019. arXiv:1912.04827 .[47] A. Ohnishi,
Approaches to QCD phase diagram; effective models, strong-couplinglattice QCD, and compact stars , J. Phys. Conf. Ser. (2016), no. 1 012004,[ arXiv:1512.08468 ].[48] M. Hirasawa, A. Matsumoto, J. Nishimura, and A. Yosprakob,
Complex Langevinanalysis of 2D U(1) gauge theory on a torus with a θ term , arXiv:2004.13982 .[49] R. Sommer, A new way to set the energy scale in lattice gauge theories and itsapplications to the static force and alpha-s in SU(2) Yang-Mills theory , Nucl. Phys.
B411 (1994) 839–854, [ hep-lat/9310022 ].1650] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu,
Adaptive stepsize andinstabilities in complex Langevin dynamics , Phys. Lett.
B687 (2010) 154–159,[ arXiv:0912.0617 ].[51] J. Engels, R. Joswig, F. Karsch, E. Laermann, M. Lutgemeier, and B. Petersson,
Thermodynamics of four flavor QCD with improved staggered fermions , Phys. Lett. B (1997) 210–216, [ hep-lat/9612018 ].[52] M. Scherzer,
Phase transitions in lattice gauge theories: From the numerical signproblem to machine learning . PhD thesis, U. Heidelberg (main), 2019.[53] H. Matsuoka and M. Stone,
Thermal distribution functions and finite size effects forlattice fermions , Phys. Lett. (1984) 204–208.[54] S. Hands, T. J. Hollowood, and J. C. Myers,
QCD with chemical potential in a smallhyperspherical box , JHEP (2010) 086, [ arXiv:1003.5813 ].[55] T. Yokota, Y. Asano, Y. Ito, H. Matsufuru, Y. Namekawa, J. Nishimura,A. Tsuchiya, and S. Tsutsui, work in progress.[56] S. Hands, T. J. Hollowood, and J. C. Myers, Numerical study of the two colorattoworld , JHEP (2010) 057, [ arXiv:1010.0790 ].[57] Y. Ito and J. Nishimura, The complex Langevin analysis of spontaneous symmetrybreaking induced by complex fermion determinant , JHEP (2016) 009,[ arXiv:1609.04501 ].[58] P. Amore, M. C. Birse, J. A. McGovern, and N. R. Walet, Color superconductivity infinite systems , Phys. Rev. D (2002) 074005, [ hep-ph/0110267hep-ph/0110267