Complex Langevin Dynamics in Large $N$ Unitary Matrix Models

Using complex Langevin dynamics we examine the phase structure of complex unitary matrix models and compare the numerical results with analytic results found at large $N$. The actions we consider are manifestly complex, and thus the dominant contribution to the path integral comes from the space of complexified gauge field configuration. For this reason, the eigenvalues of unitary matrix lie off the unit circle and venture out in the complex plane. One example of a complex unitary matrix model, with Polyakov line as the unitary matrix, is an effective description of a QCD at finite density and temperature with $N$ number of colors and $N_f$ number of quark flavors defined on the manifold $S^1 \times S^3$. A distinct feature of this model, the occurrence of a series of Gross-Witten-Wadia transitions, as a function of the quark chemical potential, is reproduced using complex Langevin simulations. We simulate several other observables including Polyakov lines and quark number density, for large $N$ and $N_f$ and found excellent match with the analytic results.

Using complex Langevin dynamics we examine the phase structure of complex unitary matrix models and compare the numerical results with analytic results found at large N . The actions we consider are manifestly complex, and thus the dominant contribution to the path integral comes from the space of complexified gauge field configuration. For this reason, the eigenvalues of unitary matrix lie off the unit circle and venture out in the complex plane. One example of a complex unitary matrix model, with Polyakov line as the unitary matrix, is an effective description of a QCD at finite density and temperature with N number of colors and N f number of quark flavors defined on the manifold S 1 × S 3 . A distinct feature of this model, the occurrence of a series of Gross-Witten-Wadia transitions, as a function of the quark chemical potential, is reproduced using complex Langevin simulations. We simulate several other observables including Polyakov lines and quark number density, for large N and N f and found excellent match with the analytic results.

I. INTRODUCTION
A nonperturbative study of the phase structure of QCD at finite temperature and nonzero baryon chemical potential still remains an outstanding problem [1,2]. This is due to the fact that the fermion determinant becomes complex and the theory has a sign problem. The standard methods to study the theory, lattice QCD algorithms based on importance sampling, fail to produce reliable simulations. There have been recent developments in tackling this problem. One method is the use of complex Langevin dynamics with stochastic quantization [3,4]. This method is not based on importance sampling but instead on a stochastic exploration of an enlarged (complexified) field configuration space. Another recently proposed method is the Lefschetz thimble method [5][6][7][8][9][10], which is also based on complexification of the original real field variables.
The complex Langevin method was proposed in the early 1980s by Klauder [3,11,12] and Parisi [4]. Though it became popular in the beginning certain problems were found immediately after. First one was the problem of runaways, where the simulations would not converge and the second one was the problem of convergence to a wrong limit. In recent years the complex Langevin method has been revived, with sometimes cases of impressive success [13][14][15][16][17][18]. It has been shown recently that complex Langevin simulations produce seemingly correct answer, even when the fermion sign problem is severe, for one-, three-and four-dimensional field theories with nonzero chemical potential [19][20][21][22]. There have also been studies of supersymmetric matrix models based on complex Langevin dynamics. See Refs. [23][24][25].
In this paper, we consider a large N unitary matrix model at low temperature with a finite quark chemical potential and quark mass. This model is obtained from the one-loop formulation of QCD on S 1 × S 3 at finite temperature with finite quark chemical potential µ, quark mass m, and with N number of colors and N f number of quark flavors. After integrating out the quark and gauge degrees of freedom we obtain the model of our interest -a conventional unitary matrix model with a complex action. The unitary matrix U in this model is the holonomy (Wilson loop) of the gauge field around the thermal time circle in Euclidean space. We can use the expectation value of the trace of Polyakov line in the fundamental representation as order parameter for the phase transitions. It is zero in the confined phase and non-zero in the deconfined phase. The model is interesting as it exhibits a rich thermal phase structure. When the chemical potential passes one of the quark energy levels there is a third order Gross-Witten-Wadia (GWW) transition from a confined to a deconfined phase and back again. This model also exhibits another interesting feature known as the Silver Blaze behavior. When the quark mass is nonvanishing the bulk observables of the model are nearly zero until the onset transition to the deconfined phase, which occurs when the chemical potential reaches the value of the lightest quark mass.
In the matrix model with complex action, the dominant contributions to the functional integral come from complexified gauge field configurations. Due to this reason, the saddle point eigenvalues of the unitary matrix U lie off the unit circle, on a contour in the complex plane. The eigenvalues of U can be written as exp(iθ i ) with θ i the angle variables and i = 1, · · · , N . We can make a change of variables such that the functional integral reduces to an integral over {θ i }. At large N , the functional integral is dominated by a single saddle point but since the action is complex this saddle point configuration lies out in the complex plane where the θ i are no longer real. As a consequence, the Polyakov line and the inverse Polyakov line are not equal, that is, P = P −1 . Through complex Langevin simulations we indeed confirm this behavior. In fact the behavior of inverse Polyakov line precedes that of the Polyakov line as a function of chemical potential. This feature was observed analytically in an earlier work by Hands et al. in Ref. [26].
In this paper, we examine this large N unitary matrix model using complex Langevin simulations. It is possible to generate representative field configurations by integrating a stochastic differential equation, known as the complex Langevin equation. The drift terms arising from the complex action force the field variables to evolve in an extended (complexified) field space, in which the large regions where the observables are plagued by phase fluctuations are avoided [17].
When N is large, we can consider the gauge field, corresponding to the angles of the Polyakov line, as a distribution on a contour. From the equation of motion, the saddle point distribution of the Polyakov line eigenvalues can be calculated analytically and plotted by mapping the angles from an arc on the unit circle to a contour over the same range of angles in the complex plane [26]. The theory is said to be in a confined phase when the contour on which the Polyakov line eigenvalues are distributed is closed. The contour opens up in between quark energy level transitions giving rise to a deconfined phase in the theory. The third derivative of the grand potential is discontinuous at each energy level crossing. These are characteristic features of a third order, GWW transition [27][28][29].
This paper is organized as follows. In Sec. II we give a brief outline of the complex Langevin dynamics and stochastic quantization. In Sec. III we discuss a simple yet nontrivial matrix model called the ab-Model, which is a complexified version of the Gross-Witten-Wadia (GWW) model. This model has two phases, confined and deconfined, and it exhibits a third-order phase transition. In Sec. IV we discuss another interesting large N unitary matrix model, which arises in the one-loop formulation of QCD on compact spaces. This model possess a tower of quark energy levels due to compactification and is defined for positive and negative chemical potential values. We then focus on to a truncated cousin of this model -a single quark energy level matrix model with positive chemical potential. This model also has a complex action and captures the physics we are interested in without loss of generality. We can define a transition parameter (which is function of the temperature and chemical potential) in this model and as we change this parameter, the model exhibits confinement/deconfinement phase transitions. We show the eigenvalue distributions corresponding to the confined (closed) and deconfined (gapped) phases of the theory using complex Langevin simulations. We also simulate the behaviors of Polyakov lines and fermion number density as a function of the transition parameter. We simulate the model for a range of temperatures and chemical potentials to study its phase structure. We also show the phase diagram of the model, at low temperature, on the (µ, β) plane, in the vicinity where quark energy level equals the chemical potential. We then simulate the model at large quark mass and show that the bulk observables exhibit the Silver Blaze behavior -the observables are roughly zero until the onset transition to the deconfined phase, which occurs when the chemical potential equals quark mass. We then move on to discuss the single-level model with a simple nontrivial gauge interaction turned on. We study the behavior of observables as a function of the interaction parameter. We see that the model prefers to stay in the confined phase as the interaction strength is increased. In Sec. V we provide conclusions and discussions. In Appendix. A we use complex Langevin dynamics to simulate QCD on S 1 × S 3 at finite chemical potential and low temperature. We are able to reproduce the series of GWW transitions, as a function of the chemical potential, as described in Ref. [26]. Our simulations also reproduce the level structure feature of the bulk observables -fermion number density, pressure and energy -of the model. In Appendix B we investigate the reliability of complex Langevin method by studying the probability distribution for the magnitude of the drift term and the Langevin runtime history of the unitarity norm. We note that the probability distribution for the magnitude of the drift term falls of (possibly) with a power law even though the simulations show excellent agreement with analytical results. We think that these diagnostics need further investigations and we save it for future work.

II. COMPLEX LANGEVIN DYNAMICS
The central idea of stochastic quantization is that expectation values of observables are obtained as equilibrium values of a stochastic process [30,31]. In order to achieve this we evolve the system in a fictitious time τ , subject to a stochastic noise. That is, the system evolves according to Langevin dynamics. When the action is complex it is still possible to consider Langevin dynamics. The force (gradient of the action) becomes complex in this case making the fields also complex during the evolution.
In this work we make use of complex Langevin dynamics with stochastic quantization to study large N unitary matrix models with complex actions. They exhibit sign problem due to the fact that the action is complex. Standard Monte Carlo methods fail to produce the correct equilibrium distributions of these models. We can use discretized complex Langevin equation with Euler method (which is a first order algorithm) to find the equilibrium field distributions of these models.
We note that in unitary models with real action the domain of the angular variables θ i , with i = 1, · · · , N , is [0, 2π). After complexification the domain becomes a strip with the the domain [0, 2π) along the real directions and (−∞, ∞) along the imaginary directions. The range of e iθi , that is, the complexified eigenvalues of U has the whole complex plane as the range. Let us take θ i (τ ) as the complexified angle variables of the gauge link U (τ ) at a Langevin time τ . (From now on we take θ i to be complex, in this paper, unless otherwise specified.) We have the discrete Langevin evolution equation where ∆τ is the Langevin time step, and η i (τ ) is a Gaussian random variable satisfying the conditions If the action S is of the order N 2 , then strictly at infinite N the fluctuation term in Eq. (1) could be safely dropped. Moreover, to reduce excursions in the imaginary directions of the field configurations, which would spoil the validity of the method, we should use real Gaussian random variables [32][33][34].
We also need to impose the SU (N ) constraint on the complexified angular variables after each Langevin time step. That is, we need This can be easily implemented by subtracting the average value θ av (τ ) from each θ i (τ ) variable, i.e.
Note that this condition is implemented in a holomorphic way. That is, both of the real and imaginary parts of θ av (τ ) are subtracted. Ideally, one should eliminate one variable (say θ 1 ) using the constraint Eq. (3) and stochastically quantize the remaining variables. To proceed we need to justify that our method of imposing the constraint after each time step leads to the same result.
A set of stochastic flow equations involving the gradient of the action like the one given in Eq. (1) is invariant under the orthogonal transformation of variablesθ where O is an orthogonal matrix. In terms of the transformed variables, the set of equations is where we have used the orthogonality of matrix O. Orthogonality also guarantees that new random variablesηs satisfy the condition Eq. (2). Now, we can always choose an O such thatθ 1 = 1 √ n i θ i . In terms of the transformed variables it is easy to understand why our method works. The constraint Eq. (3) is now rewritten simply asθ 1 = 0. If we start with a set of variables which already satisfies this constraint then a valid Langevin time evolution step may be performed by simply discarding any evolution inθ 1 . This is precisely our method of imposing constraint after each time step, rewritten in terms of the new variables. To emphasize, one can straight forwardly argue that in terms of old variables, this step is same as Eq. (4). Our argument works for any arbitrary linear constraint.
We note that there also exists another complementary method in which one could implement complex Langevin dynamics directly on the matrix variables U (τ ). In this case the evolution equation takes the form where the matrix R is a stochastic unitary matrix. We note that this method can be used for studying similar models in higher spacetime dimensions.
In this paper, we use the first method described above where the link field U is diagonalized and the SU (N ) constraint has been imposed.
We note that the complexification of the dynamical variables in the theory can change the Langevin evolution drastically. There can be unstable directions on the complexified field configuration space and the Langevin evolution can converge to wrong limits. One should be aware that the numerical integration must be performed carefully when the Langevin trajectory makes a large excursion into imaginary directions. One could, in principle, use a small step size but it still has two problems: (i) it does not solve instabilities in all directions and (ii) it will result in a slow evolution, which can be computationally very inefficient. In order to take care of both of these problems we follow the algorithm given by Aarts et al. in Ref. [35]. We consider an adaptive step size in the discretized complex Langevin equations. We compute the absolute value of the maximum drift, K max , at a given Langevin time τ and the stepsize for the next evolution step is taken to be where γ is a number chosen according to the model we want to simulate. In our simulations we typically take γ to be O(1).

III. AB-MODEL
To demonstrate the effectiveness of Complex Langevin Dynamics, we begin by studying a simple, yet nontrivial model -a complexified version of Gross-Witten-Wadia (GWW) Model [27][28][29]36]. We refer to our model as ab-Model. It has two phases, confined and deconfined, exhibiting a third-order phase transition. The action is given by where a, b ∈ C, U is an element of SU (N ), and when a = b it becomes the Gross-Witten-Wadia model. Before proceeding further let us make a few generic comments. A linear term in Tr U breaks the center symmetry. Furthermore, the above action (or other polynomial generalization of it) is complex. If a = b, then the Z 2 symmetry U → U † is broken. This implies Tr U = Tr U † . One may ask, that what it means in terms of manifestly gauge invariant operators. This means that the contribution from baryon and anti-baryon is different. Another related observation is one may naively expand Eq. (12) in a series Here we have separated the "mesonic" and "baryonic" contributions. Due to the center symmetry only a center symmetry invariant combination of Tr U and Tr U † contributes. By mesonic contribution we mean product of traces for which sum of powers all the occurrence of unitary matrix and its inverse sum to zero. For a baryonic operator, the sum is only zero up to modulo N , i.e., proportional to a non-zero integral power of N . If baryonic contributions are neglected then Eq. (12) is equivalent to a model with parameters, a = b = (ab). We will later see that for center symmetry invariant operators, this equivalence is actually held in the ungapped phase.
Expressing the action in diagonal gauge, the effective action becomes where the first term is the Vandermonde piece and M is the Lagrange multiplier which ensures that det(U ) = 1. At large N , the theory is dominated by the saddle-point equation which gives the equation of motion On substituting z i = e iθi the equation of motion becomes and M is given by In the saddle point, M may have a nonzero value and could be thought as effective baryon number.
At N → ∞ limit, we can replace the summation by an integral over a nondecreasing function and performing a change of variables from s to complex variables z(s) the equation of motion becomes and P implies we are taking the principal value of the integral.

A. Ungapped Phase
In the GWW model, it is known that for small potential, i.e., a < 0.5, the theory is in an ungapped phase. Assuming a similar picture also holds for the ab-model, we solve it by taking an ansatz for ρ(z) in ungapped phase as, Comparing with the left hand side of Eq. (22) we have Therefore ρ becomes, We also find which indicates that the theory is in an ungapped phase. Demanding normalization of ρ(z) we fix A 1 = 1. Therefore, We can solve for the contour, where ρ(z) is positive definite, by integrating Eq. (21) Since s is purely real, and assuming that z = r(θ)e iθ , a = |a|e iφ1 and b = |b|e iφ2 , the above equation is satisfied only if the real part of the right hand side is zero. That is, To fix c, we invoke the condition that det(U ) = 1, i.e., where the branch-cuts are taken from z = 0 to the point z(±π). Replacing ln(z) using Eq. (30), the above equation becomesˆd Hence the contour is got by solving the transcendental equation Now we can compare the distribution of eigenvalues from complex Langevin dynamics with the analytic result for any (a, b) combination. In Fig. 1 we show the analytical result and the data obtained through complex Langevin simulations without noise for parameters a = 0.35, b = 0.2 and N = 100. In Fig. 2 we show the result with Gaussian noise turned on. We see an excellent agreement between the analytical and numerical results. The data are obtained through complex Langevin simulations without noise. We used a fixed Langevin step size ∆τ = 0.00001 and evolved the system for 45000 steps. The dashed unit circle is guide to the eye. We also note that the complex Langevin simulations show excellent agreement with analytical results when the parameters are also complex. In Fig. 3 we show the analytical result and the data obtained through complex Langevin simulations without noise for parameters a = 0.2 + i0.2, b = −0.1 + i0.1 and N = 100. In Fig. 4 we show the result with Gaussian noise turned on. The data are obtained through complex Langevin simulations without noise. We used a fixed Langevin step size ∆τ = 0.00001 and evolved the system for 45000 steps. The dashed unit circle is guide to the eye.

B. Gapped Phase
In the gapped phase, similar to GWW model, the eigenvalues lie on an open contour C.
To study this phase, we employ resolvent/spectral-curve method used in Ref. [26], and reviewed in Ref. [37]. The resolvent is defined as At large N limit, ω(z) is analytic everywhere in the complex plane, except along a square-root branch cut running along C, and expressed as For a given potential V (z), the equation of motion (similar to Eq. (22)) can be expressed in terms of ω(z) using the Plemelj formulae where z ± lies on either side of the branch cut and → 0 limit is taken. We can also express ρ(z) as the discontinuity of ω(z) across the cut C as The expectation value of any function G(z) can be found aŝ For ab-model wherez,z * are the end points of branch cut C and f (z) is an unknown function, which remains to be fixed. Since ω(z) has to be regular over the entire plane except along C and the origin we can fix the form of f (z) as Therefore ω(z) becomes (substitutingz = Re iφ ) Normalization of ρ(z), from Eq. (37), translates to and lim |z|→∞ ω(z) = −1.
This fixes f (z) as We also get two more relations between R, M and cos(φ) and To fix the three unknowns completely, we need a third equation, which comes from invoking the det(U ) = 1 condition, from Eq.
whereC is a contour encircling the branch cut C, and the branch cut of ln(z) ranges from (−∞, 0). Deforming the contour Fig. 5 to the one in Fig. 6 and evaluating in → 0 and Γ → ∞ limits, we find that the divergences arising from the cutoffs Γ and cancel separately and we arrive at the following condition From Eq. (27), we can numerically compute M, both in ungapped and gapped phases, and compare it against analytical results. Choosing b = 2.0a and varying a from 0 to 1.2, we find that it matches very well both in ungapped and gapped regimes -see Fig. 7. (Gap opening point can be found from Fig. 10.) Similarly we compare other observables, Tr (U ) and Tr (U −1 ) . Analytically Tr (U ) is given by, and Tr (U −1 ) is given by In Fig. 8 and Fig. 9 we show the observables Tr (U ) and Tr (U −1 ) , respectively. We see that the analytical and numerical results show excellent agreement.

C. Phase transition of ab-Model
The eigenvalue density Eq. (29) on contour Eq. (35), is proportional to ds, which in terms of r(θ) is given by which is not positive definite for all (a, b) combinations. It fails to do so, when the function inside the brackets, [. . . ], becomes negative. Restricting to a, b ∈ R, the condition simplifies as the gap opens about θ = 0 From Eq. (35) r(0) is given by The phase diagram of the model is shown in Fig. 10. It would be interesting to know how quantities change across the gap opening transition and also the order of the phase transition. To study that we first restrict ourselves to a special case, b = 0 in our model. Then from Eqs. (56) and (57) the gap opens about a = 1 e , R = e, and since the ungapped phase has no branch cuts in the eigenvalue distributions, φ should start from zero, about the gap-opening point. And the conditions Eqs. (48) and (49) The observable Tr (U ) becomes Since the first derivative of free-energy F [a] is the expectation value of Tr (U ) we find that it is continuous across the gap. Upon expanding about a = 1 e + δa, cos(φ) = 1 − 2δp and R = e + δR (62) the variation of δ Tr (U ) is given by From Eqs. (58) and (59) we get and δp ln(δp) = δR e . (65) Eliminating δR from above two equations we get the equation To invert the above equation let us substitute δp → e k . Then we have The above equation is of the form, xe x = y, which can be inverted to express x as a function of y and it is known as the Lambert-W function [38]. (It is often expressed as W c (y).) This function is in general a multivalued-complex function, where c ∈ Z, chooses each branch. Since δa > 0 and δp ∈ R we have two real valued branches: W 0 (y) (the principal branch) and W −1 (y). Therefore, For small values of δa we know that Now the second derivative of free energy goes to zero as δa → 0 and is continuous across the gap. However, the third derivative diverges as δa → 0. Hence it has a third order phase transition. It can also be shown that similar arguments hold in the generic case b = 0. Thus we conclude that the ab-model displays a third order phase transition.

IV. GAUGE THEORY TO UNITARY MATRIX MODEL
A unitary matrix model arises in a one-loop formulation of QCD [and analogous SU (N ) gauge theories] on compact spaces (often S 1 × S 3 ). This was originally derived in Refs. [39][40][41][42] for theories with more general matter content.
The one-loop effective action of QCD on S 1 × S 3 with inverse temperature β, chemical potential µ and quark mass m has the following form [26], with thermal Polyakov line as the unitary matrix model where R is the radius of S 3 and N f is the number of flavors of fundamental fermions. The quadratic term in Polyakov loop is the contribution from adjoint fields and the linear term is the contribution from the fundamental matter fields. Here, we have taken the adjoint contribution to be bosonic and the the contribution from fundamental fields to be fermionic.
To be noted is that in the free theory the effective action is determined in terms of single particle (bosonic and fermionic) partition functions and Also note that we will be using dimensionless variables β/R, µR and mR in numerical simulations. An analogous action, for the simpler 0 + 1 dimensional case would be, and where the parameter m is the mass of the fundamental fermions.
In the low temperature limit, β → ∞, we have z b (∞) = 0 and so the gluonic contribution is negligible. Thus the action is where S Vdm is the Vandermonde piece of the action and S f is the fundamental fermionic contribution. The fermionic part could be summed in a logarithm where We would like to simulate the action given in Eq. (79) using complex Langevin method. We can study several interesting observables in this model. We briefly describe them below 1. Polyakov line P and inverse Polyakov line P −1 These are the most natural set of observables to study the confined/deconfined phases in the theory.
3. Fermion number f N It gives the number of fermions minus the number of anti-fermions in a given volume In the model we study here we have a single chemical potential µ. In general there can be chemical potential for each fermion flavor.
The quark number susceptibility χ f measures the response of the fermion number density to infinitesimal changes in the chemical potential, This observable follows the behavior of the Polyakov line. Thus, it also serves as an indicator of confinementdeconfinement transitions for nonzero chemical potential.

Pressure
with V 3 denoting the spatial volume.

Energy E
It can be constructed from pressure and fermion number density It is also possible to compute the chiral condensate and average phase, though we will not compute them in this work. The chiral condensate ψψ is given by and the average phase e iφ pq has the form where pq refers to the phase quenched theory.

B. Single Level Model with Positive Chemical Potential
We can truncate the action given in Eq. (79) in a double scaling limit: where 0 is a fixed quark energy level and we call ξ the transition parameter.
Only contribution from a single level survives here and the action takes the form The effective action on the complexified angle variables include the Vandermonde piece and a Lagrange multiplier.
In the large N limit, the integral over the angles is dominated by a saddle point obtained by solving the equation of motion that follows from the effective action involving Eq. (89) Here also the action is not hermitian, giving rise to the sign problem in the presence of a chemical potential. As a result the saddle point configuration will lie out in the complex plane. If we define z i = exp(iθ i ) then in the presence of the non-real potential the z i will move off the unit circle in the z-plane.
We can explore the nature of eigenvalue distribution in the complex plane for various values of transition parameter ξ. We find that when ξ is either very small or large, the potential vanishes and so we expect the {z i } to be uniformly distributed around the unit circle. Thus, when µ varies from µ to µ the quark energy level becomes occupied and the effective fermion umber jumps by factor σ. In Ref. [26] the authors provide a detailed description of this transition.
Let us look at the various regimes of ξ and see how it affects the eigenvalue distribution, following the analytical study given in Ref. [26].

The small ξ confined phase
In the small ξ confining phase the effective fermion number vanishes, N = 0, and the Polyakov line expectation values are Thus we have P = P −1 , as a result of the complex action.
As ξ is increased the contour of eigenvalue distribution opens into an arc, just as the matrix model solved by Gross and Witten [27] and Wadia [28,29].
The line of phase transitions in the (µ, T ) plane corresponds to the straight line Note that is approximation is valid only in the low temperature (β → ∞) limit.

The large ξ confined phase
In this phase the effective fermion number is indicating that the level is now occupied.
The Polyakov line expectation values are Comparing with the previous case the behavior of P and P −1 swaps over along the replacement ξ → ξ −1 .
The large ξ confined phase persists until the value For smaller values of ξ the contour of eigenvalue distribution is not closed and the phase does not exist. The points of transition ξ = ξ 1 and ξ = ξ 2 satisfy ξ 1 ξ 2 = 1.
In the (µ, T ) plane the boundary lies along the straight line again valid in the low temperature limit.

The deconfined phase
In the region ξ 1 ≤ ξ ≤ ξ 2 , experience with GWW matrix model suggests that the eigenvalue distribution exhibits the shape of an open contour.
In this regime we get a condition This equation determines N as a function of ξ.
From the above equation it follows that across the transitions at ξ = ξ 1 and ξ = ξ 2 , fermion number density N and its first derivative ∂N /∂µ are continuous, however higher derivatives are discontinuous. Since N is the effective fermion number, the first derivative of the grand potential, it follows that the transitions are third order, just as in the original GWW model.
For a single winding, the Polyakov lines are Using complex Langevin dynamics we have simulated the single level matrix model given by the action in Eq. (89). In Fig. 11 we show the eigenvalue distributions of the Polyakov line in the confined and deconfined phases as a function of the logarithm of the transition parameter, log ξ, for SU (N ) case with N = N f = 500 and quark mass m = 0. We see that the eigenvalue distributions start with a closed contour (confined phase), passes through an open contour (deconfined phase) and again goes into a closed contour. (This figure can be compared with Fig. 12 in Sec. 4.1 of Ref. [26], where it was obtained through analytical methods.) In Fig. 12 we provide the (normalized) effective fermion number f N , and in Fig. 13 the Polyakov line expectation value P and the inverse Polyakov line expectation value P −1 across a pair of GWW transitions from the small ξ confined phase through the deconfined phase to the large ξ confined phase. The transitions from confined/deconfined phases occur when either P or P −1 vanish. The parameters used are: N = N f = 3 and 500 and quark mass m = 0. The simulations show excellent agreement with the analytical results in the large N .
In Figs. 14 and 16 we show the Polyakov lines and fermion number density for a range of simulation parameters of the single level matrix model (see Eq. (89) for the form of the action): β = {10, 15, · · · , 100} and µ = {3.0, 3.025, 3.05, · · · , 4.0}. The quark energy level of the model is fixed to the third level ≡ (l=3) = 3.5. The Polyakov loops peak around µ = 3.5 in this model. In Fig. 14 we show the behavior of Polyakov and inverse Polyakov loops for β = {25, 50, 75, 100}. It is clear that the widths of the Polyakov loops decrease as the temperature is reduced (large β) and the behavior of inverse Polyakov line precedes that of the Polyakov line as a function of µ. In Fig.  15 we show the Langevin evolution history of the Polyakov loop observable in this model for β = 50, 75 and with µ = 3.0, 3.3, 3.5 for each β value. We note that the observables saturate to their equilibrium values rather quickly in this model. In Fig. 16 we show the behavior of the (normalized) fermion number density f N as a function of chemical potential and inverse temperature. The transition in fermion number becomes sharper as the temperature is decreased (high β). The model is in a deconfined phase when 0 < f N < 1.
When the quark mass is non-vanishing in QCD, the expectation values of bulk observables such as the fermion number density, Polyakov lines and energy, exhibit the 'Silver Blaze' behavior. The bulk observables are nearly zero until onset [43] to a deconfinement transition, which occurs when the chemical potential increases to the value of the lightest quark mass. We simulate the model given by the action in Eq. (79) to see this phenomenon. In this model the onset occurs at µ = m. The Polyakov line is given in Fig. 17 (Left) as a function of chemical potential for large quark mass, near the onset µ = m = 25 for N = N f = 500 and β = 25 (low T ). In the large m limit, similar to the m = 0 case, the behavior of inverse Polyakov line P −1 precedes that of P as a function of µ. The transition in µ occurs around onset at m. In Fig. 17 (Right) we show the effective fermion number density as a function of chemical potential. As we can see in the figures the bulk observables are close to zero until the onset transition at µ = m. The observables rise smoothly from the onset and as µ is increased further from m the observables behave as they would for m = 0. This is reflected in the oscillations that appear in the observables at larger µ. The oscillations in the Polyakov and inverse Polyakov loops are clearly visible. In order to see the prominent nature of oscillations in the fermion number density one has to normalize this observable by its Stefan-Boltzmann value. (See Ref. [26] for a discussion on this.)

C. Single Level Model with U and U †
In this section we consider the phase diagram of the model given by the following action where Such a model naturally arises from 0 + 1-dimensional gauge theory with a fundamental fermion. In Fig. 18 we provide the phase diagram of this model on the (µ, β) plane for the level l = 1. (Corresponding to quark energy level = 1.5 and σ = 4.) From the behavior of the expectation value of the fermion number density we see that the phase transition from confined to deconfined phase is smooth on the (µ, β) plane even at high temperature (0.1 ≤ β ≤ 2.0).

D. Single Level Model with Interaction
It would be interesting to consider the single-level matrix model with a nontrivial interaction turned on. We take a Polyakov line interaction term of the form Here g denotes a coupling parameter. Thus we have  Here N = N f = 500 and β = 25 (low T ). The data are obtained through complex Langevin simulations with an adaptive Langevin step sizes ∆τ ≤ 0.00005, thermalization steps N therm = 5000, generation steps Ngen = 5000 and with measurements performed with an interval of 50 steps. Here also we take the quark energy level to be fixed at ≡ (l=3) = 3.5. The action is again not hermitian, giving rise to the sign problem in the presence of a chemical potential. In Figs. 19 and 20 we plot the fermion number density and the Polyakov lines of the interacting model for various values of the coupling g = 0, 5, 20, 100. It is evident that the confinement/deconfinement transition becomes sharper as the interaction strength is increased. The behavior of the Polyakov lines show that the model is in a confined phase for most of the values of the chemical potential.

V. CONCLUSIONS AND DISCUSSIONS
In this work we have successfully used complex Langevin dynamics with stochastic quantization to simulate the thermodynamics of large N unitary matrix models with complex actions. We started with a simple matrix model called the ab-model and investigated its phase structure analytically and numerically. The numerical simulations show excellent match with analytical results. We also studied a model obtained from the effective theory of QCD on S 1 ×S 3 at low temperature and finite quark chemical potential. At zero quark mass and low temperature our simulations showed a series of GWW confinement-deconfinement phase transitions as a function of the chemical potential. The phases are characterized by the distribution of eigenvalues of the Polyakov line in the complex plane. In the large quark mass regime we were also able to observe the Silver Blaze behavior in that the bulk observables are roughly zero until the onset transition to the deconfined phase, which occurs at µ = m. We also simulated the model with a simple nontrivial Polyakov loop interaction turned on. The model prefers to live in the confined phase as the interaction strength is increased.
We also note that each confinement-deconfinement transition in the Polyakov loop is associated with a quark energy level transition. It is interesting to note that the non-monotonic behavior of Polyakov loops have been observed in lattice simulations of QCD with gauge group SU (2) near its saturation density in Ref. [44].
We successfully applied complex Langevin dynamics to QCD on S 1 ×S 3 with finite chemical potential and computed several bulk observables. We provided our simulation results on this in Appendix A.
There are several interesting future directions. One could consider complex Langevin simulations of the model with several quark flavors with masses m f and different chemical potentials µ f . One could also add other types of nontrivial interaction terms into the model and look for cross-over transitions on the (µ, β) plane [45]. It would also be interesting to see if there exists an AdS/CFT type gravitational dual of the models we studied here. One could ask the question whether the infinite sequence of GWW transitions that we observe in the matrix model can be seen in the dual gravitational description.  The data are obtained through complex Langevin simulations with adaptive Langevin step sizes ∆τ ≤ 0.000005, thermalization steps N therm = 5000, generation steps Ngen = 5000 and with measurements performed after every 50 steps. The solid lines are guid to the eye. The plots indicate that the model prefers to stay in a confined phase as the interaction strength g is increased.

Fermion number fN
In Fig. 22 we show f N as a function of µ for low temperatures for m = 0. The presence of an occupation level structure is evident. The transitions occur when l − µ changes sign, that is, when µ passes a quark energy level.
It is interesting to compare with the results obtained in Ref. [26]. We also note that in Ref. [46] Banerjee and Chandrasekharan observed the same level structure in the particle number in the nonlinear O(2) sigma model.
The fermion number can be used as an order parameter of the confinement-deconfinement transitions in the large N theory. The first and second derivatives of the grand potential, f N and ∂f N /∂µ are continuous as a function of the chemical potential but the third derivative ∂ 2 f N /∂µ 2 is discontinuous. This indicates that the transitions are third order, of the GWW type.

Polyakov Lines P and P −1
When the chemical potential is zero the Polyakov line P and the conjugate Polyakov line P −1 coincide and it is no longer the case for non-zero chemical potential. In Fig. 22 we show P and P −1 as a function of µ. Each spike in P and P −1 corresponds to a level transition in f N . They exhibit similar behavior as a function of µ however, the the behavior of P −1 always precedes that of P at the start and finish of each level transition. We note that the lines peak at µ = 1.5, 2.5, · · · . We also note that the widths of deconfined regions increase as µ is increased.

Pressure p and Energy E
In Figs. 23 and 24 we provide the pressure multiplied by the 4-volume and energy E = − p + µ f N . We note that the pressure exhibits a level structure. The energy levels are not horizontal. The factor µ in front of the fermion number causes the levels to rise linearly with µ.
In Fig. 25 we show the eigenvalue distributions in the confined and deconfined phases as a function of the quark chemical for N = N f = 30 and barious µ values.

Appendix B: Reliability of Complex Langevin Dynamics
We would like to justify the use of complex Langevin dynamics for the matrix models we simulated in this work. In Ref. [47,48] the authors suggested a possible criterion to determine the correct convergence of the complex Langevin method -the probability distribution of the magnitude of the drift term should fall off exponentially or faster. This criterion can, in general, be violated if the complexified fields develop large imaginary parts (the excursion problem). In Fig. 27 we show the probability distributions P (u) for the magnitude of the drift term of the single level SU (N ) matrix model. However, in our case the plots hint that the probability distribution falls off like a power law with u even though we have excellent agreements with analytical results. Figs. 12 and 13 show excellent agreement between simulation and analytical data in this model. We also observed a similar fall off behavior in the ab-model. We think this needs further investigations and we save it for future work.
It is desirable to have a well localized distribution of dynamical variables of the theory in the complexified field configuration space. A convenient measure of the size of the distribution in imaginary directions of the field variables is the unitarity norm [49] defined as of Langevin time for the single level SU (N ) matrix model with N = N f = 500, quark mass m = 0 and inverse temperature β = 30 (low T ). We see that the unitarity norm remains bounded in the simulations. In Fig. 28 we show the Langevin evolution of the Polyakov line observable for the same set of parameters.