Gauge cooling for the singular-drift problem in the complex Langevin method — a test in Random Matrix Theory for finite density QCD

Recently, the complex Langevin method has been applied successfully to finite density QCD either in the deconfinement phase or in the heavy dense limit with the aid of a new technique called the gauge cooling. In the confinement phase with light quarks, however, convergence to wrong limits occurs due to the singularity in the drift term caused by small eigenvalues of the Dirac operator including the mass term. We propose that this singular-drift problem should also be overcome by the gauge cooling with different criteria for choosing the complexified gauge transformation. The idea is tested in chiral Random Matrix Theory for finite density QCD, where exact results are reproduced at zero temperature with light quarks. It is shown that the gauge cooling indeed changes drastically the eigenvalue distribution of the Dirac operator measured during the Langevin process. Despite its non-holomorphic nature, this eigenvalue distribution has a universal diverging behavior at the origin in the chiral limit due to a generalized Banks-Casher relation as we confirm explicitly.


Introduction
Investigating the QCD phase diagram at finite temperature and finite density is an important subject in theoretical physics since it will provide us with microscopic understanding of various types of matter including nuclear matter. The standard Monte Carlo simulation, however, cannot be applied due to the sign problem, which is caused by the complex-valued fermion determinant in the presence of the quark chemical potential. The sign problem is also of general interest since it appears in many important cases including investigations of Chern-Simons gauge theory, supersymmetric gauge theories and the real-time dynamics of quantum systems. While there are several attempts within the framework of conventional Monte Carlo simulation such as reweighting, Taylor expansion, analytic continuation from the imaginary chemical potential and the canonical approach, they all suffer from increasing uncertainties with the increasing chemical potential. (See refs. [1,2] for reviews.) Recently, there was remarkable progress in a new type of approach based on complexifying the dynamical variables. In particular, the complex Langevin method (CLM) [3,4] was not only developed theoretically [5][6][7][8][9][10] but also shown to work practically in various interesting cases [11][12][13][14][15]. As a closely related method, Monte Carlo simulation on the Lefshetz thimble [16][17][18][19][20][21][22] has also been studied intensively. Comparison of the two methods is considered useful in deepening our understanding in this approach and in improving it [23][24][25].
The CLM [3,4] is an attempt to solve quantum systems with a complex-valued action using the idea of stochastic quantization based on the Langevin equation [26]. Since the JHEP07(2016)073 stochastic quantization does not rely on the probability interpretation of the Boltzmann weight, the method has a chance to be applicable to the case of complex action. In the case of real action, the Langevin equation is shown to converge to the correct limit using the Fokker-Planck equation (See ref. [27], for instance.). In the case of complex action, however, there are some subtleties in the proof of the correct convergence. Indeed, the CLM works well in some cases but fails in the other cases [28][29][30].
In refs. [6,7], the conditions for the CLM to work have been clarified. In solving the Langevin equation for the complex action, dynamical variables are necessarily complexified due to the complex drift term derived from the action. It is important here that the observable as well as the drift term is extended to complexified variables in a holomorphic fashion. The expectation value of the observable is defined with the probability distribution of the complexified variables. Then, under certain conditions, the expectation value thus defined is shown to be equal to the expectation value defined in the original theory for real variables with the complex weight.
One of the problems that can arise in this method is the runaway problem, which refers to the instability of the simulation. This problem, however, was shown to be avoided by using an adaptive step-size (rather than a fixed one) for the discretized Langevin time [5]. A more serious problem is the convergence to wrong limits. In ref. [6,7], one of the causes of this problem was identified to be the insufficient falloff of the probability distribution in the imaginary direction, which causes the violation of the aforementioned conditions. We call it the excursion problem in this paper since it occurs when the simulation makes a long excursion into the deeply imaginary regime.
In order to cure the excursion problem, ref. [8] proposed the gauge cooling, which enabled the application of the CLM to certain parameter regions of QCD such as the heavy dense limit [8,11] and the deconfined phase [12,13]. In our previous paper [10], we provided explicit justification of the gauge cooling based on the argument given in refs. [6,7]. We wrote down how the probability distribution of the complexified variables evolves in the Langevin time under the influence of the gauge cooling procedure. Then the corresponding Fokker-Planck equation for the complex weight was shown not to be affected by the gauge cooling as long as the observables are restricted to gauge invariant ones.
Whether the CLM works also in the confined phase with light quarks is still an open question, though. This issue was addressed using chiral Random Matrix Theory (cRMT) for finite density QCD at zero temperature [31,32], where a straightforward application of the CLM led to wrong convergence when the quark mass becomes small [33]. It was speculated that the problem is caused by the branch cut associated with logarithmic singularity in the action, which appears from the fermion determinant [33][34][35]. It was also pointed out [36] that the singular drift term derived from the logarithmic action spoils the holomorphy, which is crucial in the CLM. Recently, two of us (J.N. and S.S.) showed that the problem is not restricted to the case with logarithmic singularities in the action but it occurs also in the case with higher order singularities [9]. There it was also realized that the problem occurs in general when the probability distribution of the complexified variables is nonnegligible for configurations close to the singularity of the drift term. Therefore, we refer to this problem as the singular-drift problem.

JHEP07(2016)073
According to the new insights gained above, the singular-drift problem is similar to the excursion problem in that they are both related to the property of the probability distribution of the complexified variables. Therefore, it is conceivable that the gauge cooling, which is useful in overcoming the excursion problem, is also useful in overcoming the singular-drift problem if we choose the complexified gauge transformation used in the gauge cooling procedure appropriately.
In this paper, we test this idea in the cRMT, in which the singular-drift problem was shown to occur for light quarks [33]. While the cRMT is not a gauge theory, the action and the observables have U(N ) symmetries, which become GL(N ,C) symmetries upon complexifying the dynamical variables. These complexified symmetries provide us with large enough freedom to modify the probability distribution of the complexified variables in such a way that the singularity is avoided. In the original gauge cooling for the excursion problem, the complexified gauge transformation is chosen to reduce the so-called unitarity norm, which provides an estimate on the deviation of link variables from SU(3) matrices. On the other hand, the singular-drift problem occurs when the Dirac operator (including the mass term) defined for the complexified dynamical variables has eigenvalues close to zero. Therefore, we introduce new types of norm, which are sensitive to the properties of the Dirac operator. Using the gauge cooling with the new types of norm, we show that the eigenvalue distribution of the Dirac operator measured during the simulation is indeed modified in such a way that the singularity of the drift term is avoided unless the quark mass becomes too small. Thus, exact results are nicely reproduced for quark mass, which is much smaller than that achieved by the CLM without gauge cooling.
In fact, the eigenvalue distribution of the Dirac operator for the complexified variables is invariant under GL(N ,C) transformations. However, the Langevin time evolution itself is affected nontrivially by the gauge cooling procedure due to the noise term, which does not transform covariantly under GL(N ,C) transformations. As a result, the eigenvalue distribution measured during the Langevin process is modified by the gauge cooling. Interestingly, the eigenvalue distribution obtained in the thermal equilibrium of the Langevin process depends on the choice of norm for the gauge cooling even in the cases where the exact results are reproduced. This is possible because the eigenvalue distribution under discussion is not a holomorphic quantity and hence has no direct connection to the corresponding quantity in the original path integral with the complex weight.
On the other hand, the chiral condensate, which is a holomorphic quantity, can be expressed in terms of the eigenvalue distribution [37]. This implies that the non-universal eigenvalue distribution still has a universal property, which does not depend on the norm used for gauge cooling as long as the CLM is working. In particular, we derive a generalized Banks-Casher relation, which implies that the eigenvalue distribution has a universal diverging behavior at the origin in the chiral limit, and show that it is indeed satisfied in the cRMT. As the original Banks-Casher relation is useful in QCD at zero chemical potential, the generalized version is expected to be useful in applying the CLM to QCD at finite density.
The rest of this paper is organized as follows. In section 2, we review some basic features of the CLM and its application to the cRMT. In section 3, we discuss the singular-JHEP07(2016)073 drift problem, which occurs in the CLM when the quark mass is small. We also discuss how this problem can be overcome by the gauge cooling using new types of norm. In section 4, we present our results and show that the exact results of the cRMT can be reproduced for quark mass much smaller than that achievable without gauge cooling. We also present the eigenvalue distribution of the Dirac operator measured during the Langevin process, and show that the singularity of the drift term is avoided in differently ways depending on the norm adopted. In section 5, we derive a generalized Banks-Casher relation, which connects the chiral condensate in the chiral limit to the asymptotic behavior of the eigenvalue distribution of the Dirac operator at the origin. We show that the relation is indeed satisfied in the cRMT, and discuss its implication in the context of this work. Section 6 is devoted to a summary and discussions. Some preliminary results of this work have already been presented in a proceeding contribution [38].

The method and the model
In this section, we briefly review some basic features of the CLM and discuss its application to the cRMT.

The complex Langevin method (CLM)
Let us consider a system of n real variables x i , (i = 1, · · · , n) defined by the partition function (2.1) The basic idea of the stochastic quantization is to investigate this system using the Langevin equation [26] dx where τ is a fictitious time dubbed the Langevin time, and η i (τ ) is a Gaussian random noise normalized by η i (τ )η j (τ ) = 2δ ij δ(τ − τ ). When the action S({x i }) is real, one can show that the average of an observable over sufficiently long Langevin time agrees with the expectation value of the observable in the original path integral [27]; namely where τ 0 represents the Langevin time necessary for thermalization and T represents the total Langevin time used for averaging. Since the method does not rely on the probability interpretation of the factor e −S({x i }) , it has a chance to be generalized to the case of complex action [3,4]. When the action is complex, however, the drift term − ∂S ∂x i in (2.2) becomes complex, and one has to allow the dynamical variables to take complex values as x i ∈ R → z i ∈ C during the Langevin process. The Langevin equation (2.2) should then be replaced by

JHEP07(2016)073
where the action S is now regarded as a holomorphic function of z i defined through analytic continuation. The expectation value can be calculated by , which is also defined through analytic continuation. Unlike the case of real action, there are certain conditions that should be met for the method to work in the case of complex action [6,7]. In particular, it is required that the probability distribution of the complexified variables should have a fast fall-off in the imaginary direction. For this reason, it is considered better to leave the noise term in (2.4) real [6] although it can be generalized to complex values in principle.

Chiral Random Matrix Theory (cRMT)
We consider the cRMT for N f quarks with degenerate mass m > 0 at zero temperature and finite chemical potential µ [31,32]. The partition function is defined by [31] where Φ k (k = 1, 2) are general N ×(N +ν) complex matrices. 1 The integer ν represents the topological index, which gives the number of exact zero eigenvalues of the Dirac operator D given by for later convenience. The bosonic action S b in (2.5) is given by The effective action of the system can be written as Strictly speaking, the logarithm in (2.11) has an ambiguity due to the branch cut [33,34]. This is not an issue, however, since the CLM can be formulated in terms of the weight w = [det(D + m)] N f e −S b without ever having to refer to the effective action (2.11) as was JHEP07 (2016)073 pointed out in ref. [9]. With that in mind, we keep on using the effective action (2.11) just for simplicity of terminology.
As observables in this model, one can consider the chiral condensate Σ and the baryon number density n B defined by The partition function of the cRMT can be calculated analytically using the orthogonal polynomial method [32]. It turns out that the partition function is independent of the chemical potential µ, and hence the baryon number density is exactly zero and the chiral condensate has no µ dependence. (See, e.g., ref. [33] for an explicit expression of the chiral condensate Σ suitable for numerical evaluation.) In the gauge cooling, symmetries of the system play a crucial role. Let us consider a transformation which leaves the bosonic action invariant. In order for (2.9) to be satisfied for the transformed configuration, 2 we need to have g ∈ U(N ) and h ∈ U(N + ν). The matrices X and Y in (2.7) and (2.8) transform in the same way as Φ k and Ψ k , respectively, and hence the Dirac operator D defined by (2.6) transforms as (2.14) Therefore, both the effective action (2.11) and the observables (2.12) are invariant under the transformation (2.13). Note that the Dirac operator D satisfies for any µ, which implies that the nonzero eigenvalues of D appear in pairs with opposite signs. Furthermore, when µ = 0, the Dirac operator D is anti-Hermitian D † = −D and hence its eigenvalues are purely imaginary. Thus, one can show that det(D + m) is real positive in this case. On the other hand, when µ = 0, D is no longer anti-Hermitian, and its eigenvalues take general complex values. The determinant det(D +m) becomes complex in general, which causes the sign problem.

Application of the CLM to the cRMT
In order to apply the CLM to this system, we introduce real variables (x k ) ij and (y k ) ij through taking account of (2.9). The effective action (2.11) and the observables (2.12) are functions of these real variables. Then we complexify the variables as ( redefining the action and the observables as holomorphic functions of (z k ) ij and (w k ) ij by analytic continuation. The complex Langevin equation (2.2) can be derived in the present case in a straightforward manner [33]. After complexifying the variables, the effective action and the observables may be regarded as holomorphic functions of (Φ k ) ij and (Ψ k ) ji , where the constraint (2.9) is no longer imposed. Therefore, they are now invariant under (2.13) with g ∈ GL(N ) and h ∈ GL(N + ν); namely the symmetry is doubly enhanced. We use this enhanced symmetry for the gauge cooling.

Gauge cooling for the singular-drift problem
The gauge cooling was originally proposed to cure the excursion problem in the CLM for finite density QCD [8,12]. The basic idea is to reduce the unitarity norm, which measures the deviation of the link variables from SU(3) matrices, by making a complexified gauge transformation after each Langevin step when one solves the discretized complex Langevin equation. While the gauge cooling procedure changes the Langevin dynamics nontrivially, the expectation values of gauge invariant observables are expected to be obtained correctly once the problem is cured. Recently, this has been proved explicitly [10] by extending the argument for justification of the CLM [6,7] to the case with the gauge cooling procedure. Let us recall here that the drift term involves Tr[(D + m) −1 ∂(D + m)], which becomes singular when the operator D + m has zero eigenvalues. For µ = 0, the CLM reduces to the real Langevin simulation, and D becomes an anti-Hermitian matrix, whose eigenvalues are purely imaginary. In this case, all the eigenvalues of D + m lie on a straight line parallel to the imaginary axis on the complex plane. As µ is increased, the eigenvalue distribution is broadened in the direction of the real axis, and the singular-drift problem occurs for µ larger than some critical value depending on m [33]. We propose to avoid this problem by using the gauge cooling with the complexified symmetry transformation chosen in such a way that the eigenvalues of D + m do not appear near the origin.

New types of norm for the gauge cooling
As we mentioned in section 2.3, the symmetry (2.13) of the effective action and the observables in the cRMT enhances from U(N ) × U(N + ν) to GL(N, C) × GL(N + ν, C) upon complexification of the dynamical variables. We use this complexified symmetry to perform the gauge cooling. The transformation (2.13) to be made after each Langevin step is determined by minimizing some positive definite quantity "norm", which is invariant under the U(N ) × U(N + ν) transformation but not under GL(N, C) × GL(N + ν, C) transformation.

JHEP07(2016)073
First, let us consider the "Hermiticity norm" which measures the violation of the relation (2.9). Reducing this norm has an effect of keeping the complexified variables closer to real values, and it may be regarded as an analog of the unitarity norm in QCD [8,12]. It turns out, however, that the gauge cooling with the Hermiticity norm does not help in curing the singular-drift problem. This is understandable since the eigenvalue distribution of D + m is not directly related to the property (2.9). In order to cure the singular-drift problem, we therefore need to consider other types of norm, whose reduction affects the eigenvalue distribution directly.
Here we propose two types of norm, which can be used for the gauge cooling to cure the singular-drift problem. The first one is given by where X and Y are defined by (2.7) and (2.8), respectively, and the dagger in Y † implies taking the Hermitian conjugate of Y after complexifying the dynamical variables. Since this norm vanishes if and only if D is anti-Hermitian, it provides a measure of the deviation of D from an anti-Hermitian matrix. Reducing this norm has an effect of making the eigenvalue distribution of D + m narrower in the real direction, and hence it is expected to cure the singular-drift problem. The second one is defined by Reducing this norm suppresses the appearance of small α a and tries to achieve α a > 1/ξ. Since α a > 1/ξ implies |λ a | 2 > 1/ξ, where λ a are the eigenvalues of M , the appearance of λ a close to zero is also suppressed. In actual simulations, the excursion problem and the singular-drift problem may occur at the same time. In that case, we take a linear combination of the Hermiticity norm and one of the new types of norm as where s (0 ≤ s ≤ 1) is a tunable parameter.

JHEP07(2016)073
Note that the eigenvalue distribution of the Dirac operator D is invariant under GL(N, C) × GL(N + ν, C). However, the Langevin dynamics is modified nontrivially 3 since the noise term respects only U(N ) × U(N + ν). As a result, the eigenvalue distribution for the configuration obtained in the next Langevin step is affected nontrivially by the gauge cooling even if one averages over the Gaussian noise.

Details of the gauge cooling procedure
Below we explain the gauge cooling procedure in more detail. After each Langevin step, we perform the transformation (2.13), where the transformation matrices g ∈ GL(N ) and h ∈ GL(N + ν) are determined in such a way that a given norm is reduced efficiently. Following the original proposal [8], we first calculate the gradient of the norm with respect to the complexified transformation. This can be done by considering the infinitesimal transformation where λ a and ρ a are the basis of N × N and (N + ν) × (N + ν) Hermitian matrices, respectively, normalized as tr(λ a λ b ) = tr(ρ a ρ b ) = δ ab . Note that (3.6) corresponds to an infinitesimal U(N ) × U(N + ν) transformation when a and δ a are purely imaginary.
Considering that the norm is invariant under U(N ) × U(N + ν), we assume that a and δ a are real in what follows. Under the infinitesimal transformation (3.6), Φ k and Ψ k transform as We can calculate the change of the norm as neglecting higher order terms in a and δ a , where G a , H a ∈ R represent the gradient of the norm. We may reduce the norm most efficiently by choosing ( a , δ a ) ∝ −(G a , H a ) at the linearized level (3.6). As a finite transformation, we consider 4 g = exp(−αG) , G = G a λ a , (3.10) where the real parameter α is determined in such a way that the norm for the transformed configuration of Φ k and Ψ k is minimized. (In practice, this minimization is done approximately because we can calculate the norm only for a finite number of α.) The above procedure is repeated until the norm is more or less minimized with respect to the GL(N, C) × GL(N + ν, C) transformation.

Results of the CLM with or without gauge cooling
In this section we present our results of the CLM for the cRMT. The Langevin simulation is performed 5 with the step-size = 5 × 10 −5 . We discard the first 20000 steps for thermalization; i.e., τ 0 = 1 in eq. (2.3). Then we perform 80000 steps, which corresponds to T = 4 in eq. (2.3), during which we measure the observables every 200 steps. After each Langevin step, we perform the gauge cooling, which amounts to making a GL(N, C) × GL(N + ν, C) transformation (3.10) and (3.11) ten times as described in the previous section. As for the parameter s in eq. (3.5), we choose s = 0 for the norm N 1 (s) and s = 0.01 for the norm N 2 (s). The parameters in the definition (3.4) of the norm N 2 are chosen as ξ = 300 and n ev = 2.
In figure 1 we plot our results for the chiral condensate and the baryon number density obtained with or without gauge cooling againstm ≡ mN . Here we set the parameters of the cRMT as ν = 0, N f = 2, N = 30 andμ ≡ µ √ N = 2, which were used in ref. [33] to reveal the problem at small quark mass. 6 Our result obtained without gauge cooling is consistent with the one obtained in ref. [33]. As was pointed out there, the exact result is reproduced only form 10 and the result form 10 turns out to be close to the result of the phase-quenched model, in which the fermion determinant is replaced by its absolute value. On the other hand, the result obtained with the gauge cooling using the norm N 1 agrees well with the exact result all the way down tom ∼ 1. The result obtained with the norm N 2 is equally good except for the data points atm ≤ 1, which exhibit certain deviation.
The above results suggest that the gauge cooling with the new types of norm can solve the singular-drift problem of the CLM, which occurs at small quark mass. In order to see it 5 Although we used a fixed step-size during the simulation instead of an adaptive one [5], we did not encounter a runaway problem. 6 The cRMT (2.5) is equivalent to the chiral perturbation theory (the low energy effective theory of QCD) in the -domain in the large-N limit with the parametersm ≡ mN andμ ≡ µ √ N fixed [39], which is called the microscopic limit. In this section we usem andμ just to make it easier to compare our results with those in ref. [33]. When we discuss the generalized Banks-Casher relation in section 5, however, we have to take the large-N limit with m and µ fixed, which is different from the microscopic limit.  Note that the eigenvalues (±λ) + m appear in pairs in all the three cases for the reason given below eq. (2.15). In the case without gauge cooling, the eigenvalue distribution covers the singularity at the origin, which implies that the singular-drift problem occurs. However, the eigenvalue distribution is changed drastically by the gauge cooling with the new types of norm. The norm N 1 makes the eigenvalue distribution narrower in the real direction, whereas the norm N 2 makes the eigenvalues repelled from the domain near the origin. Thus, we find that the gauge cooling with the new types of norm indeed has an effect of removing the eigenvalues near the singularity. In order to see this effect more quantitatively, let us consider the radial distribution of eigenvalues [9] defined by ρ(r) = 1 2πr dxdy ρ (CL) (x, y) δ( (x + m) 2 + y 2 − r) , part of the eigenvalue, respectively. (See (5.1) for the precise definition of ρ (CL) (x, y).) In figure 2 (d), we plot the results for the radial distribution in the three cases discussed above. We find that the radial distribution is strongly suppressed near the singularity by the effect of the gauge cooling.
Let us then discuss what happens if we increase the chemical potential µ further. The exact result for the cRMT is known to be independent of µ as we mentioned below (2.12). When the gauge cooling is not performed, the eigenvalue distribution of D + m measured during the complex Langevin simulation tends to become wider in the real direction for increasing µ. The question is whether the gauge cooling with the new types of norm can remove the eigenvalues of D + m near the origin even in that case. In order to answer this question, we measure the width ∆ of the eigenvalue distribution of D on the real axis. As long as ∆/2 < m, the eigenvalues of D + m do not appear near the origin and the CLM works. In figure 3 we plot ∆/2 obtained by the CLM withm = 5 andμ = 1, 2, 3, 4 with or without gauge cooling. We find that the CLM without gauge cooling works only for µ 1.2, whereas the gauge cooling with the norm N 1 and N 2 makes the CLM work for µ 3.3 andμ 2.2, respectively. We should recall, however, that the norm N 2 has two parameters ξ and n ev , which can be optimized to make the gauge cooling more efficient. Also, the range of applicability may be further enlarged by increasing the number of gauge cooling procedures after each Langevin step or by using a smaller Langevin step-size.

The generalized Banks-Casher relation
In the previous section, we have seen that the gauge cooling with the new types of norm cures the singular-drift problem caused by the eigenvalues of D + m near the origin. Interestingly, we find that the eigenvalue distribution of the Dirac operator measured during the complex Langevin simulation depends on the choice of the norm even in the parameter region where the singularity is avoided. This is possible because the eigenvalue distribution JHEP07(2016)073 of the Dirac operator is non-holomorphic, and hence the distribution measured during the complex Langevin simulation is not directly related to the corresponding quantity defined in the original path integral. (Note, for instance, that the former is real non-negative by definition, whereas the latter is complex in general due to the complex weight.) This is in contrast to the observable such as the chiral condensate, which is holomorphic, and hence the quantity measured during the complex Langevin simulation should agree with the corresponding quantity defined in the original path integral.
In fact, the chiral condensate can be written in terms of the eigenvalue distribution of the Dirac operator. This leads to the well-known Banks-Casher relation [40] at µ = 0, which expresses the chiral condensate in the chiral limit in terms of the eigenvalue distribution at the origin. At µ = 0, the eigenvalue distribution of the Dirac operator defined in the path integral becomes complex due to the complex weight, and the violent oscillation of the complex eigenvalue distribution in an extended region is responsible for the nonzero chiral condensate in the chiral limit [39]. On the other hand, the eigenvalue distribution of the Dirac operator measured in the complex Langevin simulation is real non-negative, and therefore its relation to the chiral condensate should be analogous to the original one at µ = 0. However, since the eigenvalue distribution spreads out in the complex plane due to the violation of the anti-Hermiticity of D, the accumulation of eigenvalues toward the origin should occur in the chiral limit in order to reproduce the non-zero chiral condensate [37]. What we have seen in the previous section is qualitatively consistent with this argument. In what follows, we make this argument quantitative by writing down the generalized Banks-Casher relation, which relates the chiral condensate with the diverging behavior of the eigenvalue distribution at the origin in the chiral limit. We confirm this relation numerically in the cRMT, and discuss its implication in the context of this work.
First let us define the eigenvalue distribution of the Dirac operator measured in the CLM by where λ i are the eigenvalues of the Dirac operator D, n is the number of the eigenvalues, and · · · CL represents an ensemble average over the thermalized configurations generated by the complex Langevin simulation. Note that the quantity (5.1) for fixed x and y (before taking the ensemble average) cannot be a holomorphic function of the configuration (Φ k , Ψ k ) since its imaginary part vanishes identically while its real part is not a constant. Therefore, the distribution (5.1) does not have to agree with the eigenvalue distribution of the Dirac operator defined in the original path integral. In fact, the latter distribution is generally a complex-valued function due to the complex weight of the path integral. As is done in deriving the original Banks-Casher relation at µ = 0, we express the chiral condensate 7 using the eigenvalue distribution ρ (CL) (x, y) as [37]

JHEP07(2016)073
Note that even after the complexification of the dynamical variables, the Dirac operator D satisfies (2.15), which leads to ρ (CL) (x, y) = ρ (CL) (−x, −y). Therefore, the eigenvalues with |x + iy| m cancel pairwise in (5.2), and hence small eigenvalues with |x + iy| m are responsible for the nonzero condensate.
Let us introduce the polar coordinates x = r cos θ and y = r sin θ, where −∞ < r < ∞ and 0 ≤ θ < π, and the corresponding distributionρ (CL) (r, θ), which has the symmetryρ (CL) (r, θ) =ρ (CL) (−r, θ). Using this symmetry, we can rewrite (5.2) with the polar coordinates as In the m → +0 limit, the terms in the parentheses can be replaced with 2πie −iθ δ(r). Here and henceforth, the chiral limit represented by lim m→+0 is assumed to be taken after taking the thermodynamic limit (n → ∞).
To proceed further, let us recall that the chiral condensate (5.2) obtained by the CLM is guaranteed to be real by symmetry. Note that the complex Langevin equations are invariant under complex conjugation (Φ k , Ψ k ) → (Φ * k , Ψ * k ), under which the Dirac operator becomes D → D * . This ensures that the eigenvalue distribution in the thermal equilibrium of the Langevin process has the symmetry ρ (CL) (x, y) = ρ (CL) (x, −y), from which it follows that (5.2) is real. In the case of polar coordinates, the symmetry translates toρ (CL) (r, θ) =ρ (CL) (−r, π −θ), which implies that (5.4) is real and leads to the generalized Banks-Casher relation lim m→+0 Σ = π lim The relation (5.5) implies that the quantity lim m→+0ρ (CL) (r) should diverge as when lim m→+0 Σ is nonzero; namely when the chiral symmetry is spontaneously broken. When the singular-drift problem occurs, the chiral condensate obtained by the CLM becomes smaller than the exact result ( figure 1 (left)). The gauge cooling with appropriate norm cures the singular-drift problem by making the eigenvalue distribution of D + m suppressed strongly near the singularity ( figure 2 (b) and (c)). This has an effect of making the eigenvalue distribution of D (not D + m) larger near the origin, and hence the chiral condensate becomes larger. Note that the CLM may work for different choices of norm as we have shown in the cRMT. In that case, the eigenvalue distribution of D itself may depend on JHEP07(2016)073   the choice of norm as we emphasized above, but the asymptotic behavior (5.7) ofρ (CL) (r) in the chiral limit is universal since it is fixed by the generalized Banks-Casher relation (5.5).
Let us confirm the generalized Banks-Casher relation (5.5) in the cRMT. Comparing (2.12) with (5.2), where n = 2N , we obtain the corresponding relation in the cRMT as . For N = 10, we double the statistics to reduce the statistical errors. 8 The gauge cooling is performed either with the norm N 1 (s) or N 2 (s). 9 For either choice of the norm, we observe a clear plateau, which extends towards the origin as N increases. This is consistent with our expectation that the quantityρ (CL) (r) diverges as ∼ r −1 for r → 0 in the chiral limit. The height of the plateau, which is almost independent of N , gives an estimate for the r.h.s. of (5.8), which is close to the exact result lim m→+0 Σ = 4 for the cRMT in the chiral limit. Certain deviation depending on the norm adopted can be understood as finite m effects.

Summary and discussions
In this paper we proposed a new method to overcome the singular-drift problem in the CLM, which occurs, for instance, in finite density QCD at small quark mass in the confined phase. This problem occurs when the drift term in the complex Langevin equation has JHEP07(2016)073 singularities, and the probability of obtaining configurations near the singularities during the Langevin process is not suppressed sufficiently. In the case of finite density QCD, the drift term becomes singular when the Dirac operator including the mass term has zero eigenvalues.
Our method is based on the gauge cooling, which was originally proposed to overcome the excursion problem. Unlike the original proposal, however, we use new types of norm which are sensitive to the eigenvalue distribution of the Dirac operator. Performing the gauge cooling with the new types of norm after each Langevin step, we can remove the eigenvalues close to zero in such a way that the calculation of gauge invariant observables is not affected. We tested the method in the cRMT, which is a simplified model of finite density QCD at zero temperature, and confirmed that the method allows us to reproduce the exact result even in the small quark mass regime, which was not accessible without gauge cooling due to the singular-drift problem.
We emphasized that the eigenvalue distribution of the Dirac operator measured during the Langevin process is not a holomorphic quantity, and hence it does not have a unique counterpart in the original path integral with the complex weight. We derived the generalized Banks-Casher relation, which expresses the chiral condensate in terms of this non-universal eigenvalue distribution in the chiral limit, and confirmed it explicitly in the cRMT. While the eigenvalue distribution for the two successful choices of norm is different, the asymptotic behavior at the origin, which is linked to the chiral condensate, is shown to be universal albeit with certain finite N and finite m effects.
Our method can be applied to QCD in a straightforward manner. In particular, the norm (3.4) can be used for any lattice Dirac operator. It is technically important here that one only needs low-lying eigenvalues of a large Hermitian matrix, which can be obtained efficiently by the Lanczos method. On the other hand, the norm (3.3) can be used for the staggered fermion since the corresponding Dirac operator has the property D(µ) † = −D(−µ) as in the cRMT, but it cannot be used for the Wilson fermion. The generalized Banks-Casher relation can also be used for the staggered fermion but not for the Wilson fermion. It would be interesting to see to what extent our method can enlarge the range of applicability of the CLM in finite density QCD.