Complex Langevin Simulation of a Random Matrix Model at Nonzero Chemical Potential

In this paper we test the complex Langevin algorithm for numerical simulations of a random matrix model of QCD with a first order phase transition to a phase of finite baryon density. We observe that a naive implementation of the algorithm leads to phase quenched results, which were also derived analytically in this article. We test several fixes for the convergence issues of the algorithm, in particular the method of gauge cooling, the shifted representation, the deformation technique and reweighted complex Langevin, but only the latter method reproduces the correct analytical results in the region where the quark mass is inside the domain of the eigenvalues. In order to shed more light on the issues of the methods we also apply them to a similar random matrix model with a milder sign problem and no phase transition, and in that case gauge cooling cooling solves the convergence problems as was shown before in the literature.


Introduction
A first principles study of the QCD phase diagram in the plane of temperature (T ) and baryon chemical potential (µ) is one of the most challenging problems of modern high energy physics. Its understanding will lead to profound answers ranging from cosmology and the early universe to the physics of neutron stars. Analytical approaches tend to fail because the theory is strongly interacting, and only for extreme values of the temperature and/or the baryon chemical potential can the theory be studied perturbatively due to asymptotic freedom. Lattice numerical simulations have contributed tremendously to the understanding of the vacuum properties of the theory and have also firmly established that with physical quark masses the deconfinement transition at zero baryon chemical potential is a crossover. Despite all these celebrated results the situation at finite baryon density is very different [1] and our knowledge is mainly based on models for QCD or lattice simulations at small values of µ (more precisely µ/T < 1 and µ < m π /2 [2]).
It is well known that the culprit behind this lack of results is the infamous sign problem which is prohibiting numerical simulations when µ/T > 1 or µ > m π /2. The determinant of the Dirac operator becomes complex for SU(N c ) Yang-Mills theories with N c ≥ 3 and quarks in the fundamental representation. Consequently, standard Markov Chain Monte Carlo (MCMC) methods, which require a real and positive probability weight, cannot be applied. There have been many attempts at tackling this problem by the QCD community. Some of them try to circumvent the sign problem, while others study related theories which have no sign problem. To circumvent the sign problem one can perform a Taylor expansion around µ = 0 [3] or use reweighting methods [4], however, these methods cannot go beyond µ/T ≥ 1 at physical quark masses due to serious problems such as the limited radius of convergence of the Taylor series or the exponentially small reweighting factor. Alternatively, one can study QCD with imaginary baryon chemical potential [5,6], or perform simulations of two-color QCD or of QCD with adjoint quarks [3,7,8], which have no sign problem at all. However, as these theories have a different phase diagram from the one of QCD, one can at best extract qualitative information regarding the QCD phase diagram.
A method that has attracted a great deal of attention recently, and which is not based on MCMC methods, is the method of stochastic quantization, also called Langevin method. For the case of complex actions, the complex Langevin (CL) method was pioneered independently by Parisi [9] and Klauder [10] more than 30 years ago. Despite the fact that stochastic quantization yields the same results as path integral quantization for systems with a real action, this is, unfortunately, not always the case when the action is complex. One of the most serious problems is that the method sometimes converges towards the wrong limit. Convergence criteria have been established [11], however, these are not fulfilled in realistic QCD simulations, at least for the range of parameters that are of interest for mapping the unknown part of the QCD phase diagram [12,13]. Nevertheless, the CL algorithm has given correct results in many non-trivial systems for which we know the solution, and it seems to be quite successful for QCD simulations in the deconfined phase [12,13], as well as in simulations for heavy quarks [14][15][16][17], where the results can be validated by other methods.
In this article we are attempting to understand the properties of the algorithm very close to the chiral limit in the cold and dense regime. To achieve that, we are studying a random matrix theory (RMT) model which shares many key features of QCD such as spontaneous breaking of chiral symmetry, a finite density phase transition, as well as a complex fermion determinant, which causes a strong sign problem. The model that we have been studying was introduced by Stephanov [18] based on a random matrix model for the finite temperature chiral phase transition [19]. There is a significant literature studying the convergence properties of the CL algorithm in RMT but all the existing studies are based on a finite density model introduced by Osborn [20] (or an improved version thereof [21]), which possesses many similarities with the one by Stephanov but also has big differences, most notably the lack of a phase transition to a nonzero baryon density phase. Sensu stricto the Osborn model is only a model of QCD in the confined phase at small chemical potential.
A great deal of analytical knowledge for non-perturbative aspects of QCD came from RMT studies. These include among others finite density results [22], lattice spacing effects on the lowest eigenvalues of the Dirac operator for Wilson fermions [23][24][25][26][27][28], and the effect of topology on the Dirac spectrum [29]. In this article we are addressing the convergence properties of the CL algorithm for a model of continuum QCD at nonzero chemical potential. The model has a known analytic solution, and by simulating it numerically we can get an explicit handle on the various issues of the algorithm.
This article starts out with the definition of the random matrix models that will be studied by the CL algorithm which is introduced in section 3. In section 4, we discuss the fermion determinant and the spectrum of the Dirac operator. The CL reweighting method is analyzed in section 5, while the shifted representation, in which the chemical potential is shifted to the bosonic part of the action, is discussed in section 6. Cooling methods are investigated for two different random matrix models in section 7. As a last attempt to fix the convergence problems of the CL algorithm, we study the deformation method in section 8. Concluding remarks are made in section 9, and analytical results for the phase quenched partition function are worked out in Appendix. A preliminary account of some of the results in this paper appeared as conference proceedings [30,31].

Random Matrix Model
In this section we discuss a random matrix theory inspired model [19,32] for QCD at finite baryon density originally proposed by Stephanov [18]. The model's partition function reads (2.1) The Dirac operator D has the form where a term containing the baryon chemical potential µγ 0 has been coupled to the chRMT Dirac operator proposed in [29,32]. The N × (N + ν) matrix elements of W are complex numbers, N is the size of the block matrix W and the index ν of the Dirac Matrix is the analogue of the topological charge. This model was first introduced for imaginary chemical potential [19] to study the QCD chiral phase transition at nonzero temperature (which in the model appears as an imaginary chemical potential). In this article we choose ν = 0, since the topological charge does not have a significant effect on the quantities of interest. For µ > 0, the eigenvalues of D become complex, and are roughly distributed homogeneously inside a strip of width ∼ µ 2 (for finite N it is an ellipse). Similarly to QCD, numerical simulations of this random matrix theory have an exponentially hard sign problem, especially when the quark mass is inside the cloud of eigenvalues of the Dirac operator.
Our attention will be focused mainly on two observables, the mass dependent chiral condensate (note that the physical chiral condensate is −Σ) defined by and the baryon number density given by There is no unique way of introducing a chemical potential in a random matrix model. Some alternatives turn to be advantageous from a symmetry point of view [20]. In particular, the Osborn model [20] has a U(N ) × U(N ) symmetry which makes it possible to obtain analytical results for the joint probability distribution function of the eigenvalues, which is the starting point of many powerful random matrix methods. The Stephanov model has only a U(N ) invariance, and it is not possible to obtain an analytical solution for the joint eigenvalue density. In the case of the Osborn model, by extending the method of orthogonal polynomials to bi-orthogonal polynomials, all n-point spectral correlators can be obtained. The Osborn model is a two-matrix model that has a form similar to the Stephanov model. In the chiral basis it is given by where the Dirac operator D has the form Remarkably, the partition function at finite baryon density can be related to the one at zero baryon density by introducing a trivial multiplicative factor and a mass rescaling as follows [20,33], Consequently, it is natural to expect that the Osborn model does not possess the rich phenomenological structure of the Stephanov model, which exhibits a phase transition separating a phase with zero baryon density from a phase with nonzero baryon density.
Strictly speaking the Osborn model should only be considered as a model for QCD at small chemical potential, precisely due to the absence of a phase transition to a phase with nonzero baryon density. In addition, one can conclude that the sign problem of the Osborn model is of a weaker nature and therefore may be remedied by some clever techniques [21,[33][34][35][36][37]. Both random matrix models possess the same global symmetries with the same spontaneous symmetry breaking pattern as in QCD and yield the epsilon limit of the QCD chiral Lagrangian. It is noteworthy that in case of QCD with three colors in the fundamental representation this chiral Lagrangian does not have a dependence on the baryon chemical potential. The reason is that the Goldstone bosons, i.e., the pions, do not carry baryon charge. The chiral Lagrangian of phase quenched QCD, where µ becomes the isospin chemical potential, has a nontrivial µ-dependence which, at the mean field level, or in the ε domain, is given by the µ dependence of the large N limit of the partition functions (2.1) or (2.5). We could contemplate other random matrix models where the dependence of the chemical potential is integrated out in the evaluation of the partition function. For example, the Dirac operator where ϕ is uniformly random in [−π, π], and the partition function is defined by It is clear that the partition function does not depend on µ while the eigenvalues of D are complex. Since the chemical potential can be eliminated by changing the integration contour of the ϕ integral, the CL algorithm should be able to solve this problem correctly. We will, however, not study this model in this paper. The unquenched partition function of the Stephanov model was cast analytically in a form that allows for either an easy numerical evaluation at finite N , or that allows for a complete analytical solution via a saddle point approximation in the thermodynamic limit where N → ∞ [18,22]. For the N f = 1 case, the partition function, in units where the chiral condensate Σ = 1, takes the following σ-model form via bosonization methods where σ is the bosonized version ofψ L ψ R . A change of variables to polar coordinates renders the angular integral calculable analytically and yields a modified Bessel function, such that the partition function can be written as a one-fold integral, A saddle point analysis of the partition function can be performed in the thermodynamic limit. This was analyzed in detail in [22] but we will repeat some of the main steps here for the convenience of the reader. The saddle point equation reads The Stephanov model exhibits a first order phase transition which takes place when |Z u=u b | = |Z u=ur |, with u b and u r being two different solutions of the saddle-point equation giving the same free-energy. One can rewrite this condition This is a transcendental equation that, in the chiral limit, has the solutions u r = 0 and u b = 1 + µ 2 . Therefore in this limit one has the critical curve Re 1 + µ 2 + log µ 2 = 0, (2.14) which for the case of a real baryon chemical potential leads to the critical value µ c = 0.527 . . . in the chiral limit. This critical curve is also valid for the case of a complex chemical potential. In particular, for an imaginary chemical potential, there is a second order phase transition to a restored phase at µ = 1. We will also compare our results to the two-flavor phase quenched RMT partition function [38] or the partition function at nonzero isospin chemical potential [39,40]. The two-flavor partition function is an eight dimensional integral and is much more complicated than the one-flavor partition function, which is only a two-dimensional integral. However, in the large N limit, it can be evaluated by a saddle-point approximation, see Appendix A. It has a pion-condensation phase for µ > m π /2 corresponding to the parameter domain when the quark mass is inside the support of the eigenvalues. In Fig. 1 we show the phase diagram in the plane of the chemical potential and quark masses (which are taken to be equal for the two flavors). For nonzero mass and increasing chemical potential, we find a phase transition to a pion condensation phase at µ = m π /2, and for larger µ, a second phase transition to a chirally restored phase. In the region between the curves, the quark mass is in the domain of the eigenvalues and CL is expected to fail. In the outside region, the mean field result for full QCD and phase quenched QCD coincide, and CL is expected to work. For QCD we expect a similar forbidden region.
In order to study the properties of the Langevin algorithm we will perform numerical simulations of the Stephanov model employing the CL algorithm and test its convergence properties by comparing the obtained numerical data for the chiral condensate and the baryon density with analytical results computed using the partition function (2.11). In several cases we will also simulate the Osborn model in order to display potential issues that might arise for the CL method when switching from a model without a phase transition to one where the sign problem triggers a phase transition.

Complex Langevin
Stochastic quantization and the Langevin equation form a natural bridge between quantum field theory (QFT) and statistical mechanics. In the case of a real action, expectation values of the path integral can be obtained by averaging over an ensemble of configurations that have been generated by the Langevin evolution. Here, in order to set the stage and define our notation, we will consider the one degree of freedom, trivial "QFT" whose partition function has the following path integral form Z = e −S(x) dx. The discretized real Langevin equation for updating the dynamical variable x is where the noise term ∆ξ is a stochastic variable with zero mean and variance given by 2 √ ∆t. Generalizing the concept of stochastic quantization to the case of complex actions, requires us to promote every real degree of freedom to its complex counterpart. This complexification will naturally occur when evolving the degrees of freedom according to the Langevin equation, as the derivative of the action, usually coined as the drift term (∂ x S(x(t))∆t), is complex and will push the dynamical variables into the complex plane. In this case, x will give its place to z = x+iy, whose evolution as a function of the Langevin time t will be given by the following update equation The Langevin equation thus generates a probabilistic ensemble {z(t)} where observables are calculated by averaging along the Langevin trajectory. One can quite easily generalize the Langevin equation from systems with one degree of freedom to more complicated systems, such as field theories with an infinite number of degrees of freedom. In our case we need to modify this formalism for the case of an RMT model, which can be done in a straightforward way as is shown below. The complex random matrix W in the random matrix model (2.2), in its Cartesian representation, can be decomposed as W = A + iB where W † = A − iB with A and B both real. In this case the measure of integration dW dW † becomes dAdB. The action corresponding to the partition function (2.1) reads At finite chemical potential the matrices A and B will take on complex values due to the complex Langevin flow. We therefore introduce the complexified matrices X = A + iB and Y = A − iB which will replace W and W † in the following expressions. At µ = 0 we have X † = Y , however this will not be the case as A and B become complex. The matrices A and B will have the following Langevin evolution where we have simplified the notation by introducing the matrix G, The first step in our simulations is to establish that in the absence of a baryon chemical potential, when the action is real, the real Langevin simulations give the correct analytical answer. This is shown in Fig. 2 where the baryon number (left) and the chiral condensate (right) are plotted as a function of the mass m. We also check the N dependence of our results to confirm that the results are independent of the matrix size. This is also important when comparing our results to those from the phase quenched analytic results, as these are computed in the large N limit. Fig. 3 shows the baryon number (left) and the chiral condensate (right) at two different matrix sizes, and we conclude that N = 48 is sufficient. Finally, we show that our results do not depend on the step size of the discretized Langevin equation. As is demonstrated in Fig. 4 for the baryon number (left) and the chiral condensate (right), a choice of ∆t = 10 −4 is sufficient to eliminate discretization errors.
Having convinced ourselves that the CL algorithm has been implemented correctly we now simulate the random matrix model at nonzero chemical potential and compare Quite surprisingly, we find that our numerical CL results agree with the analytical phase quenched results, and only see agreement with the dynamical one flavor results when these coincide with the phase quenched results. This is the case when the quark mass is outside the domain of the eigenvalues of the Dirac operator, or equivalently, when the chemical potential is outside the domain of the eigenvalues of γ 0 D. In Fig. 6 we show the baryon density and the chiral condensate as a function of the quark mass for µ = 0.2 (top row) and µ = 1 (bottom row), from which we draw similar conclusions. The shaded regions in Figs. 5 and 6 and further down in the paper denote the region where the analytical one-flavor mean field results do not agree with the phase quenched mean field results. We thus conclude that the CL algorithm fails in the region when the baryon number density and the chiral condensate are not holomorphic functions of the matrix elements. We will discuss this in more detail in the next section where we discuss the fermion determinant and the Dirac spectrum.

The Dirac Spectrum and the Fermion Determinant
One of the requirements for the correct convergence of CL is that the "operator", in our case is a holomorphic function of the complexified variables. This is not the case for the chiral condensate when the quark mass is inside the two-dimensional locus of the eigenvalues of D, or equivalently, for the baryon density if the chemical potential is inside the spectral support of γ 0 D, which is also a two-dimensional domain. So for the CL to be convergent we need to require that det(D +µγ 0 +m) = det(γ 0 (D +m)+µ) > ε with ε a finite constant. In Fig. 7 we show scatter plots of the determinant in the complex plane, obtained during the CL simulation for zero and nonzero mass, respectively. Indeed, we will find that simulations converge well when the flow of the determinant avoids the origin, as has also been observed, for example, for the Osborn RMT model [36] and for two-dimensional QCD [41]. In random matrix theory the effect of the fermion determinant on the global distribution of the eigenvalues is a 1/N -correction, and will arrange only a small number of eigenvalues near zero. Also for a nonzero imaginary chemical potential the quenched and the dynamical eigenvalue distribution are the same for large N . For real chemical potential the eigenvalue distribution is complex because of the phase of the fermion determinant, but for large N the spectral support is still the same as for the quenched or phase quenched theory. Since a fermion determinant does not change the overall spectral density to leading order in 1/N , we expect that also for real chemical potential the distribution of the eigenvalues will not be affected significantly by the CL evolution. Indeed, as can be seen in Fig. 8, where we plot the eigenvalues of the Dirac operator, for various values of the chemical potential, the support of the Dirac eigenvalues on the CL trajectory is still given by the quenched result (green curve). Since the CL algorithm is probabilistic we thus necessarily have that the chiral condensate and the baryon number, which are determined by the distribution of Dirac eigenvalues, will be given by the (phase-)quenched result, as we have seen in Figs. 5 and 6. In the next section we will show that the correct result can be obtained using a reweighting algorithm.

Reweighted Complex Langevin
After having established that the CL algorithm fails to reproduce the known analytical results of the random matrix model (2.2), the obvious question to be asked is if something can be done to fix the pathologies of the algorithm in regions of the parameter space where it fails. We apply the reweighted complex Langevin (RCL) method [42,43], and we will show that one can significantly improve the convergence properties of the algorithm. This will lead to correct results for values of the parameters for which a naive implementation of the CL algorithm was giving wrong results. Our efforts to test the algorithm will be focused on the region close to the phase transition. The motivation of the method mainly comes from the expectation that reweighting CL trajectories might work better than other traditional forms of reweighting, mainly because the target ensemble with parameters (ξ = m, µ) and the auxiliary ensemble with parameters (ξ 0 = m 0 , µ 0 ) are expected to have larger overlap. In reweighting one computes However, contrary to the traditional forms of reweighting, the weight w(x; ξ 0 ) = e −S(x;ξ 0 ) is complex, and thus, we need to employ the CL algorithm to sample this auxiliary ensemble. Of course, this is performed through a judicious selection of the reweighting parameters, chosen from the region where the algorithm satisfies the CL convergence properties. In this case the reweighting equation becomes, after complexification of the variables, where P (z; ξ 0 ) is the real probability in the complex variables z = x + iy, generated by the CL trajectory. In practice, the configurations of the auxiliary ensemble are sampled according to their probability P (z, ξ 0 ) in the complexified variables by evolving the CL equations for the auxiliary action. An expectation value in the target ensemble is then computed as a ratio of the average effective observable and the average reweighting factor, both measured along the auxiliary CL trajectories.  Figure 9. Results for the RCL method for the random matrix model (2.2) with matrix size N = 6. The mass scan for µ = 1 (upper panels) is generated from an auxiliary ensemble with m 0 = 4 at the same µ, while the µ-scan for m = 0.2 (lower panels) uses an auxiliary ensemble at µ 0 = 2 and same mass. We show the baryon number density on the left and the chiral condensate on the right.
We first test the method with small matrices (N = 6) and we see in Fig. 9 that reweighting the CL trajectories can fix all the problematic issues of the algorithm. It is interesting to observe that for this relatively small matrix size the analytical answer can be reproduced for the whole range of parameters, as can be seen from the scans of the mass and of the chemical potential. It is important to stress that this is already quite intriguing since the naive implementation of the algorithm was failing to reproduce the correct answer in the region where the operators are non-holomorphic. Moreover, it is interesting to observe that, even though the auxiliary ensemble is chosen at one side of the phase transition, i.e., with large m 0 for the mass scan or large µ 0 for the µ-scan, the RCL data at the other side of the phase transition still agree very well with the analytical results. Nevertheless, we already notice that the error bars start to grow in the phase transition region, as is expected if the sign problem grows in that region.
Of course in order to claim to have solved the sign problem, which is exponentially hard with respect to the volume of the system, we need to show that the number of matrices needed to achieve the sought precision does not scale exponentially with the matrix size N . For this reason we increased the matrix size, in order to investigate if this method of reweighting actually works and how it scales with respect to the matrix size N . In µ and the lower graphs a µ scan for fixed mass. To keep the error under control in the phase transition region the number of configurations had to be increased by a factor of 100 compared to N = 6. This allows us to have small error bars in the mass scan for µ = 1.0, for all mass values. In this scan the RCL always gives the correct value even though CL obviously fails as soon as m < 1.0. In the µ-scan we see that RCL is able to improve on the CL method outside the phase transition, however, inside this region, i.e., for 0.6 < µ < 0.8, the sign problem clearly reappears. The figure plainly shows that RCL performs qualitatively better than the CL method in all cases, but with the caveat that the phase transition region is still difficult to access. Undoubtedly, the same is true for N = 48 as can be seen in Fig. 11. The RCL works reasonably well in the mass scan, performed for µ = 1, whereas the CL does poorly over most of the mass range. Unfortunately, the µ scan distinctly shows that the reweighting only works above and below the phase transition region. Again we observe the salient feature that even though the auxiliary ensemble is taken at one side of the phase transition, the reweighting procedure reproduces the data well also at the other side of the transition. This seems to point to the absence of an overlap problem in the RCL method, even though the sign problem is clearly present in the phase transition region. One disturbing point of this investigation is that it confirms how bad the original CL performs for a very large range of parameters, even away from the phase transition. Indeed, the CL seems to fail in regions where reweighting still works quite well and is not yet hampered by the sign problem.
Using the data for N = 6, 12, 24, 48 we find a naive volume scaling for the RCL method that is proportional to exp(0.3 × N ) for the number of configurations necessary to get the same accuracy for all values of N . This shows that, even though this reweighting method rectifies the failing of the CL method for a large range of parameter values, it still is exponential in the volume in the phase transition region.

Shifted Representation
In an attempt to mend the problems due to the phase of the Dirac operator we now shift the effect of the chemical potential away from the fermionic term. This can be done with a simple shift of variables. Written out in the Cartesian representation, the Dirac operator is We can absorb µ into A with a simple change of variables, A = A − iµ. The action in terms of the matrices A and B is where X = A + iB and Y = A T − iB T . In this representation the µ dependence has been shifted from the fermionic to the bosonic term. Computing the CL force term results in where G = (m 2 + Y X ) −1 is defined in terms of the shifted fields. The advantage of the shifted representation is that it starts in an anti-Hermitian state, and due to the fact that CL is non-deterministic, the configurations could potentially evolve to a different minimum.
To analyze the dynamics of the shifted representation we analyze the elements of the matrices A and A during CL evolution; the real and imaginary part of their average diagonal entry are shown in Fig. 12. Although the two matrices start out very differently, they are similar after thermalization. This seems to indicate that and thus they converge to the same solution. Since the Dirac operator in the shifted representation starts the CL evolution at a chiral condensate and a baryon number density for µ = 0, one might expect better convergence properties at least below the critical value of the chemical potential. In the next section we will use the shifted representation when analyzing the effect of gauge cooling.

Gauge Cooling
The complexified action takes on redundant degrees of freedom which is evident from the fact that the action is invariant under an enhanced symmetry group as compared to the original RMT. One can utilize this enlarged symmetry in an attempt to steer the Langevin flow towards more physical configurations due to the fact that although the action is invariant under these transformations, the flow itself is not. This method is commonly referred to as gauge cooling, and has been used to great effect in a plethora of models [14,44]. Most relevant to our study is its successful application to the Osborn RMT model (2.6) [44], which we will refer to for comparison. The original RMT is invariant under the U(N ) transformation However the complexified action is invariant under the enlarged GL(N, C) transformation We stress that the cooling transformation does not change the eigenvalues of the Dirac operators D and γ 0 (D + m), and that the effect of cooling occurs in tandem with the Langevin updates. Next we will look at how to choose h in an advantageous way.

Cooling norms
The transformation matrices h are chosen such that a cooling norm is reduced. These norms are constructed to quantify an undesirable property of the matrix configurations. The most basic of these is the Hermiticity norm [44] which measures the deviation of the CL configuration from a valid RMT configuration. It is zero when X † = Y , and grows when the matrices A and B acquire imaginary parts. We also introduce an eigenvalue norm [44] N ev = nev i=1 e −ξγ i (7.4) where γ i are the n ev lowest eigenvalues of the positive definite matrix D † D, and ξ is a real positive parameter. This norm suppresses configurations with Dirac eigenvalues close to zero.
Finally, we will also use a generalization of the anti-Hermiticity norm, which was introduced in [44] for p = 1. The matrices ψ and ϕ are the off-diagonal elements of D. For the Stephanov model they are given by ψ = iX + µ and ϕ = iY + µ, so that the norm becomes For p = 1 the µ-dependent terms do not depend on the similarity transformation h, and the Dirac operator is generally not anti-Hermitian at the minimum of the norm. Therefore will use the p = 2 anti-Hermiticity norm below.
As the different norms try to fix different problems one can also combine them in aggregate norms. One useful choice is to combine the Hermiticity norm, which quantifies how much the configurations drift into the imaginary plane, with either the anti-Hermiticity or the eigenvalue norm, both of which handle problematic configurations related to a singular behavior of the drift (7.7)

Computing h
We follow the procedure outlined in [44] to compute the transformation matrix h. We can write h in terms of the U(N ) generators, Because the RMT is invariant under U(N ) transformations, we can choose a i ∈ R to only pick out the GL(N, C)/U(N ) transformations. Assuming the norm is a function N (X, Y ) we want to solve the following equatioñ This can be reduced to a one dimensional minimization problem by first computing the gradient descent vector of the transformation through , (7.10) and then solve the one parameter minimization problem where β is a real positive quantity. β is computed by applying Brent's method [45], to which we add an upper bound to avoid a numerically unstable minimization. We take this upper bound to be 0.1. The derivative of the norm with respect to a i can be computed either numerically or analytically depending on the norm. After applying the similarity transformation X → hXh −1 and Y → hY h −1 , the derivative of the Hermiticity norm is where [A, B] is the standard commutator. After applying the similarity transformation to ϕ and ψ in (7.5), the derivative of the p = 2 anti-Hermiticity norm is found to be Finally, the derivative of the eigenvalue norm is computed numerically.

Results
Below we present the results obtained by applying the gauge cooling method to the Stephanov model. We will also show results for the Osborn model using the gauge cooling procedure outlined in [44]. For the runs we have used a block size N = 24, a Langevin step size ∆t = 10 −4 , and a total Langevin time t end = 1. Whenever cooling is involved, we apply 10 cooling transformations between every Langevin update. For this investigation, the Stephanov model is simulated with parameters {m = 0.2, µ = 0.5}, while the Osborn model is simulated with {m = 0.1, µ = 0.25}. The two sets of parameters were chosen in a region where the full and the phase quenched results deviate by an intermediate amount, and the two models have a comparably severe sign problem. First, we discuss the effect of cooling on the distribution of the Dirac eigenvalues. We start with results for the eigenvalue norm, see Fig. 13. On the right hand side, we present results for the Osborn model [44], which show that applying gauge cooling using the N ev norm results in the eigenvalue distribution developing a "wedge" that excludes zero. In contrast, this does not happen for the Stephanov model, as can be seen in the left figure.
In this case the distribution of the cooled CL evolution is even wider than that of the uncooled CL evolution.
Since it is surprising that cooling with the eigenvalue norm results in a wider eigenvalue distribution for the Stephanov model, we have studied the dependence on the cooling parameters ξ and the number of eigenvalues included in more detail. In Fig. 14 we show scatter plots of the Dirac eigenvalues for small ξ (left) and large ξ (right). For small values of ξ the eigenvalue distribution turns into a spherically symmetric ring, which gives rise to a vanishing chiral condensate; see Fig. 15 for the chiral condensate and the baryon number density as a function of ξ. For increasing ξ the eigenvalue distribution becomes more elongated along the imaginary axis, and at ξ ≈ 60, the central hole in the eigenvalue distribution disappears. For large ξ, see Fig. 14 right, the eigenvalue distribution approaches the spectral domain of the quenched theory, albeit with many more outlying eigenvalues. The chiral condensate and the baryon number density in Fig. 15 show a continuous dependence on ξ up to ξ ≈ 60, and take on approximately constant values beyond this point. The chi-  ral condensate approaches its quenched value, while the baryon number remains different from the quenched result. A possible interpretation of these results is that for small ξ the cooling process moves the CL trajectories to a Lefschetz thimble that does not give the correct dynamical result, while for large ξ, many Lefschetz thimbles that contribute each with their own sign, wove the result in the direction of the quenched result because CL does not take this phase factor into account.
Results for the eigenvalue distributions obtained by cooling with the anti-Hermiticity norm N AH are shown in Fig. 16. Once more we observe that the Osborn model is susceptible to the effects of gauge cooling, while it has no effect on the fermionic eigenvalues of the Stephanov model.
We also looked at the evolution of the norms N AH and N ev as a function of the Langevin time, see Figs. 17 and 18, respectively. The plots show that the evolution of the norm reflects the eigenvalue situation. Whereas for the Osborn model the norm is clearly reduced by the corresponding cooling algorithm, no such improvement is seen for the Stephanov model. Even the shifted representation, which could leverage its more advantageous initial condition (the Dirac operator is anti-Hermitian) for cooling to work, simply falls back to that of the uncooled, unshifted CL. The difference between the Osborn model and the Stephanov model is reminiscent of the difference in convergence of the CL algorithm between U(N ) and SU(N ) one-dimensional lattice QCD models [41,46,47].

Deformation technique
Another procedure which attempts to fix the issues CL has for simulating systems at finite chemical potential was proposed in [48]. The basic idea is to deform the Dirac operator such that its eigenvalues are removed from the region around the origin, and then extrapolate the deformation parameter to zero. We will deform the random matrix model by a finite temperature term which in essence is given by the two lowest Matsubara frequencies ±πT [19,49], and α can be thought of as the lowest Matsubara frequency. Following [50], we measure the physical quantities in question as a function of α, and then extrapolate α → 0. Beyond a critical value of α the eigenvalue spectrum opens up in the imaginary direction at which point chiral symmetry is restored. We can thus extrapolate from higher values in α for which there are no eigenvalues at the origin. This behavior is clearly demonstrated in Fig. 19 for (m, µ) = (0.2, 0.5), where we see a gap opening at α ≈ 1.0. This is however a fairly large range to extrapolate over, and what is more, these parameter values correspond to a different phase of the model. Since N is finite, the latter is not a fundamental problem though. The extrapolation problem unfortunately does not really improve if we choose values of (m, µ) where the sign problem is milder. In Fig. 20 we show a similar scatter plot for (m, µ) = (0.2, 0.35). As can be seen from the location of the origin with respect to the eigenvalue cloud this is a relatively mild case, as the origin is close to the edge. For a more quantitative approach we can also analyze the behavior of the force norm as suggested by [51]. It is postulated that if P (|F |), which is the density of the norm of the Langevin force, falls off at an exponential rate (or faster), the Langevin algorithm will give the correct result. However, if it falls off as a power law (or slower), we do not expect Langevin to converge to the right answer [51]. Therefore we define α c , from which one may extrapolate, as the first α for which the Langevin force decays as a power law or slower. This is plotted in Fig. 21 and demonstrates that the value of α c does not change much as we move from a hard to a mild problem, depending on the value of µ; we also saw this in Figs. 19 and 20 which demonstrates that the gap does not open until α ≈ 1.0. In Fig. 22 we plot the analytic solution for the random matrix theory (8.1) in the thermodynamic limit, for masses m = {0.1, 0.2} [49]. We also plot the corresponding CL results. As predicted by the histogram study of the previous paragraph we see agreement with the analytic curve for α 1.0. There are however two more crucial observations to be made. First, looking at the m = 0.2 data, we can conclude that although theoretically possible, it is infeasible in practice to extrapolate the values for the condensate and the baryon number to α = 0 due to the rapid change of these quantities in the region α ∈ [0, 1]. Second, looking at the m = 0.1 data, we observe that there is a phase transition separating the α = 0 and α ≥ 1 region, meaning that the method has a limited range of convergence in mass. This means that even if the issue of precision and statistics can be overcome to solve the first issue, there is only a limited mass range this can work for.

Conclusions
In this work we have analyzed the complex Langevin algorithm for the Stephanov model, which is a random matrix theory model for QCD at finite baryon density. This model possesses a rich structure due to a phase transition that takes place at a finite critical value of the baryon chemical potential, and separates two distinct phases with zero and nonzero baryon density, respectively. The main issue that was discussed in this paper is the convergence of the complex Langevin algorithm which is particularly problematic in the cold and dense regime as one is approaching the chiral limit. We observed that a naive implementation of the complex Langevin algorithm yields phase quenched results for the chiral condensate and the baryon number density, whose analytical expressions were also derived in this article. The issues of the wrong convergence was addressed by implementing several methods that have been suggested in the literature to rectify the pathologies of the complex Langevin algorithm, a complex Langevin reweighting method, several cooling methods and an extrapolation method. In order to shed more light on the properties of some of these methods, we also performed a direct comparison with a relatively similar matrix model, the Osborn model, which however has a milder sign problem.
We were able to recover the correct solution with a novel reweighting technique that uses complex Langevin trajectories chosen from the parameter regime where the algorithm converges to the correct solution. A striking result of the reweighting procedure was that by choosing an auxiliary ensemble at one side of the transition we could reproduce the correct results on the other side of the transition too. However, the cost of this method is exponential in the volume and therefore it does not solve the sign problem. Second, we tested the gauge cooling method, where one utilizes the enhanced gauge symmetry of the complexified action to modify the complex trajectories with the hope of retrieving the correct solution. While this method works remarkably well for the Osborn model, as was already mentioned in the literature, it fails for the Stephanov model. We carefully studied the effect of cooling using different norms on the Dirac spectrum, the chiral condensate and the baryon number density. Two of the norms yield the phase quenched results for these observables, while the so called eigenvalue norm which depends on two parameters, gives results that depend on these parameters and do not agree simultaneously with the correct analytical result for any value of the parameters. Another attempt to assist cooling was done by shifting the entire µ-dependence from the fermion determinant to the "gauge" part of the action and that was unsuccessful as well. The hope was that a different complexification, which actually starts off with the correct value of the "anti-Hermiticity norm", could potentially alleviate the convergence problem, but to no avail. Finally, we tested the so called deformation technique, which is a novel idea that was introduced only very recently, and which has produced some promising first results for QCD in small volumes. The idea is rather intriguing because the deformation parameter can be interpreted as an imaginary chemical potential or a finite temperature in the matrix model language, and it is well established by now that complex Langevin has far less problems at high temperatures due to the much milder singular drift term problem. However, it extrapolates from a parameter domain where the quark mass is outside the spectral domain of the Dirac operator, to a parameter domain where the quark mass is inside this domain, and it is not surprising that the extrapolation to zero deformation parameter cannot be made in a controlled and reliable way.
In conclusion, we have shown that the complex Langevin algorithm including cooling and deformation techniques cannot solve the sign problem of the random matrix model originally proposed by Stephanov in the domain where the quark mass is inside the spectrum of the Dirac operator. The only method that gives the correct solution is a complex Langevin reweighting method, but since the cost of this method remains exponential in the volume we cannot claim that this method solves the sign problem. What distinguishes random matrix theory from QCD is that it is a much stronger coupled theory, making it much harder for the drift term, to evolve to the correct Langevin trajectory if it indeed exists. A plausible explanation of our results is that there is no such trajectory, and that correct results can only be obtained by taking into account multiple thimbles, each with its own phase. These results do not necessarily generalize to QCD which is a much weaker coupled theory, but do raise serious concerns that the complex Langevin method does not work when the quark mass is inside the spectral domain of the Dirac operator.

A The Phase Quenched Partition Function
In this appendix we evaluate the mean field result for the free energy of the phase quenched random matrix partition function (2.1) with Dirac operator (2.2). Since in random matrix theory the number of flavors only enters as O(N f /N ) we evaluate the large N limit of this partition function for the simplest case which is N f = 2. Some of the results on this appendix also appeared in [18,52] The phase quenched two-flavor partition function (2.1) can be rewritten identically as [18] Z pq (m, µ) = dD e −N (|a−m| 2 +|b−m * | 2 +|c| 2 +|d| 2 ) det n D (A.1) The integration dD is over the real and imaginary parts of a, b, c and d. For large n, this partition function can be evaluated by a saddle point approximation. The determinant can be evaluated as det D = (|a| 2 − z 2 )(|b| 2 − z * 2 ) + a * b * c * d + abcd * + |z| 2 (|c| 2 + |d| 2 ) + |c| 2 |d| 2 . F 5 ≡ c * det D − abd * − c * |d| 2 − c * |µ| 2 = 0, (A.8) The solution of these equations occur in two different phases, the normal phase with c = d = 0 and the pion condensation phase with c = 0 and d = 0. In the first case, the equations decouple and are easier to solve.
Note that we did not use complex conjugation to prove the first relation. The equations F 5 = 0, F 6 = 0, F 7 = 0, F 8 = 0 are not linearly independent, Combining this with a * F 1 − aF 2 = 0 and b * F 3 − bF 4 = 0 we obtain a = a * , and b = b * . (A.14) As mentioned above, this does not imply that a and b are real -in fact, they are not as we will see below.

A.1 Solution for the Condensed Phase
In the condensed phase we have c and d nonzero.

A.2 Normal Phase
In the normal phase we have that c = d = 0 and saddle point equations for a, a * and b, b * decouple. We still have a = a * and b = b * with a and b given by the solutions of (a − m)(a 2 − µ 2 ) = a, For m > 0 and µ > 0 these equations have three real solutions and the correct solution is given by the one that minimizes the free energy at the saddle-point. Below the shaded area in Fig. 1