Complex Langevin analysis of the spontaneous symmetry breaking in dimensionally reduced super Yang-Mills models

In recent years the complex Langevin method (CLM) has proven a powerful method in studying statistical systems which suffer from the sign problem. Here we show that it can also be applied to an important problem concerning why we live in four-dimensional spacetime. Our target system is the type IIB matrix model, which is conjectured to be a nonperturbative definition of type IIB superstring theory in ten dimensions. The fermion determinant of the model becomes complex upon Euclideanization, which causes a severe sign problem in its Monte Carlo studies. It is speculated that the phase of the fermion determinant actually induces the spontaneous breaking of the SO(10) rotational symmetry, which has direct consequences on the aforementioned question. In this paper, we apply the CLM to the 6D version of the type IIB matrix model and show clear evidence that the SO(6) symmetry is broken down to SO(3). Our results are consistent with those obtained previously by the Gaussian expansion method.


abstract
In recent years the complex Langevin method (CLM) has proven a powerful method in studying statistical systems which suffer from the sign problem. Here we show that it can also be applied to an important problem concerning why we live in four-dimensional spacetime. Our target system is the type IIB matrix model, which is conjectured to be a nonperturbative definition of type IIB superstring theory in ten dimensions. The fermion determinant of the model becomes complex upon Euclideanization, which causes a severe sign problem in its Monte Carlo studies. It is speculated that the phase of the fermion determinant actually induces the spontaneous breaking of the SO(10) rotational symmetry, which has direct consequences on the aforementioned question. In this paper, we apply the CLM to the 6D version of the type IIB matrix model and show clear evidence that the SO(6) symmetry is broken down to SO (3). Our results are consistent with those obtained previously by the Gaussian expansion method.

Introduction
Monte Carlo methods have been playing a crucial role in nonperturbative studies of quantum field theories and statistical systems relevant to particle, nuclear and condensed matter physics. However, in many interesting cases, it happens that such methods cannot be applied straightforwardly because the effective "Boltzmann weight" can become negative or even complex. A brute-force method would be to use the absolute value of the weight in generating configurations and to take into account the sign or the phase in calculating the expectation values. This reweighting method indeed works for small systems, but the computational cost grows exponentially with the system size due to huge cancellations among configurations, which is commonly referred to as the sign problem. This problem occurs, for instance, in studying finite density systems, including that of QCD, the real time dynamics of quantum systems, supersymmetric theories and strongly correlated electron systems.
In recent years there has been major progress in evading the sign problem by complexifying the dynamical variables, which are supposed to be real in the original system. One such approach is the generalized Lefschetz-thimble method [1][2][3][4], which amounts to deforming the integration contour in such a way that the sign problem becomes mild enough to be handled by the reweighting method. Another approach is the complex Langevin method (CLM) [5,6], which attempts to define a stochastic process for the complexified variables so that the expectation values with respect to this process are equal to the expectation values defined in the original system extending the idea of stochastic quantization [7]. In both approaches, holomorphy plays a crucial role. The advantage of the CLM as compared with the other is that it is computationally less costly, which enables its application to much larger system size. The disadvantage, on the other hand, is that the equivalence to the original system does not always hold. This has been a serious issue for a long time because one can obtain explicit numerical results which are simply wrong without even noticing it. The progress in this direction was made by clarifying the conditions for the equivalence [8][9][10][11][12][13] and by inventing new techniques that made it possible to meet these conditions for a larger space of parameters [14][15][16][17][18][19][20][21].
Here we aim at applying the CLM to the so-called type IIB matrix model [22], which is conjectured to provide a nonperturbative definition of superstring theory. In its Euclidean version, this model suffers from a severe sign problem due to the complex Pfaffian, which appears after integrating out the fermionic matrices. The phase of the Pfaffian fluctuates less violently for configurations with lower dimensionality [23], and this effect is speculated to cause the spontaneous breaking of the SO(10) rotational symmetry for ten bosonic matrices that represent the spacetime [24,25]. While the original expectation was that the SO(10) symmetry is spontaneously broken to SO(4) in order to account for the appearance of four-dimensional spacetime [26,27], explicit calculations based on the Gaussian expansion method (GEM) suggested that it is broken down to SO(3) instead [28]. In order to address this issue by first principle calculations, we clearly need to overcome the sign problem, which plays a crucial role in the phenomenon itself.
A first step in that direction has been taken by two of the authors (Y. I. and J. N.) in ref. [18], where a toy model [29] in which SO(4) symmetry is expected to break spontaneously due to the complex fermion determinant has been studied by the CLM. In particular, a new technique has been proposed to avoid the so-called singular-drift problem [10] caused by the eigenvalues of the Dirac operator close to zero [30]. The idea is to avoid this problem by deforming the Dirac operator with a kind of mass term. This makes it possible to satisfy the condition for justifying the method as one can confirm by probing the histogram of the drift term appearing in the complex Langevin equation [11]. After extrapolating the deformation parameter to zero, the CLM was able to reproduce the results of the GEM [31] including the pattern of the spontaneous symmetry breaking (SSB), which was not possible [32,33] by a Monte Carlo method [25] based on reweighting.
In this paper, we extend this work to the 6D version of the type IIB matrix model, which can be obtained by dimensionally reducing 6D super Yang-Mills theory to a point in the same way as one obtains the type IIB matrix model by dimensionally reducing 10D super Yang-Mills theory. The 6D version also suffers from the sign problem due to a complex determinant, which appears after integrating out the fermionic matrices, and its phase is speculated to cause the SSB of SO(6) rotational symmetry to SO(3) according to the GEM [34]. We show that the CLM indeed enables us to address this issue from first principles with the aid of the deformation technique, and our results turn out to be consistent with those of the GEM. Here again, the success of the CLM is remarkable compared with the marginally successful results of the reweighting-based method [35]. This gives us a big hope that the original type IIB matrix model can also be studied by the CLM, and that one can determine the pattern of the SSB, which seemed to be extremely hard in the reweighting-based method [36].
The rest of this paper is organized as follows. In section 2 we briefly review the known results for dimensionally reduced super Yang-Mills models. In section 3 we discuss how we apply the CLM to the 6D version of the type IIB matrix model and present the results in section 4. Section 5 is devoted to a summary and discussions. In appendix A we provide some details of our complex Langevin simulation.

Brief review of dimensionally reduced super Yang-Mills models
As is well known, one can define N = 1 pure super SU(N) Yang-Mills theories in D = 3, 4, 6, 10 dimensions. 6 By dimensionally reducing these theories to a point, one obtains matrix models with D bosonic matrices A µ and their superpartners ψ α , which are called dimensionally reduced super Yang-Mills models in the literature. In particular, the D = 10 case corresponds to the type IIB matrix model [22], which has been proposed as a nonperturbative definition of superstring theory in the same sense as the lattice formulation provides a nonperturbative definition of quantum field theories. Of particular interest is the fact that the spacetime is represented by the eigenvalue distribution of the bosonic matrices A µ in this model [26], and hence even the spacetime dimensionality should be determined dynamically.
In what follows, we consider the Euclidean version of the dimensionally reduced models, in which the spacetime indices are contracted using the Kronecker delta instead of the Minkowski metric. 7 The D = 3 model is ill-defined since the partition function is divergent [41,42]. The D = 4 model has a real nonnegative fermion determinant, and Monte Carlo simulation suggested that the SO(4) rotational symmetry is not spontaneously broken [43]. The D = 6 and 10 models have a complex fermion determinant/Pfaffian, whose phase is expected to play a crucial role in the SSB of SO(D) [23-25, 35, 36].
In this paper we therefore focus on the D = 6 model defined by the partition function as a simplified model of the type IIB matrix model. The bosonic part S b and the fermionic part S f of the action are given, respectively, as 2) The bosonic matrices A µ (µ = 1, · · · , 6) are traceless Hermitian N × N matrices, while the fermionic matrices ψ α (α = 1, · · · , 4) are traceless N × N matrices with Grassmann entries. The action is invariant under SO(6) transformations, under which A µ transforms as a vector, while ψ α transforms as a Weyl spinor. The 4 × 4 gamma matrices Γ µ obtained after the Weyl projection are given, for instance, as in terms of the Pauli matrices σ i (i = 1, 2, 3).
Integrating out the fermionic matrices ψ α , we obtain acting on the linear space of traceless complex N ×N matrices Ψ α . The determinant det M takes complex values in general, and we define its phase Γ by det M = | det M| e iΓ . If one 7 The Lorentzian version of the dimensionally reduced models is also found to have interesting dynamical properties [37][38][39][40]. The relationship between the two versions is not clear, though. omits the phase Γ, Monte Carlo studies become straightforward, and it is found that the SSB of SO(6) does not occur [35]. (Similar results are obtained in the 10D case as well [36].) Let us here define "d-dimensional configurations" by those configurations which can be transformed into a configuration with A d+1 = · · · = A 6 = 0 by an appropriate SO(6) transformation. One can then prove the following properties of the determinant [23]. For d = 5 configurations, det M is real. For d = 4 or d = 3 configurations, we obtain which implies that the phase Γ becomes more stationary for d = 3 than for d = 4. For d = 2 configurations, we have det M = 0 and the phase Γ becomes ill-defined. While these properties suggest that the SO(6) symmetry may be broken spontaneously down to SO (3), identifying the actual pattern of the SSB is a nontrivial issue, which should be addressed by taking into account the competition between the effect of the phase Γ discussed above and the entropic effect that favors configurations with higher dimensionality.
In order to address this issue, the free energy of the SO(d) symmetric vacuum (d = 2, 3, 4, 5) was obtained by the GEM up to the fifth order, and it was found that the SO(3) vacuum has the lowest free energy [34]. This implies that the SO(6) symmetry is spontaneously broken to SO(3). In the SO(d) vacuum, the extent of spacetime λ µ = 1 N tr(A µ ) 2 in each direction has d large values and (6 − d) small values, which implies that the dynamically generated spacetime has d extended directions and (6 − d) shrunken directions. This quantity λ µ was calculated up to the fifth order in the GEM as well, and the result for d = 3 turned out to be λ µ ≃ 1.7 for the three extended directions, 0.2 for the three shrunken directions. (2.8) Let us also mention that the free energy for the SO(2) vacuum turned out to be substantially higher than that for the vacua with higher dimensionality [34]. This is consistent with the fact that d = 2 configurations are suppressed by the fermion determinant. (See also ref. [28] for similar results in the 10D case.) 3 Applying the CLM to the 6D type IIB matrix model In this section, we apply the CLM to the 6D version of the type IIB matrix model (2.1) following an analogous study on a simplified model [18]. In particular, we discuss how we probe the SSB of SO(6) symmetry and explain the important techniques such as the gauge cooling and the deformation, which are crucial in making the CLM work.

complex Langevin equation
where the effective action S = S b − log det M for the bosonic matrices is complex in general due to the complex fermion determinant. In the CLM, one considers a fictitious time evolution of the bosonic matrices defined by the complex Langevin equation [5,6] where η µ (t) are traceless Hermitian matrices whose elements are random variables obeying the Gaussian distribution ∝ exp − 1 4 tr {η µ (t)} 2 dt . The first term on the right-hand side of eq. (3.2) is called the drift term, which is given explicitly as where Tr represents the trace for a 4( is not Hermitian as a result of the fact that the fermion determinant is complex. Therefore, when we consider the fictitious time evolution based on the complex Langevin equation where t 0 is the time required for thermalization, and T should be large enough to achieve good statistics. Let us recall here that the observable O[A] and the drift term (3.3) are originally defined for Hermitian A µ . In the above procedure, their definitions need to be extended to complexified A µ (t) by analytic continuation, which plays a crucial role in the argument for justification [8,9,11]. When we solve the complex Langevin equation (3.2) numerically, we have to discretize the fictitious time t. The discretized version of (3.2) reads where the probability distribution of η µ (t) is given by exp − 1

how to probe the SSB
As is commonly done in probing the SSB, we introduce a term in the action with the order 0 < m 1 ≤ · · · ≤ m 6 , which breaks the SO(6) symmetry explicitly. After taking the thermodynamic limit, which amounts to taking the large-N limit in our case, we send the coefficient ε in (3.6) to zero.
As an order parameter for the SSB, we consider 8 where no sum over µ is taken. Note that λ µ calculated for a configuration A µ generated by the CLM is not necessary real because A µ is no longer Hermitian. However, it becomes real after taking an ensemble average due to the symmetry property of the drift term (3.3) under A i → (A i ) † (i = 1, · · · , 5) and A 6 → −(A 6 ) † . For this reason, we take the real part of λ µ when we define the expectation values λ µ . Due to the chosen ordering of m µ , we have an inequality λ 1 ≥ · · · ≥ λ 6 for positive ε. When we take the large-N limit followed by the ε → 0 limit, it can happen that the expectation values λ µ are not equal. In that case we conclude that the SO(6) symmetry is spontaneously broken. Throughout this paper, we use m µ = (0.5, 0.5, 1, 2, 4, 8) , which retains SO(2) ⊂ SO(6) instead of breaking the SO(6) symmetry completely. The reason for making this compromise is that having m 1 = m 2 necessarily results in a wider spectrum of m µ , which makes the ε → 0 extrapolation more subtle. This is expected not to harm anything since it is unlikely that the SO(6) symmetry is broken completely according to the discussion below (2.7).

gauge cooling for the excursion problem
In order to justify the CLM, the probability distribution of the magnitude of the drift term 8 Unlike the previous studies [35,36] using a method based on reweighting, we cannot use the eigenvalues of the tensor T µν = 1 N tr(A µ A ν ) as an order parameter since they are not single-valued with respect to complexified A µ and hence the relationship (3.4) does not hold. should fall off exponentially or faster [11]. There are two sources for the violation of this property. One is the "excursion problem", which occurs when A µ develops a large anti-Hermitian part. The other is the "singular-drift problem", which occurs because of the appearance of M −1 in (3.3) when some eigenvalues of M come close to zero frequently.
In order to avoid the excursion problem, we use a technique called gauge cooling [14], which keeps A µ as close to Hermitian matrices as possible. This amounts to minimizing the "Hermiticity norm" defined by after each step of solving the discretized Langevin equation (3.5). Here, G is the gradient of N H with respect to the SL(N, C) transformation, which is derived in ref. [18]. We choose the real positive parameter α in such a way that the Hermiticity norm N H is minimized. In refs. [11,15], it was proven that adding the gauge cooling procedure in the CLM does not affect the argument for its justification.

deformation for the singular-drift problem
In order to avoid the singular-drift problem, we deform the fermionic action by adding where m f is the deformation parameter. This technique was proposed originally in ref. [18], where the singular-drift problem was indeed overcome in an SO(4) symmetric matrix model with a complex fermion determinant. The fermionic mass term (3.12) modifies the linear transformation M in (2.6) to M given by using Γ 6 = 1 ⊗ 1, which implies that the eigenvalue distribution of the matrix M is shifted in the direction of the real axis for the same A µ . We will see that this effect enables us to avoid the eigenvalues coming close to zero. Note that we cannot construct an SO(6) symmetric mass term for ψ α since ψ α transforms as a Weyl spinor under an SO(6) transformation. The term (3.12) we have chosen is a kind of mass term, which breaks the SO(6) symmetry minimally to SO(5). We first investigate whether this SO(5) symmetry of the deformed model is spontaneously broken for various values of m f , and then discuss what happens as m f is decreased.
Let us also note that the m f → ∞ limit of the deformed model is nothing but the SO(6) symmetric bosonic model since fermionic matrices ψ α obviously decouple in that limit. Thus, the deformation (3.12) may be regarded as an interpolation between the dimensionally reduced super Yang-Mills model and its bosonic counterpart, in which the fermionic matrices are omitted. It is known that the SO(D) symmetry of the D-dimensional version of the bosonic type IIB matrix model is not spontaneously broken [44]. We use this fact to test the validity of our method.
As the deformation paramter m f is decreased, the effects of fermionic matrices are gradually turned on. From the trend for decreasing m f , we try to draw some conclusions on the undeformed model, assuming that nothing dramatic occurs in the vicinity of m f = 0. This strategy has been quite successful in the simplified model [18], where the spontaneous breaking of SO(4) symmetry to SO(2) was confirmed as predicted by the GEM [31].
To summarize, the model we investigate by the CLM is defined by the partition function  While the singular-drift problem is cured for large enough m f , the problem occurs for moderate m f with ε smaller than some value depending on m f . Therefore, when we make an ε → 0 extrapolation, we have to choose carefully the data points that do not suffer from this problem. For that purpose, we probe the probability distribution p(u) of the magnitude of the drift term u defined by eq. (3.9). In Fig. 1 we show the log-log plot of p(u) for N = 24 with m f = 0.65, 0.8, 0.9. We find at m f = 0.65 that p(u) falls off exponentially or faster for ε ≥ 0.15, while a clear power-law tail develops for ε ≤ 0.125. Therefore, we can trust only the results for ε ≥ 0.15 according to the criterion of ref. [11]. Similarly, at m f = 0.80, we find that p(u) falls off exponentially or faster for ε ≥ 0.05. At m f = 0.90, we observe such a behavior for all values of ε investigated.
In order to understand how the singular-drift problem is avoided by the deformation (3.12), we show in Fig. 2 the scatter plot of the eigenvalues of the 4(N 2 − 1) × 4(N 2 − 1) matrix M obtained with a thermalized configuration for N = 24 and m f = 0.65 with ε = 0.1 (Left) and ε = 0.25 (Right). Note that the eigenvalue distribution is shifted in the direction of the real axis compared with that for m f = 0 as one can deduce from (3.13). From this figure, we find that eigenvalues close to zero appear for ε = 0.1, but not for ε = 0.25.

Results
In this section we present our results obtained in the way described in the previous section. Let us recall that the deformation (3.12) breaks the SO(6) symmetry to SO (5). In order to probe the SSB of SO(5) that remains for m f = 0, we need to take the N → ∞ limit first, and then the ε → 0 limit. Finally, in order to obtain results for the undeformed model, we need to extrapolate m f to zero.
First we discuss how we take the N → ∞ limit. In Fig. 3 we plot the expectation values λ µ for ε = 0.25 and m f = 0.65 with N = 24, 32, 40, 48 against 1/N. For λ 1 and λ 2 , we plot the average to increase the statistics since we know theoretically that they should be equal due to the chosen parameters (3.8) for the symmetry breaking term. We find that our data can be fitted nicely to straight lines, which implies that our data follow the large-N asymptotic behavior a + b/N. Similar behaviors were observed for other (ε, m f ) as well. Thus we can make extrapolations to N = ∞ reliably. Let us denote the extrapolated values obtained for each (ε, m f ) as λ µ ε,m f . The straight lines represent fits to a + b/N, which enable us to make reliable extrapolations to N = ∞.
Next we make an extrapolation to ε = 0. For that purpose, it is convenient to consider the ratio [18] ρ µ (ε, m f ) = λ µ ε,m f SO(5) to SO (4), whereas at m f = 1.4, we find that the SO(5) symmetry of the deformed model is not broken spontaneously. Let us recall that in the large-m f limit, the deformed model reduces to the SO(6) symmetric bosonic model, in which the SSB of SO(6) is known not to occur [44]. Indeed our results for m f = 1000 in Fig. 4 (Bottom-Right) show that the fitting curves for (ρ 1 + ρ 2 )/2 and ρ µ (µ = 3, · · · , 6) all merge at ε = 0 as expected. This confirms that the SSB observed for smaller m f is a physical effect, which cannot be attributed to some artifacts in the ε → 0 extrapolations.
In Fig. 5 (Left), we plot the ε → 0 extrapolated values of ρ µ (ε, m f ) against m f 2 for m f = 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4. We find that an SO(3) vacuum is realized for m f 0.9, while an SO(4) vacuum is realized for 1.0 m f 1.3. Judging from the pattern of the SSB with decreasing m f , it is reasonable to consider that the symmetry that survives for m f < 0.65 is at most SO(3). Note also that an SO(2) vacuum is suppressed by the fermion determinant as we discussed below (2.7). Combining this argument with our results, we conclude that the SO(6) rotational symmetry is spontaneously broken to SO(3) in the undeformed model corresponding to m f = 0, which agrees with the prediction from the GEM. Note that the results (2.8) obtained by the GEM [34] can be rephrased in terms of the ratios ρ µ as These values are represented by the filled squares at m f = 0 in Fig. 5. As m f decreases, the data points for (ρ 1 + ρ 2 )/2 and ρ 3 are coming closer to 0.3, whereas the data points for ρ 4 , ρ 5 and ρ 6 are coming closer to 0.035. For the sake of clearer comparison, we plot in Fig. 5 (Right), the averages (ρ 1 + ρ 2 + ρ 3 )/3 and (ρ 4 + ρ 5 + ρ 6 )/3 for m f ≤ 0.9 corresponding to the SO(3) symmetric phase. The asymptotic behavior of these quantities for m f → 0 is expected to be a power series with respect to m f 2 due to the symmetry under m f → −m f . We therefore fit our data points to the form a + bm f 2 + cm f 4 with the fitting range 0.65 ≤ m f ≤ 0.9, and the extrapolation to m f = 0 yields which are close to the values (4.2) predicted by the GEM. In fact, we can make reasonable fits passing through the points at m f = 0 predicted from the GEM as shown by the dashed lines. Note, however, that the GEM involves a truncation, which necessarily yields certain systematic errors. Hence, precise agreement is not anticipated.

Summary and discussions
In this paper we have discussed the SSB of rotational symmetry conjectured to occur in dimensionally reduced super Yang-Mills models for D = 6 and D = 10. In particular, the D = 10 case is relevant to the dynamical generation of four-dimensional space-time in a nonperturbative formulation of superstring theory in ten dimensions. It is known that the phase of the complex fermion determinant should play a crucial role, which implies that a first principle investigation of this issue based on Monte Carlo methods necessarily faces a severe sign problem.
Extending the previous work [18] on a simplified model, we have investigated the D = 6 case using the CLM. In particular, we use the deformation technique to overcome the singular-drift problem, which occurs in many interesting systems with a complex fermion determinant including finite density QCD at low temperature. In the data analysis, it is important to probe the probability distribution of the drift term so that we can determine the parameter region in which the CLM is valid. After taking the large-N limit, we are able to observe that the remaining SO(5) symmetry of the deformed model is spontaneously broken to SO(4) or to SO(3) as the deformation parameter m f is decreased. Combining this result with the argument that an SO(2) vacuum is suppressed by the fermion determinant, we conclude that an SO(3) vacuum is chosen in the undeformed model, which is consistent with the prediction of the GEM. Extrapolations to m f = 0 for the extent of spacetime in each direction also give results consistent with the values predicted by the GEM. These results go far beyond those obtained by a reweighting-based method [35].
An application of the same method to the D = 10 case is ongoing [47]. In this case the GEM results [28] suggest that the SO(10) symmetry is spontaneously broken to SO(3) as opposed to the original expectation that the SO(4) symmetry survives. Whether we can observe an SO(3) vacuum for nonzero m f already gives us an important clue on this issue. The fact that the CLM with the deformation technique turned out to be useful in investigating such a physically interesting issue that has been hard to address due to the severe sign problem encourages us to apply it also to other interesting complex action systems such as finite density QCD [48].

A Details of the complex Langevin simulation
In this section we provide some details of our complex Langevin simulation.
First we discuss the idea of adaptive step-size [46] used in solving the discretized complex Langevin equation (3.5) numerically. The point here is that the drift term in (3.5) can become large occasionally, and the associated discretization artifacts can make the simulation imprecise or even unstable in the worse case. In order to avoid these problems, we probe the magnitude of the drift term (3.9) at each step, and when it gets larger than a certain threshold value u 0 , we decrease the step-size as ∆t = (∆t) 0 × u 0 /u (for u ≥ u 0 ). A fixed step-size (∆t) 0 = 10 −5 is used during the thermalization process, and the threshold value u 0 is determined by taking an average of (3.9) during that process.
Next we discuss how we estimate the second term in the drift term (3.3), which is the most time-consuming part of our calculation. Since the matrix M has the size 4(N 2 − 1) × 4(N 2 − 1), direct calculation of M −1 would require O(N 6 ) arithmetic operations. We can reduce this cost to O(N 3 ) by using the so-called noisy estimator. 10 The idea is based on the identity where the average · χ is taken with respect to the Gaussian noise χ, which represents a 4(N 2 − 1)-dimensional vector whose elements are complex Gaussian variables normalized as χ * k χ l = δ kl . In practice, we generate the Gaussian noise only once and estimate the trace using it instead of taking the average on the right hand side of (A.1). The use of this approximation does not yield any systematic errors in the CLM in the ∆t → 0 limit since the associated Fokker-Planck equation remains the same [49].
The quantity ζ = M −1 χ in (A.1) can be calculated by solving the linear equation where M † M is a Hermitian positive semi-definite matrix, using the conjugate gradient method. This method is an iterative one, in which one acts M and M † on a 4(N 2 − 1)dimensional vector many times. If one does this directly using an explicit representation of M, it requires O(N 4 ) arithmetic operations. In fact, one can perform this calculation using (2.6) for M and for M † , which requires only O(N 3 ) arithmetic operations. In the parameter region in which the CLM is valid, the number of iterations required for convergence of the conjugate gradient method is almost independent of N, and it is typically of the order of 100. Thus the computational cost of our simulation grows only as O(N 3 ) with the matrix size N.