Complex Langevin Method on Rotating Matrix Quantum Mechanics at Thermal Equilibrium

Rotating systems in thermal equilibrium are ubiquitous in our world. In the context of high energy physics, rotations would affect the phase structure of QCD. However, the standard Monte-Carlo methods in rotating systems are problematic because the chemical potentials for the angular momenta (angular velocities) cause sign problems even for bosonic variables. In this article, we demonstrate that the complex Langevin method (CLM) may overcome this issue. We apply the CLM to the Yang-Mills (YM) type one-dimensional matrix model (matrix quantum mechanics) that is a large-$N$ reduction (or dimensional reduction) of the $(D+1)$-dimensional U$(N)$ pure YM theory (bosonic BFSS model). This model shows a large-$N$ phase transition at finite temperature, which is analogous to the confinement/deconfinement transition of the original YM theory, and our CLM predicts that the transition temperature decreases as the angular momentum chemical potential increases. In order to verify our results, we compute several quantities via the minimum sensitivity method and find good quantitative agreements. Hence, the CLM properly works in this rotating system. We also argue that our results are qualitatively consistent with a holography and the recent studies of the imaginary angular velocity in QCD. As a byproduct, we develop an analytic approximation to treat the so-called ``small black hole"phase in the matrix model.


Introduction and Summary
Rotating objects are ubiquitous in our nature. If the objects are macroscopic, they may obey thermodynamics. However, in statistical mechanics, investigating rotating interacting many body systems at thermal equilibrium are problematic. The chemical potential for the angular momentum, which is equivalent to the angular velocity for point particles, makes the Euclidean actions complex even for bosonic systems, and the standard Monte Carlo method (MC) does not work because of the sign problem (see Sec. 3). Hence overcoming this issue is a quite important challenge in theoretical physics.
There are various approaches to the sign problem, such as the Lefschetz-thimble method [1,2] and the complex Langevin method (CLM) [3,4]. The Lefschetz-thimble method handles the sign problem by deforming the integration contour to mitigate the sign problem. The CLM, which we focus on in this paper, is a stochastic process for the complexified variables. The advantage of the CLM is that it allows us to study large systems. At first, the CLM has suffered from the problem that CLM results in a wrong result without noticing it. Recent studies [5][6][7][8][9] have clarified the conditions that the CLM results converge to the correct result equivalent to the path integral. The sufficient condition to obtain the correct result is that the probability distribution of the drift norm falls off exponentially or faster (see Sec. 5). The drift norm is a byproduct of solving the Langevin equation numerically, which can be calculated with negligible extra CPU cost.
To test whether the CLM works in rotating systems, we study the U(N) matrix quantum mechanics whose Euclidean action is given by [10,11], Here X I (t) (I = 1, · · · , D) are N × N hermitian matrices. D t := ∂ t − i[A t , ] is a covariant derivative and A t (t) is the gauge field of the U(N). g is a coupling constant, and we take the 't Hooft limit N → ∞ and g → 0 with a fixed 't Hooft coupling λ := g 2 N [12]. The model is invariant under the local U(N) gauge transformation X I → UX I U † and D t → UD t U † . Since we will study the system at finite temperature, we have introduced the Euclidean time t and the inverse temperature β = 1 T . Note that A t is a vector and X I are scalars in this one dimension.
This model has an SO(D) rotation symmetry X I → O I J X J , and the angular momentum is conserved. Then, we analyze this model at thermal equilibrium with a finite angular momentum chemical potential through the CLM. Although the model is strongly coupled and the standard perturbative analysis does not work, the model without rotation has been analyzed through the MC method [13][14][15][16][17][18][19][20][21][22][23][24][25] and other non-perturbative methods such as the minimum sensitivity [26][27][28][29] and the 1/D-expansion [18,19,[30][31][32][33]. Particularly, the minimum sensitivity and the 1/D-expansion would work even in the presence of the angular momentum chemical potential. Thus, we can test the CLM by comparing these methods. This is one advantage of studying this model. (Indeed, some earlier results through the 1/Dexpansion have been reported in Ref. [32].) Actually, we will show that the CLM reproduces the results of the minimum sensitivity quantitatively, when the chemical potential is not so large. (We investigate the model at D = 9 and D = 16, and we find that the minimum sensitivity is slightly better than the 1/D-expansion there.) Although the CLM does not provide reliable results at larger chemical potentials, it is enough to observe non-trivial phase transitions. Hence, the CLM partially overcomes the sign problem of the angular momentum chemical potential in the model (1.1). As far as we know, this is the first example that quantum many body systems in thermal equilibrium with a finite angular momentum chemical potential are solved through the first principle computation.
Before going to explain the details of our analysis, we present the importance of the model (1.1). This model appears in various contexts of high energy physics. Here we list some of them related to the current work: • This model is a large-N reduction [34] (or dimensional reduction [35]) of the (D + 1)dimensional U(N) pure Yang-Mills (YM) theory to one dimension. In this view, X I are the dimensional reductions of the spatial components of the original (D+1)-dimensional gauge fields. It is known that this model in the large-N limit is confined at low temperatures and shows a confinement-deconfinement transition at finite temperature [13,14,26,31]. This phase structure is similar to that of the original YM theory. Hence, this model is important as a toy model of the original YM theory. Particularly, if we turn on the angular momentum chemical potential in our model, it may be regarded as a toy model of a rotating pure gluon system [36][37][38][39][40][41], which is actively being studied to reveal natures of neutron stars and quark gluon plasma in relativistic heavy ion colliders [42].
• This model at D = 9 appears as a low energy effective theory of N supersymmetric D-particles [11] on Euclidean S 1 β × S 1 × R 8 in type IIA superstring theory [14]. 1 In this picture, a gravity dual [43,44] of the model (1.1) is given by black brane geometries. There, the aforementioned confinement/deconfinement transition corresponds to a Gregory-Laflamme (GL) transition [14,45,46].
• This model is also a toy model of N D-particles, which may be regarded as a microscopic description of a black hole [47]. 2 In this picture, the diagonal components of the matrix X I ii (i = 1, · · · , N) represent the position of the i-th D-particle on the I-th coordinate and the off-diagonal components X I ij (i = j) represent the open strings connecting the i-th and j-th D-particles, which induce interactions between the D-particles. In this way, the model (1.1) describes quantum mechanics of N interacting D-particles in the D-dimensional space. The interactions cause attractive forces between the D-particles and they compose a bound state. This bound state is chaotic and is regarded as a 1 Here S 1 β is the temporal direction of the model (1.1) and another S 1 is a Scherk-Schwarz circle whose radius is taken large. This system is described by a supersymmetric version of the model (1.1), but the Scherk-Schwarz circle breaks the supersymmetry and makes the fermions on the D-particles massive. When the radius of the Scherk-Schwarz circle is taken large, the fermion's masses become large and they can be ignored at low energy. Note that the diagonal components of X I represent the positions of the D-particles in the S 1 × R 8 space, and the scalar for the Scherk-Schwarz circle, say X 1 , is not distinguishable to other X I due to the large circle limit. Then, the low energy effective theory becomes the model (1.1). toy model of a black hole. 3 Hence, when the system rotates, it may correspond to a rotating black hole.
In this article, we mainly focus on the last D-particle picture, since X I ii represents the particle position and the angular momenta for these particles are easily understood intuitively.

Summary of this work
We summarize our findings on the model (1.1) through the CLM and the minimum sensitivity analysis.
• We study the model with a finite angular momentum chemical potential through the CLM, and the results agree with the minimum sensitivity analysis quantitatively, as far as the chemical potential is not so large. It indicates that the CLM works properly in the rotating system. We also observe that the transition temperature for the confinement/deconfinement transition decreases as the chemical potential increases. See Sec. 4. This decreasing critical temperature is consistent with the previous studies in rotating black holes in holography [53][54][55][56][57][58][59] and rotating pure gluons [39,41].
• We develop an approximation method in the minimum sensitivity analysis, which enables us to explore the so-called "small black hole" solution [51]. 4 This solution has a negative specific heat similar to Schwarzschild black holes, and is important in the context of quantum gravity. See, for example, Fig. 2.
• By using the minimum sensitivity, we study the properties of the model with an imaginary angular momentum chemical potential, which has been employed to evade the sign problem in the MC computations [37][38][39][40][41]. We find a stable confinement phase at high temperature when D = 3. This is consistent with the recent study in the four-dimensional pure Yang-Mills theories [39]. We also argue a condition for the existence of the stable high temperature confinement phase in our model (1.1). Besides, we compute the critical temperature for the real chemical potential through the analytic continuation of the imaginary chemical potential, and find that the analytic continuation quantitatively works for finite chemical potentials. See Sec. 7.
3 Indeed, it is known that the dynamics of the model (1.1) is similar to that of the N = 4 supersymmetric YM theory (SYM) on S 1 β × S 3 . For example, this SYM theory shows a large-N confinement/deconfinement phase transition related to the model (1.1) [48,49]. Particularly, the SYM theory at strong coupling has a gravity dual given by AdS geometries [50,51]. There, at low temperature, a thermal AdS geometry is stable while an AdS black hole geometry is favored at high temperature, and a phase transition between these two geometries is called the Hawking-Page transition [52]. These geometries correspond to the confinement and deconfinement phases in the SYM theory, and hence they are also related to those of the model (1.1). Note, however, that these phase structures differ from the supersymmetric D-particles in the type-IIA superstring theory [11]. In particular, they do not show confinement at low temperatures [44]. 4 Analyses of the model (1.1) through the minimum sensitivity were also done in Refs. [26] and [28]. Particularly, Ref. [28] confirmed the existence of the small black hole solution in the model. However, this study merely pointed out the existence of the solution, and the solution itself was not derived. They also could not derive the gapped solution, which corresponds to a large black hole. The derivation of these solutions through the new approximation is one of our results in this article.
gapped distribution Figure 1: Three typical eigenvalue distributions ρ(α) of A t (2.4). The uniform distribution (left) is favored at low temperatures, while the gapped distribution (right) is favored at high temperatures. The non-uniform distribution might appear at middle temperatures.
This paper is organized as follows. In Sec. 2, we review the thermodynamical properties of the model (1.1) without angular momentum (readers familiar with this topic can skip this section). In Sec. 3, we introduce angular momentum to the model (1.1), and show that the model has a sign problem. In Sec. 4, we present our main results that the CLM successfully predicts non-trivial phase structure of the model (1.1), which agrees with the minimum sensitivity analysis. In Sec. 5, we present the details of the application of the CLM to the model. In Sec. 6, we explain the minimum sensitivity analysis. In Sec. 7, we study the imaginary chemical potential in our model by using the minimum sensitivity, and compare the results obtained from four-dimensional YM theories. Sec. 8 is devoted to a discussion.
In this article, we take an unit c = = k B = 1. For numerical computations, we take λ = 1.

Review of the previous results on the non-rotating model
In this section, we briefly review the properties of the model (1.1) without rotation. At finite temperatures, the model has a large-N confinement/deconfinement transition. The order parameters of this transition are the Polyakov loop operators If u n = 0 ( ∀ n), it indicates a confinement, and u n = 0 ( ∃ n) signals a deconfinement. In order to investigate the phase transition, it is useful to take the static diagonal gauge Here α k (k = 1, · · · , N) is independent of the Euclidean time ∂ t α k = 0. It also satisfies α k = α k + 2π. Then, the Polyakov loops u n become We also introduce the density function for {α k }, As we see soon, the profile of this function characterizes the phases of our system. Roughly speaking, the phase structure of this system can be explained as follows. If the scalars X I are not present, "repulsive forces" work between {α k }, and they uniformly distribute on the configuration space that is a circle (α k = α k + 2π). Then, the density function ρ(α) becomes ρ(α) = 1/2π as plotted in Fig. 1 (left). We call this solution a uniform solution. 5 From Eq. (2.4), u n = 0 is satisfied, and the system is in a confinement phase.
Once we turn on the scalars X I , they induce "attractive forces" between {α k }, which strengthen as temperature increases. Thus, while the system remains confined at sufficiently low temperatures, the attractive forces would overcome the repulsive forces at high temperatures, and {α k } collapse to a cluster. Then, the density function is depicted as shown in Fig. 1 (right). There, a gap exists in the distribution of {α k } around α = π, and this solution is called a gapped solution. 6 In this solution, u n = 0 and it corresponds to a deconfinement phase.
There is another solution that connects the uniform solution and the gapped one as shown in Fig. 1 (center). In this solution, any gap does not exist, and it is called a non-uniform solution. We can easily check that u 1 = 0 in this configuration, and it is in a deconfinement phase. (Hence, we have two deconfinement phases in large-N gauge theories.) Between these three solutions, large-N phase transitions may occur. Note that, whether the non-uniform solution arises as a stable phase depends on the dynamics of the system, and, for example, it is always unstable at D = 9 as we will see soon.
In this way, the model (1.1) has a confinement/deconfinement transition similar to higherdimensional YM theories. One feature of this transition is that the order of the transition would change depending on the dimension D [19]. For small D, it is first order, while it would be second order for large D. Indeed, MC computations have established that it is first order up to D = 25 [19,22]. On the other hand, the 1/D-expansion predicts that the transition is second order at sufficiently large D [31]. Also, the minimum sensitivity analysis at three-loop order indicates that it is first order up to D = 35 and it becomes second order from D = 36 [28]. (Note that the minimum sensitivity analysis at two-loop order predicts the phase transition is always first order independently of D. Also, it is not clear whether the three-loop computation is sufficiently reliable, and it has not been established if the order of the transition changes at D = 36.) Such a D-dependence of the transitions is also consistent with a holography. 7 5 As we will see, in the minimum sensitivity analysis, phases of the model (1.1) are obtained as saddle point solutions of an effective action. Hence, we may call phases "solutions". 6 In this solution, we have taken a gauge that the peak of the density ρ(α) is at α = 0. Then, all u n would become real. We use this gauge throughout our minimum sensitivity analysis. However, in the MC and CLM, taking this gauge is difficult, and we evaluate |u n | instead of taking this gauge. Note that |u n | and u n in the gauge are different in general, but the deviations are highly suppressed at large N , and we do not observe any issue when we compare these two quantities. 7 As we mentioned in the introduction, the confinement/deconfineement transition of the model (  In this article, we will study the model with a finite angular momentum chemical potential through the CLM. There, we investigate D = 9 and 16, where the first order transitions occur. Hence we will employ the minimum sensitivity analysis at two-loop order rather than the 1/D-expansion to test the CLM, since it also predicts the first order transitions. (The minimum sensitivity analysis at three-loop order with the finite chemical potential is much more complicated than the two-loop analysis and we do not do it in this article.) In order to test the minimum sensitivity analysis at two-loop order, we compare the results at zero chemical potential obtained through this analysis and MC (not CLM), and we find quantitative agreements. 8 (We will show the details of the minimum sensitivity analysis in Sec. 6.) See Fig. 2 for u 1 at D = 9. Hence, we expect that the minimum sensitivity analysis at two-loop order may be reliable even with a finite chemical potential, and we use this method to test the CLM.
We also plot the temperature dependence of free energy through the minimum sensitivity at D = 9 in Fig. 2 (left). It shows a first-order phase transition as we mentioned above. related to a Gregory-Laflamme (GL) transition in gravity [45]. Here, D = 9 is taken and the dual gravity theory is in Euclidean S 1 β × S 1 × R d space with d = 8. We can formally extend this correspondence to a general D by taking d = D − 1 [19]. Then, the order of the GL transition depends on the dimension D. It is first order for D ≤ 11 and it is second order for D ≥ 12, and they are qualitatively similar to our matrix model results [60,61]. Interestingly, a similar D-dependence has been observed in Rayleigh-Plateau instability in a fluid model, too [62,63]. 8 The agreements become not so good as temperature increases. This is because we use an approximation (A.22), and it is not reliable at higher temperatures.

Introducing chemical potentials for angular momenta
We introduce angular momentum chemical potentials to the model (1.1). For simplicity, we first consider the chemical potential in a single particle quantum mechanics in two dimensions. We take the Hamiltonian of this system as Here, x and y are the two-dimensional position operators, and p x and p y are their conjugate momenta. It is convenient to employ a complex coordinate z := (x + iy)/ √ 2. Then, the angular momentum is given by Now we introduce the angular chemical potential µ, and the Hamiltonian is modified as Then, the the Euclidean action for this Hamiltonian, which is used in the path-integral formalism at finite temperature, is derived as Here, the first term is complex. Hence, a finite angular momentum causes a sign problem in MC in general. Let us introduce the angular momentum chemical potentials to the matrix model (1.1) through the same procedure. Note that the model (1.1) is a quantum mechanics in D dimension, and we consider rotations onD planes (2D ≤ D). Hence, we introduce the chemical potential µ I (I = 1, · · · ,D) for each plane. Then, the action becomes [32,59] where we have defined Z I := X I + iXD +I / √ 2, (I = 1, · · · ,D). In the following analysis, for simplicity, µ I is taken to be a common value µ > 0 to all µ I .

Overview of our main results
In this section, we show our main results for the model (3.5) derived through the CLM and the minimum sensitivity. The details of the CLM and the minimum sensitivity are presented in Sec. 5 and Sec. 6, respectively. We mainly investigate D = 9 and D = 16. We have chosen them because D = 9 is the critical dimension of superstring theories and D = 16 is suitable to be compared with the minimum sensitivity, since they agree better for larger D in the µ = 0 case as demonstrated in Ref. [28].

Phase diagrams
We draw the µ − T phase diagrams at D = 9 withD = 1 andD = 3 in Fig. 3. These are obtained through the minimum sensitivity analysis (see Sec. 6.1.4). These phase diagrams show that the uniform phase is favored in a low temperature and chemical potential region, and the gapped phase is favored in a high temperature and chemical potential region. A first-order transition occurs between them. In the "unknown" region depicted in Fig. 3, the minimum sensitivity analysis does not work, and we presume that this region may be thermodynamically unstable due to a large chemical potential. (As far as we attempt, we obtain similar phase diagrams when we change D and/orD.) These phase diagrams show that a larger µ or a largerD makes the transition temperature lower. Similar properties are expected in four-dimensional pure YM theories and black holes as we will argue in Sec. 7 and Sec. 8.1.

Results for observables
By using the CLM and the minimum sensitivity, we investigate µ dependence of the four observables: the Polyakov loop |u 1 | and |u 2 |, the angular momentum J and the expectation values of the square of the scalars (X I ) 2 . The results are shown in Figs. 4, 5, 6 and 7, respectively. There, we take D = 9 and D = 16, and investigate them by changingD and temperature T . The values ofD and T that we have taken are summarized in Table 1. In this table, the phase transition point µ * computed by the minimum sensitivity, at which the uniform phase turns into the gapped phase, is also listed.
In the CLM, N=16 and 32 are taken at D = 9, and N = 16 is taken at D = 16. Note that we present the numerical results obtained by the CLM only in the parameter region of µ where the data are acceptable. (The criterion for the acceptable data is based on the drift norms argued in Sec. 5.) We find that obtaining acceptable results for larger µ is harder.
In the following subsequent subsections, we argue the details of the obtained observables. Our results for |u 1 | and |u 2 | against µ are plotted in Figs. 4 and 5. The minimum sensitivity predicts the first-order transition from the uniform phase (u 1 = u 2 = 0) at low µ to the gapped phase (u 1 = 0 and u 2 = 0) at high µ. The CLM results agree with the minimum sensitivity results at large N in the gapped phase. In the uniform phase, |u n | has a finite-N effect of order O 1 N [19], and we compare the CLM results of |u n | with the minimum sensitivity result with the 1/N correction (6.25), which also agree with each other.
Among our results, the deviation between the CLM and the minimum sensitivity at T = 0.90 in the D = 9 withD = 1 case is larger. We presume that the fluctuation near T = 0.90 is large in the CLM, since the phase diagram in Fig. 3 indicates that T c ≃ 0.9 for µ 1.
Although we observe the phase transitions in the CLM, we do not attempt to determine their order. In order to do it, we may need to evaluate, for example, the susceptibility for u 1 [19], but it requires computation at larger N in the CLM, and we leave it for a future work.

Angular momentum J
Our results for angular momentum J is shown in Fig. 6. Note that we have introduced the common angular momentum chemical potential µ for theD planes in the model (3.5), and we evaluate the angular momentum for the single plane. The definition of J in the CLM is given in Eq. (5.12). In the minimum sensitivity, we derive J through Eq. (6.28).
We observe that J is close to zero in the uniform phase in the CLM. This is a feature of the confinement phase (see Sec.6.2.2), and the CLM correctly reproduce it. In the gapped phase, we observe slight discrepancies between the CLM results and the minimum sensitivity results. We investigate the lattice space dependence and find that J is more sensitive than other quantities. Thus, the discrepancies may include artifacts of the CLM.

Expectation values of scalars
The square of the scalar (X I ) 2 represents the distribution of the D-particles on the I-th direction. Since we have the rotational and the non-rotating directions, we define the averages and investigate them separately. (We have assumed that any spontaneous symmetry breaking on the rotation direction or the non-rotating directions does not occur. Indeed, we do not observe any signal for it in the CLM. We have also assumed that they are time independent.) In Fig. 7, the results of R 2 Z and R 2 X are plotted. The discrepancies between the CLM results and the minimum sensitivity results are very small in the uniform phase and a few percent in the gapped phase. They are getting worse as µ increases, where the approximation (A.22) in the minimum sensitivity is not reliable. Again, the discrepancy is larger in the D = 9 andD = 1 at T = 0.90 case due to the phase transition.
The both analyses capture important features of the rotating system. In the uniform phase, there is no separation between R 2 Z and R 2 X . In the gapped phase, R 2 Z gains a larger value, and R 2 Z and R 2 X begin to separate with each other. It means that the D-particles spread to the rotation directions as J increases, and it is a natural consequence as rotating objects.        Z and R 2 X (4.2) are plotted against µ. We present R 2 Z and R 2 X obtained by calculating (5.13) and (5.14) via the CLM respectively, at N = 16, 32 for D = 9 and N = 16 for D = 16. The lines denote the results of the minimum sensitivity (6.29) and (6.30) at N = ∞. The dashed, dotted and solid lines represent those of the uniform, non-uniform and gapped phase, respectively.

15
In this section, we present the details of the application of the CLM [3,4], which is a promising method to simulate the system with sign problem, to the action (3.5), and explain the derivation of the results shown in Sec. 4.

Complex Langevin equation
In solving the complex Langevin equation for the action (3.5), we rescale it as where λ = g 2 N is the 't Hooft coupling. This gives . We omit ′ in the following. At µ = 0, this action is invariant under the transformations 4) where I N is an N × N unit matrix. To put the action (5.1) on a computer, we adopt a lattice regularization, where the number of the lattice sites is n t and the lattice space is (∆t) = β nt . We also adopt the periodic boundary condition X µ (n t + 1) = X µ (1) and A(n t + 1) = A(1). This yields the lattice-regularized action where V (n) = e iA(n)(∆t) . We find it convenient to take a static diagonal gauge (2.2). In this gauge, we have Together with the gauge fixing term, as derived in Refs. [49,64,65], we work on the action The CLM consists of solving the complexified version of the Langevin equation. The complex Langevin equation is given by Here, σ is the fictitious Langevin time, and the white noises η I kℓ (n, σ) and η where σ 0 is the thermalization time, and σ T is the time required for statistics, both of which should be taken to be sufficiently large. In Refs. [5,9,66], it was found that the holomorphy of the observable O plays an essential role in the validity of Eq. (5.9). When we solve the Langevin equation (5.8) with a computer, we need to discretize it as X I kℓ (n, σ + ∆σ) = X I kℓ (n, σ) − (∆σ) ∂S eff ∂X I ℓk (n, σ) Here, ∆σ is the step size, which we take to be ∆σ = 10 −5 . The factor √ ∆σ stems from the normalization of the discretized version of the white noisesη I kℓ (n, σ) andη

Observables
In the CLM, we calculate the Polyakov loop u n , which is defined as Eq. (2.3), and the following observables. 9 In a rotating system, we cannot apply an angular momentum to the center-of-mass, which leads us to remove the trace part and implement the constraints 1 Nβ We solve the Langevin equation (5.8) without imposing the constraints (5.15), which facilitates the numerical calculation. By using these methods, we perform the numerical computations at n t = 60 and obtain the results shown in Figs. 4-7. As mentioned in Sec. 4, these results nicely agree with the minimum sensitivity results for finite µ. However, for larger µ, we encounter several troubles in the CLM. We present them in the next subsection.

Testing the validity of our CLM
The CLM faces the following two typical problems. One is the "excursion problem", which occurs when X I and α k are far from Hermitian matrix and real number, respectively. 10 The 9 In the CLM, if observables are not holomorphic, we may not obtain correct answers. A subtle observable in our analysis is |u n | which is not holomorphic with respect to {α k }. However, at large N , |u n | 2 = |u n | 2 = u n u −n is held for real {α k }, and, particularly, u n u −n is holomorphic. It turns out that this relation is approximately satisfied in our CLM in the relevant parameter region where we accept the CLM result based on the criterion described in Sec. 5.3. Hence, evaluating |u n | would not be problematic. 10 The gauge cooling [6][7][8] is a standard technique to suppress the excursion problem. However, in our complex Langevin studies we already fix the gauge as Eq. (2.2), which prevents us from applying the gauge cooling.
other is the "singular drift problem", which occurs when the drift terms are too large. It is found in Ref. [9] that a sufficient condition to justify the CLM is that the probability distribution of the drift norms fall off exponentially or faster. If we look at the drift term, we get the drift of the CLM, and we can easily test this criterion. We work on the D = 9,D = 1, 3, T = 0.85, 0.90, N = 16, 32 and D = 16,D = 1, 5, T = 0.80, 0.85, N = 16 cases. In these cases, we take n t = 60. To probe the region of µ we can study by the CLM, we present the log-log plots of the probability distribution of the drift norms u X and u α , which we denote as p(u X ) and p(u α ), respectively. As a typical example, we show the D = 9,D = 1, T = 0.90 case in Fig. 8. p(u X ) falls exponentially or faster for µ ≤ 1.2 at N = 16 and µ ≤ 0.9 at N = 32, respectively. On the other hand, p(u α ) falls in a power law even for small µ. As we present in Appendix B, the power-law decay of p(u α ) is observed even at µ = 0, which has no sign problem. At µ = 0, without the static diagonal gauge (2.2), the probability distribution p(u A ), where u A is the drift norm without the static diagonal gauge as defined by Eq. (B.7), falls exponentially or faster. At µ = 0, we confirm the agreement of the observables between the cases with and without the static diagonal gauge (2.2). Also, in Ref. [67], the static diagonal gauge (2.2) is taken to study the unitary matrix model, and the consistency between the CLM and analytic results is reported despite the power-law behavior of the drift norms as presented in Appendix B. Hence, we presume that the power-law decay of p(u α ) in the static diagonal gauge is harmless (Recall that the criterion in Ref. [9] is a sufficient condition).
To save the CPU time, we take the static diagonal gauge (2.2) and accept the results for µ ≤ 1.2 at N = 16 and µ ≤ 0.9 at N = 32, where p(u X ) falls exponentially or faster while p(u α ) does not. When µ is so large as to be outside the parameter region where we accept the CLM result, the simulation gets unstable and crashes. Similar trends are observed for other D,D and T . In Figs. 4 -7, we present the numerical results in the parameter region of µ, where we accept the CLM result with this criterion.
In the CLM, the hermiticity of X I (n) and the reality of α i are lost, and the observables (5.12), (5.13) and (5.14) are not real in general. In the parameter region where we accept the CLM result, we present their real part of the expectation values of the observables (5.12), (5.13) and (5.14) obtained by the CLM, as the ensemble average of their imaginary part turns out to be close to 0. Also, the ensemble average of the imaginary part of α i turns out to be close to 0. These disappearances of the imaginary parts indicate the validity of our analysis.

Derivation of Free energy
We study the thermodynamical properties of the model (3.5) at large N via a minimum sensitivity analysis [68]. For this purpose, we introduce trial masses m Z and m X for Z I and X I , respectively, and deform the action (3.5) as [28] S κ :=S 0 + κS int , (6.1) Here κ is a formal expansion parameter. If we set κ = 1, the mass terms are canceled and this action reproduces the original one (3.5). However, if we perform a perturbative expansion with respect to κ up to a certain order and take κ = 1 after that, then the obtained quantities would depend on the trial mass parameters m Z and m X . 11 The idea of the minimum sensitivity is that we fix m Z and m X so that the dependence of a certain physical quantity on these parameters is minimized. It has been demonstrated that this prescription works in various models 12 . Note that the obtained result through this method would depend on of which quantity we minimize the parameter dependence. In our study, we investigate free energy at two-loop order, and minimize its m Z and m X dependence. By integrating out X I and Z I perturbatively, we obtain the effective action for {α k } at two-loop order as shown in Appendix A.1, Here we have taken κ = 1 and dU is an integral measure defined by [49] dU := k dα k e −S g.f. , In order to derive the free energy F , we need to evaluate the dU integral in Eq. (6.4). This integration at large N has been studied in Ref. [69], and the result depends on the sign of f 1 . If f 1 > 0, the saddle point for the uniform solution ( Fig. 1 [left]) dominates, while when f 1 < 0, the saddle point for a gapped solution (Fig. 1 [right]) does, and the free energy is given by [69] βF (T, µ, m X , m Z ) = (6.6) (6.7) Here we have defined Note that we treat f 1 = 0 case separately as we will explain in Sec. 6.1.3. Correspondingly, the Polyakov loop (2.1) at large N is computed as [70] (6.9) (6.10) Here, we have taken the gauge mentioned in footnote 6. Similarly, for u n (n ≥ 2), we obtain [71,72] where P (1,2) n−2 (z) denotes the Jacobi polynomial. These results show that, when f 1 > 0, u n = 0 for all n and the density function ρ(α) becomes uniform through Eq. (2.4). Hence the system is confined. On the other hand, when f 1 < 0, the gapped solution satisfies u n = 0 and the system is deconfined.
We have so far evaluated the dU integral in the partition function (6.4). Now we fix the trial masses m X and m Z so that the dependence of the free energy F (T, µ, m X , m Z ) on these masses is minimized. Hence, we solve the equations In the subsequent subsections, we evaluate these equations in each f 1 > 0, f 1 < 0 and f 1 = 0 cases. It will determine the free energies of each solution, and they tell us their stabilities and phase structures as drawn in Fig. 3.

Uniform solution
In this subsection, we evaluate Eq. (6.13) when f 1 > 0, which is in the confinement phase. In this case, the free energy becomes F = N 2 f 0 through Eq. By using this result, the free energy is given by Hence, the free energy in this solution depends on neither temperature nor the chemical potential. See Figs. 2 and 10 for the D = 9 case. Note that, when we derived the effective action (6.4), we had assumed m Z < µ. Hence, the uniform solution is not reliable in the region µ ≥ m 0 . For example, in the D = 9 with D = 3 case shown in Fig. 3, such a region appears in the uniform phase, and the fate of the system there is not obvious. It is likely that the system is unstable in this region and any stable phase does not exist.

Gapped solution
We discuss the f 1 < 0 case. Here, we cannot solve Eq. (6.13) analytically, and we evaluate it numerically. For a fixed T and µ, we may find several solutions of m X and m Z . See   9 for D = 9 withD = 1 andD = 3. However, the condition m Z > µ had been assumed when we derived the effective action (6.4), and the solutions that do not satisfy it are not reliable. As far as we investigated, only the solutions connected to the non-uniform solution at µ = µ GWW (T ), which we will argue in the next subsection, are reliable. As we increase µ with a fixed T , even these reliable solutions reach the point m Z = µ, which we call µ = µ unstable , and the fates of the systems beyond it are unclear. These regions are presented as "unknown" in the phase diagrams in Fig. 3. We presume that the systems are unstable in these regions.

Non-uniform solution
We have seen that the uniform solutions and the gapped solutions appear when f 1 > 0 and f 1 < 0, respectively. Thus, a transition would occur at f 1 = 0, and we investigate this case in details. For this purpose, it is useful to rewrite the partition function (6.4) as Here, we have changed the integral variables from {α k } to {u n }. (Such a change is possible in the large-N limit [49].) Note that {u n } are not completely independent, since they have to satisfy the condition that the density function ρ (2.4) is non-negative. shows the first-order phase transition between the uniform and gapped phase at µ * . Our analysis for the gapped solution is reliable in the region µ < µ unstable , where the condition µ < m Z is satisfied. The purple curve is for the gapped solution shown in Fig. 9, although it appears in the unreliable region.
The action (6.16) is quadratic in {u n }, and their coefficients are all positive for n ≥ 2. Thus, u n = 0 (n ≥ 2) is a stable solution. Hereafter, we assume u n = 0 (n ≥ 2) and focus on u 1 .
By differentiating the action (6.16) by u 1 , we obtain an equation Obviously, one solution is given by u * 1 = 0, which represents the uniform solution studied in Sec. 6.1.1. The other possible solution is f 1 (T, µ, m X , m Z ) = 0. (6.18) This is what we are interested in. Besides, we have the conditions (6.13) that the dependence of the free energy on m X and m Z is minimized, By combining these equations and Eq. (6.18), we obtain three equations, The last two equations determine m X and m Z , and the first one fixes u 1 . These equations can be solved numerically, and the results for D = 9 withD = 1 andD = 3 are shown in Fig. 9.
On the other hand, the non-uniform solution also merges to the uniform one at u 1 = 0. Through Eqs. (6.14) and (6.18), it occurs when µ and T satisfies We define this solution as µ c (T ), and call it a critical point. Therefore the non-uniform solutions exist in the region µ GWW (T ) ≤ µ ≤ µ c (T ). See Fig. 2 and 10, again. Note that f 1 may be negative beyond the critical point µ = µ c (T ). Since f 1 is the coefficient of |u 1 | 2 in the effective action (6.16), the uniform solution becomes unstable in this case.

Free energy and phase structure
So far, we have obtained the three solutions: uniform, non-uniform and gapped solution. To see the phase structure, we evaluate their free energies. For the uniform solution, we have derived the free energy in Eq. (6.15). For the non-uniform and gapped solutions, we obtain their free energies by substituting the numerical solutions m X and m Z into Eqs. (6.16) and (6.7). These results are summarized as Note that we have used f 1 = 0 and u n = 0 (n ≥ 2) for the non-uniform solution. The free energy for the D = 9 withD = 1 case at T = 0.80 is plotted in Fig. 10. This figure shows that a first-order transition occurs in this system. There, the transition point is given when the free energy of the uniform solution and the gapped one are coincident. We define µ * (T ) for this point. 13 Fig. 10 also shows that the free energy of the non-uniform solution is always higher than those of the uniform solution and the gapped one. Actually, the free energy of the non-uniform solution is a concave function, and the specific heat is negative. These results imply that the non-uniform solution is always unstable in the grand canonical ensemble.
One feature of this phase transition is that neither the non-uniform solution nor the gapped one exists in the region µ < µ GWW (T ). 14 This is due to our approximated effective action (6.16), and, if the action involves higher-order terms such as u n u m u −n−m , the location of the GWW transition point µ GWW (T ) would change [70].
By combining all the results, the whole phase diagrams are obtained as drawn in Fig. 3.

Calculating observables
We have investigated the phase structure of the model. Now, we explain the derivation of the observables shown in Figs. 4 -7 in Sec.4.2. The results in this section are for the large-N limit, unless it is specified.

Polyakov loops
We evaluate the Polyakov loop operators (2.1), which are the order parameters of the confinement/deconfinement transition. In the uniform solution, u n = 0 is obtained through Eqs. (6.9) and (6.11) at large N. We can also derive the leading 1/N correction (A.23) as discussed in Appendix A.2. The result is given by , (non-uniform solution), (6.26) at large N as argued in Sec. 6.1.3. Here, m X and m Z are the solutions of Eq. (6.21).
The results for u 1 are plotted in Figs. 2 and 4, and u 2 is shown in Fig. 5. As we have seen in Sec. 4.2.1, they are consistent with the CLM.

Angular momentum
We have derived the free energy (6.24) in Sec. 6.1.4. By using this result, we can read off the angular momentum via (6.28) (As we mentioned in Sec. 4.2.2, we calculate the angular momentum for the single plane. Hence, we have divided −(∂F/∂µ) byD.) The results are compared with the CLM as shown in Fig. 6. Interestingly, J decreases as µ increases in the non-uniform solution. This property would be related to thermodynamical instabilities of the non-uniform solution.
Besides, J = 0 in the uniform phase, because the free energy (6.24) does not depend on µ, there 15 . Thus, the uniform phase does not rotate at large N, although the chemical potential is finite. (Similarly, entropy in the uniform phase is zero even at finite temperatures. It means that thermal excitations are highly suppressed in the uniform phase indicating a confinement [48,49].)

Expectation values of scalars
We evaluate the expectation values of the scalars R 2 Z and R 2 X defined in Eqs. (4.1) and (4.2). Through the one-loop computation (A.18), we obtain To compute them, the suitable solutions for m X , m Z and u n need to be substituted. The results are plotted in Fig. 7.

Imaginary chemical potentials and relation to rotating YM theory in four dimensions
Rotating quark gluon plasma (QGP) is actively being studied motivated by relativistic heavy ion colliders [42] and neutron stars. As a related problem, rotating pure YM theories are also being investigated. Particularly, one important question is how the rotation affects the confinement/deconfinement transition temperatures. (Rotating media are non-uniform in space, and the transition temperatures on the rotation axis is mainly studied.) However, these theories are strongly coupled and standard perturbative computations do not work. In addition, the sign problem prevents lattice MC computations.
To avoid these issues, the imaginary angular velocity µ Im ∈ R is considered [37][38][39][40][41]. (This imaginary angular velocity µ Im corresponds to the angular momentum chemical potential µ in our model (3.5) as µ = iµ Im through the dimensional reduction, and we call µ Im a imaginary chemical potential, hereafter.) The imaginary chemical potential does not cause the sign problem and MC works. Once we obtain the results for the imaginary chemical potential, through the analytic continuation, we may reach the results for the real chemical potential. Such analytic continuation would work as far as the chemical potential is sufficiently small.
In this section, we review some results in pure YM theories with the imaginary chemical potential. Then, we compute the corresponding quantities in the matrix model (3.5) through the minimum sensitivity, and compare them with the YM theories. We will see some similarity between the matrix model (3.5) and the YM theories, and the matrix model provides some insights into the YM theories.

Stable confinement phase at high temperatures
Recently, one remarkable result on the SU(3) pure YM theory was reported in Ref. [39]. The authors investigated the high temperature regime (T → ∞), where the perturbative computation is reliable, and found that the system is in a confinement phase for π/2 ≤ βµ Im ≤ 3π/2 and in a deconfinement phase for 0 ≤ βµ Im ≤ π/2 and 3π/2 ≤ βµ Im ≤ 2π. Hence, the system is confined, although temperature is high. See Fig. 4 in Ref. [39]. (Note that, when the imaginary chemical potential is turned on, the Boltzmann factor is multiplied by exp(iβµ Im J), and the thermal partition function is periodic with respect to βµ Im : Z(βµ Im ) = Z(βµ Im + 2π).) Then, one important question is whether this high temperature confinement phase in π/2 ≤ βµ Im ≤ 3π/2 continues to the low temperature confinement phase at βµ Im = 0. To answer this question, we need to study the strong coupling regime in the YM theory, and it has not been understood.
This result motivates us to study the imaginary chemical potential in our matrix model (3.5). Particularly, exploring the high temperature regime and investigating the fate of the confinement phase at T → ∞ would be valuable.
The analysis through the minimum sensitivity is almost straightforward. We simply need to repeat the same computations done in Sec. 6 by using the effective action (6.4) with µ = iµ Im . To investigate the phase structure at high temperatures, we evaluate the critical point (6.23)  Interestingly, in the D = 3 andD = 1 case that is the dimensional reduction of the four-dimensional YM studied in Ref. [39], the system is confined in π/2 ≤ βµ Im ≤ 3π/2. This is the same result as that of Ref. [39], although SU(3) is taken in Ref. [39] and we have taken the large-N limit. Therefore, our model might capture the phase transition of the original model.
We also derive the whole phase structures in the D = 3 withD = 1 case and the D = 9 withD = 1 case as drawn in Fig. 11. In the D = 9 withD = 1 case, the condition (7.2) is not satisfied, and the high temperature confinement phase does not appear. In the D = 3 withD = 1 case, we observe that the high temperature confinement phase continues to the conventional confinement phase at µ Im = 0. This suggests that the confinement phase may be continuous in the four-dimensional YM theory, too.

Analytic continuation of the chemical potential
Refs. [39,41] argued that the imaginary chemical potential increases the transition temperature in the pure YM theory. 17 Through the analytic continuation µ → iµ Im , this result implies that the decreasing transition temperature under the presence of the real chemical potential at least for small µ. However, it is unclear whether such an analytic continuation provides quantitatively good results for a finite chemical potential. Thus, it may be valuable to test the analytic continuation in our matrix model, since the transition temperatures for each the real and imaginary chemical potential can be computed.
In Fig. 12, we explicitly compare the critical temperature T c against the real chemical potential µ computed through Eq. (6.23) and that of the analytic continuation of the imag-  Fig. 11. In both cases, they agree for µ 0.5 × µ c | T =0 , and the analytic continuation works there.
inary chemical potential derived in Fig. 11. 18 In both D = 3 and D = 9 cases, we observe good agreement for µ 0.5 × µ c | T =0 . This result suggests that the analytic continuation of the imaginary chemical potential may work in similar ranges even in the four-dimensional YM theories and QCD, too.

Discussions
In this article, we have studied the matrix model (1.1) at finite angular momentum chemical potentials by using the CLM and the minimum sensitivity, and found the quantitative agreements. The action (3.5), with the chemical potential for the angular momentum added, suffers from the sign problem. This prevents us from using the conventional Monte Carlo methods, as we cannot regard e −S for the complex action S as a probability. This leads us to study the model numerically (3.5) using the CLM. The CLM turns out to work successfully in the parameter region of µ, wide enough to elicit the behavior of the confinement and deconfinement phases. While minor discrepancies between the results of the CLM and minimum sensitivity are observed in the results presented in Figs. 4 -7, they would be mitigated with more lattice space than n t = 60 and higher loop corrections (A.1) and higher order corrections (A.21) in the minimum sensitivity treatment. This is the very first result showing that a rotating quantum many body system at thermal equilibrium is analyzed through the first principle computation, as far as the authors know. Such rotating quantum systems are important in various topics including condensed matter 18 To obtain the analytic continuation results, we plot T c against µ Im by using the data employed in Fig. 11, and fit the obtained curve by a polynomial T c = n c n (µ 2 Im ) n , where c n are the fitting parameters. Then, we perform the analytic continuation and obtain T c = n (−1) n c n µ 2n . This is the red curve plotted in Fig. 12. and high energy physics, and our result encourages us to apply the CLM to these systems, too.
In Sec. 7, we have compared our matrix model and rotating pure YM theories in four dimensions under the presence of the imaginary chemical potentials. We found the stable confinement phase at high temperature in the matrix model akin to the YM theories argued in Ref. [39]. We also found that the increasing transition temperature consistent with Refs. [39,41]. Therefore, the natures of the matrix model (1.1) is quite similar to the YM theories. Since we can investigate the model (1.1) with the real chemical potential, it may provide insights into the YM theories. Besides, if we apply the CLM to the pure YM theories, it may shed more light on the properties of the YM theories.

Relation to gravity
We have seen that the transition temperature decreases as the chemical potential increases in our model. A similar result has been obtained in the N =4 SYM on S 1 β × S 3 [58,59]. (See also footnote 3.) This model at strong coupling would be described by dual AdS geometries through the AdS/CFT correspondence [43], and the gravity computation [53][54][55][56][57] also indicates a similar phase structure. (See Fig. 1 of Ref. [59].) There, rotating black D3brane solutions correspond to the rotating gapped solutions in the SYM theory. 19 Therefore, the decreasing transition temperature by a rotation is a common feature of these gauge theories and gravities.
In relation to the black holes, the non-uniform solution derived through the minimum sensitivity analysis is interesting. As shown in Sec. 6.1.4, this solution has the negative specific heat akin to black hole solutions such as Schwarzschild black holes and small black holes in AdS correspondence [51]. Hence, the non-uniform solution may explain the origin of the negative specific heat of the black holes through microscopic description. It would be valuable to pursue this question further.
Besides, the properties of the large chemical potential regions (the "unknown" regions in Fig. 3) in the model (3.5) may be understood through the gravity. We can compute the free energy of the black branes as a function of temperature and chemical potential by using the results of Ref. [55]. Then, we will see a similar result to our result shown in Fig. 10: Two black brane solutions appear and they merge at a large chemical potential, and, beyond this point, there is no solution. (The two black brane solutions correspond to the two gapped solutions in the matrix model.) Although our result in Fig. 10 beyond µ = µ unstable is not reliable, the gravity analysis indeed predicts a similar result. It is tempting to improve our approximation in the matrix model and verify the gravity prediction. In addition, gravity systems with angular momentum have various exotic solutions such as black rings and black Saturn, and it might be possible to find the corresponding solutions in the matrix model, too.
in part by Grant-in-Aid for Scientific Research C (No. 20K03946) from JSPS. Numerical calculations were carried out using the computational resources, such as KEKCC and NTUA het clusters.
A Details of the minimum sensitivity analysis A.1 The derivation of the effective action (6.4) We show the derivation of the effective action (6.4) through the minimum sensitivity analysis. We will integrate out X I and Z J through a perturbative calculation in the deformed action (6.1) with respect to κ, and will obtain the two-loop effective action for the Polyakov loop {u n }, To perform the perturbative calculation, we take the static diagonal gauge (2.2). Then the propagators of X I and Z I in the free part of the deformed action (6.2) become [31,32], where x := e −βm X , z := e −βm Z , q := e −βµ and u i n := e inα i , which satisfies N i=1 u i n = Nu n from Eq. (2.1). ||t|| denotes ||t + nβ|| = t for 0 ≤ t < β. In this calculation, |µ| < m Z has been assumed, otherwise Z I becomes unstable.
By using these propagators, we obtain the one-loop term in the expansion (A.1), [28,31,49] Here S G.F. is the gauge fixing term (6.5), which arises when we take the constant diagonal gauge (2.2). In the two-loop computation, we evaluate Here, the first term can be expanded as In the large-N limit, the planar diagram depicted in Fig. 13 dominates and each term in Eq. (A.7) is calculated as follows.
The terms involving the products of the propagators can be computed by using the expres- 2x n + z n (q n + q −n ) |u n | 2 + x m+n u m−n u −m u n + ∞ m,n=1 x m+n (u m+n u −m u −n + u −n−m u m u n ) . (A.14) By substituting these equations into Eq. (A.7), we reach where B n has been defined in Eq. (A.17). In principle, {u n } can be integrated out in this effective action by using the technique proposed in Ref. [75], and we may obtain the free energy. However, the computation would be complicated. In order to reduce the calculations, we use the following drastic approximation, and ignore the interaction term O(u n u m u −n−m ). In this approximation, the contribution of X I and Z I integrals to f n (n ≥ 2) are neglected, and only the terms coming from the gauge fixing (6.5) survive. Then, we obtain the effective action (6.4), which we use in the main analysis of the phase structure of the model (3.5). (The advantage of this approximation is that the analysis in Ref. [69] is available.) This approximation may be justified in the low temperature and low chemical potential regime until T ≃ T c (µ) by assuming that D is large. This is because m X = m Z ∼ (λD) 1/3 for a large D in the confinement phase (6.14), and the equation (6.23) 21 for T c (µ) leads to

B Drift norm of the bosonic BFSS model
In this section, we discuss the behavior of the drift norm in the real Langevin simulation of the bosonic BFSS model (1.1), which is free from the sign problem. In the following, we discuss the action (5.2), which is lattice-regularized as (5.5), with µ = 0. We compare the behavior of the drift norms, as defined by Eqs. (5.17) and (B.7), for the cases with and without the static diagonal gauge (2.2), respectively.

B.1 Langevin equation without the static diagonal gauge
Here, we describe how to solve the real Langevin equation for µ = 0 without the static diagonal gauge (2.2). The scalar fields X I are updated as (5.10), which makes no difference from the simulation with the static diagonal gauge (2.2). The gauge field now depends on the temporal direction, which is updated in terms of V (n) = e iA(n)(∆t) as V (n, σ + ∆σ) = exp i G a=1 λ a −(∆σ)ν a (V (n, σ)) + √ ∆ση a (n, σ) V (n, σ), (B.1) instead of (5.11). 22 Due to the invariance (5.4), we set the constraint TrA(t) = 0, and the gauge group is SU(N). λ a is the generator of the SU(N) Lie algebra, such that Tr(λ a λ b ) = δ ab , and G = N 2 − 1 is the dimension of SU(N). ν a (V (n, σ)) is the drift term defined by where τ is a real number and S lat [e iτ λ a V (n)] is defined by replacing V (n) and V (n) −1 in S lat , as defined by Eq. (5.5), with e iτ λ a V (n) and (e iτ λ a V (n)) −1 = V (n) −1 e −iτ λ a , respectively. Also, the drift norm for the case without the static diagonal gauge is defined as |ν a (V (n, σ))| 2 , (B.7) instead of u α in Eq. (5.17), while the drift term for X I is the same as u X in Eq. (5.17). 22 Here, we present the result µ = 0, which has no sign problem from the outset, and there is no need to implement the gauge cooling [6][7][8]. When we work on the µ = 0 case, the gauge cooling is useful to suppress the excursion problem by minimizing the hermiticity and unitary norm for X I and V defined by Here, H X (n) and H V (n) are the gradients of N X and N V with respect to the gauge transformation, and the real positive parameters γ X and γ V are chosen so that the norms N X and N V are minimized, respectively.

B.2 The fall-off of the drift norms
We compare the result with and without the static diagonal gauge (2.2), for D = 3, N = 16, µ = 0. With the static diagonal gauge, we add the gauge fixing term as (5.7). In Fig. 14, we compare the observables |u 1 |, |u 2 | and R 2 = 1 Nβ This is invariant under the gauge transformation V (n) → g(n + 1)V (n)g −1 (n). |u n | and R 2 are calculated by removing the trace part as (5.16). The trace part does not affect F 2 , which is expressed only in terms of the commutator of X I (t). We see that the results with and without the static diagonal gauge (2.2) agree. The histograms of the drift norms are plotted in Figs. 15 and 16 with and without the static diagonal gauge (2.2), respectively, for D = 3, N = 16, µ = 0, T = 0.5, 1.0, 1.5, 2.0. The probability distributions p(u X ) of the drift norms u X fall exponentially or faster both in Figs. 15 and 16 with and without the static diagonal gauge (2.2). However, those of the drift norms p(u α ) for the gauge field u α fall only in power law in Fig. 15 with the static diagonal gauge (2.2), while those of u A , which we denote as p(u A ), fall exponentially or faster without the static diagonal gauge (2.2). We attribute this to the singularity stemming from the derivative