Complex Langevin analysis of the spontaneous breaking of 10D rotational symmetry in the Euclidean IKKT matrix model

The IKKT matrix model is a promising candidate for a nonperturbative formulation of superstring theory, in which spacetime is conjectured to emerge dynamically from the microscopic matrix degrees of freedom in the large-$N$ limit. Indeed in the Lorentzian version, Monte Carlo studies suggested the emergence of (3+1)-dimensional expanding space-time. Here we study the Euclidean version instead, and investigate an alternative scenario for dynamical compactification of extra dimensions via the spontaneous symmetry breaking (SSB) of 10D rotational symmetry. We perform numerical simulations based on the complex Langevin method (CLM) in order to avoid a severe sign problem. Furthermore, in order to avoid the singular-drift problem in the CLM, we deform the model and determine the SSB pattern as we vary the deformation parameter. From these results, we conclude that the original model has an SO(3) symmetric vacuum, which is consistent with previous results obtained by the Gaussian expansion method (GEM). We also apply the GEM to the deformed matrix model and find consistency with the results obtained by the CLM.

space-time. Here we study the Euclidean version instead, and investigate an alternative scenario for dynamical compactification of extra dimensions via the spontaneous symmetry breaking (SSB) of 10D rotational symmetry. We perform numerical simulations based on the complex Langevin method (CLM) in order to avoid a severe sign problem. Furthermore, in order to avoid the singular-drift problem in the CLM, we deform the model and determine the SSB pattern as we vary the deformation parameter. From these results, we conclude that the original model has an SO(3) symmetric vacuum, which is consistent with previous results obtained by the Gaussian expansion method (GEM). We also apply the GEM to the deformed matrix model and find consistency with the results obtained by the CLM.

Introduction
Superstring theory has been studied intensively as a unified theory that includes quantum gravity. The theory is defined in ten spacetime dimensions and the connection to the real world, where only four dimensions are macroscopic, is realized via compactification of the extra dimensions. How this can actually occur has been investigated perturbatively by using D-brane configurations as a background, leading to tremendously many vacua giving rise to the so-called string landscape. Clearly, it is important to see if this picture remains valid when the issue is addressed in a fully nonperturbative manner.
The IKKT matrix model [1], also known as the type IIB matrix model, was proposed as a nonperturbative formulation of superstring theory. Formally, the action of the model can be obtained by dimensionally reducing the action of 10D N = 1 SU(N) super Yang-Mills (SYM) theory to 0D. In this model, spacetime emerges dynamically from the eigenvalues of the ten bosonic matrices in the large-N limit [2], which enables the scenario for dynamical compactification as a purely nonperturbative effect.
Evidence supporting such a scenario has been provided by simulating the Lorentzian version of the IKKT matrix model [3][4][5][6][7][8][9][10][11]. In these simulations, continuous time emerges dynamically and three-dimensional space undergoes expansion after a critical time, with the six extra dimensions remaining small. This is highly nontrivial since time is given by the ordered eigenvalues of the temporal matrix A 0 , and in the SU(N) basis used to diagonalize A 0 , the dominant configurations of the spatial matrices A i have a band diagonal structure, from which one can read off the time evolution of space. By taking the large-N limit, the eigenvalue distribution of A 0 extends in physical units and a sensible continuum limit can be defined, a fact that emerges also from the dynamics of the model. The results in Ref. [5] suggest that the expansion is exponential at early times, which turns into a power law at later times [6], providing evidence that a realistic cosmological scenario may also arise dynamically from this model.
Monte Carlo simulations of the Lorentzian model are hindered by a severe sign problem coming from the phase factor e iS b , where S b is the bosonic part of the action. In the early work [3][4][5][6][7][8][9][10], this problem was avoided by first integrating out the scale factor of the bosonic matrices, yielding a function of S b which is sharply peaked near the origin. Then, by approximating this function by a sharply peaked Gaussian function, the sign problem was avoided. This, however, leads to singular spatial configurations, showing that it is highly nontrivial to obtain 3D expanding space with a smooth structure [10]. Recent work, which avoids this approximation but confronts the sign problem, provided evidence that the 3D expanding space can have a smooth structure in the large-N limit [11]. This work also shows that the complex Langevin method (CLM) [12,13] can be used successfully in circumventing the sign problem in this model.
The Euclidean version of the IKKT matrix model and related models, on the other hand, have been studied numerically [14][15][16][17][18][19][20][21] since long time before numerical studies of the Lorentzian version were started. This is because it can be thought of as being the direct analog of lattice QCD for superstring theory, and moreover it was shown [22,23] to have a finite partition function without introducing infrared cutoffs unlike the Lorentzian version. However, the simulations are hindered by a severe sign problem here as well because the Pfaffian obtained after integrating out the fermionic matrices is complex in general.
In the Euclidean model, the scenario for dynamical compactification of extra dimensions is expected to be realized via spontaneous symmetry breaking (SSB) of the SO(10) rotational symmetry due to the wild fluctuations of the phase of the Pfaffian for SO(d) symmetric configurations with larger d [24,25]. A strong evidence for the SSB has been provided by studying the model using the Gaussian expansion method (GEM). This method was applied to the Euclidean IKKT model and related matrix models [26][27][28][29][30] realizing SSB to lower dimensional spaces. In particular, it was shown [30] in the IKKT matrix model that the SO(3) symmetric vacuum has the lowest free energy, which implies SSB to SO (3). The ratio between the three extended directions and the seven shrunken directions was also calculated and it was found to be finite. This SSB can be naturally attributed to the effect of the phase of the Pfaffian since the phase-quenched simulations of the model showed no SSB [16]. A direct confirmation of this scenario by numerical simulation, however, requires some new ideas to overcome the sign problem. In a series of work [18][19][20], the Euclidean version of the IKKT matrix model was studied by using a density of states based method [18,[31][32][33][34], which was successful in that it made it possible to obtain the extent of space in each direction for a given pattern of SSB although it was not powerful enough to determine the SSB pattern itself.
In fact, the sign problem occurs in Monte Carlo simulation of various interesting systems such as finite density QCD, supersymmetric theories, strongly correlated electron systems and real-time quantum field theories. If one uses the reweighting method to simulate the system, the computational cost increases exponentially with the system size. In recent years there has been major progress in evading the sign problem by complexifying the dynamical degrees of freedom of the system under study. One of such methods is the generalized Lefschetz-thimble method [35][36][37][38][39][40], 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 CLM [12,13], which extends the idea of stochastic quantization [41] and defines 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. The use of the CLM allows one to study large systems, but it has the problem of not always yielding correct results. (See, for instance, Refs. [42,43] for early work.) Recently, the conditions for the correct convergence were clarified [44][45][46][47][48][49] and various new techniques have been proposed to meet these conditions for a large space of parameters [50][51][52][53][54][55]. Thanks to these developments, the CLM has been applied successfully to many systems in lattice quantum field theory [49,[56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] and matrix models [11,21,49,53,[72][73][74][75][76] with complex actions. (For a recent review on the CLM and related methods, see Ref. [77].) In Ref. [21], we applied the CLM to the 6D version of the Euclidean IKKT matrix model, and obtained results consistent with the ones obtained by the GEM [29], i.e., SSB to SO(3). We used a technique [53] to avoid the singular-drift problem [46] caused by the eigenvalues of the Dirac operator that accumulate near zero [72]. The model is deformed by introducing a parameter m f , which corresponds to adding a mass-like term in the Dirac operator. Thus the condition [47] that ensures the correctness of the CLM can be met, and one can make an extrapolation m f → 0 to obtain the results for the original model.
In this paper we apply the same method to the original Euclidean IKKT matrix model. The simulations are more challenging than those of the 6D version because the number of the fermionic matrices increases by a factor of four. We find that SSB to lower dimensional space occurs as the deformation parameter m f is reduced and that our data are consistent with an SO(3) symmetric vacuum appearing at m f = 0 as predicted by the GEM for the undeformed model [30]. These results are in sharp contrast to the ones obtained in Ref. [20] using the density of states based method, where calculations had to be performed for each vacuum with SO(d) symmetry. We also find at larger N that the singular-drift problem becomes milder, which enables the use of a smaller deformation parameter in the simulations. This gives us a hope that a more complete understanding of the Euclidean IKKT matrix model may be possible by extending this work to larger N.
Our calculations may still suffer from finite-N effects and the results can be sensitive to systematic errors introduced by extrapolations including that for large N. We therefore perform a consistency check by applying the GEM to the deformed IKKT matrix model. We calculate the free energy up to three loops and obtain physical solutions for the selfconsistency equations with the SO(d) (d = 6, 7) ansatzes. We find that the SO(6) symmetric vacuum has smaller free energy as m f is decreased. We also calculate the extent of space and find consistency with the results obtained by the CLM.
The rest of this paper is organized as follows. In Section 2, we give a brief review of the Lorentzian and Euclidean versions of the IKKT matrix model. In Section 3, we apply the CLM to the Euclidean IKKT matrix model with a mass deformation of the fermionic action and present our results. In Section 4, we apply the GEM to the deformed Euclidean IKKT matrix model and compare the results with those of the CLM. Section 5 is devoted to a summary and discussions. In Appendix A, we show some results suggesting that the singular-drift problem vanishes at large N.

the Lorentzian version
In this section we review the Lorentzian version of the IKKT matrix model, in which the indices are contracted using the Minkowski metric η µν = diag(−1, 1, 1, . . . , 1). The model is invariant under SO(9, 1) Lorentz transformations, which act on the vectors A µ and the Majorana-Weyl spinors ψ α . The model also possesses the SU(N) symmetry which is inherited from the gauge invariance of the 10D N = 1 SYM action after reduction to 0D. The supersymmetry of the SYM theory, on the other hand, enhances to an N = 2 supersymmetry in the IKKT matrix model [1]. This allows us to interpret the eigenvalues of the matrices A µ as the N points in the target spacetime [2], which are expected to represent the continuum spacetime in the large-N limit 1 . The partition function is given by [3] Z = dA dψ e iS = dA e iS b Pf M , (2.5) where the Pfaffian Pf M comes from integrating out the fermionic matrices ψ α . The 16(N 2 − 1) × 16(N 2 − 1) anti-symmetric matrix M is defined by its action on the linear space of traceless complex N × N matrices. In fact, it turns out that the Pfaffian Pf M is real 2 in the present Lorentzian model. The bosonic action S b can be written as are Hermitian matrices and i, j = 1, . . . , 9 are spatial indices. Since the partition function (2.5) is divergent as it is, one has to introduce cutoffs in the temporal and spatial directions [3]. Using Eq. (2.4), it is possible to choose a gauge that diagonalizes A 0 as In this gauge, the spatial matrices A i turn out to have a band-diagonal structure and for an appropriate integer n, the n × n submatricesĀ i can effectively represent space at time t ν defined by where ν = 0, . . . , N − n [3][4][5][6][7][8][9]. Time emerges dynamically and it is a nontrivial dynamical question whether this leads to a continuum time in the large-N limit. In Ref. [5], it was shown to be possible to take a continuum limit in the large-N limit such that the "volume" ∆ and the lattice spacing ǫ in time can be tuned to go to ∞ and 0, respectively, keeping the product ǫ ∆ fixed. Using this definition of time, it was found [3][4][5][6][7][8][9] that there exists a critical time t c , after which three spatial directions undergo rapid expansion, whereas the other directions remain small. This happens due to spontaneous breaking of the SO(9) rotational symmetry down to SO (3). In order to see this, we define the "moment of inertia tensor" where the trace here is over the I, J indices in (2.9), and obtain its nine eigenvalues λ i (t) with the ordering λ 1 (t) > λ 2 (t) > . . . > λ 9 (t). When the expectation values λ i (t) for i = 1, 2, 3 are equal but larger than λ i (t) for i = 4, . . . , 9 in the large-N limit, we conclude that SSB to SO(3) occurs. Using simplified models that describe the qualitative behavior of the IKKT model at early and late times respectively, it was shown that at early times the large eigenvalues grow exponentially λ i ∼ e Λt with t [5], whereas at late times the expansion turns into a power law λ i ∼ t 1/2 [6]. This gives us a hope that the IKKT model has the dynamics that contain a realistic cosmology with an early time inflationary expansion and a late time FRW expansion in the radiation dominated era. The structure of space was examined recently in Ref. [10], and it was found to be dominated by rather singular spatial configurations, whose (3+1)D expanding behavior is due to submatrices that are close to the Pauli matrices. Namely, the radial distribution of spacetime points is such that two points are located very far, whereas the rest accumulate near the origin. The reason for the domination of such configurations was attributed to an approximation used in order to avoid the sign problem in the simulation. As a result of this approximation, one effectively simulates a model (2.5) with e iS b replaced by e S b .
Having the factor e S b in the partition function makes configurations with the Paulimatrix structure dominant. This can be understood by looking at Eq. (2.7). Note that the first term favors configurations such that the A i commute with the A 0 , whereas the second term favors configurations such that the A i are maximally noncommuting. The balance of these two terms gives the band-diagonal structure important in defining (2.9). In Ref. [10] it was shown that configurations with the Pauli-matrix structure maximize the second term, subject to the constraint 1 N tr (A i ) 2 = 1 coming from the spatial cutoff introduced in the model to make the partition function finite [3].
The approximation, however, may miss the important contribution of the configurations that extremize S b instead of maximizing it. This is suggested indeed in Ref. [11], where the D = 6 bosonic IKKT matrix model was studied numerically using the CLM in order to avoid the sign problem without approximations. The model was generalized by two parameters s and k, which correspond to Wick rotations on the worldsheet and on the target space, respectively. The Lorentzian model corresponds to (s, k) = (0, 0) and the Euclidean model studied in this paper corresponds to (s, k) = (1, 1). The results for (s, k) = (−1, 0), which correspond to replacing e iS b by e S b , are consistent with the results obtained by the previous simulations [3][4][5][6][7][8][9] of the Lorentzian model using the approximation, and show that the dominant configurations are singular in that the 3D expanding space has the Pauli-matrix structure. The generalized model was also studied in the vicinity of s = 0 and it was shown for a range of parameters (s, k) that it exhibits a (3+1)D expanding behavior, while the dominant configurations depart from the Pauli-matrix structure and the spacetime points are distributed more smoothly than for (s, k) = (−1, 0).
Since an infinite number of (3+1)D expanding classical solutions without the Paulimatrix structure are known to exist [78][79][80][81], it is possible to imagine that the spacetime structure becomes smooth without loosing the (3+1)D expanding behavior if one can approach (s, k) = (0, 0) in the large-N limit. Furthermore, the simulation supported a speculation that some classical solution dominates at late times due to the expansion of space. This is important because it shows the possibility to understand the late-time behavior of the model by finding classical solutions that contain a realistic cosmology [82][83][84][85][86][87][88][89][90][91][92]. It is also expected that the solution that dominates at late times can accommodate Standard Model particles as excitations around them. Early attempts to find such solutions used slightly modified models by orbifolding [93,94] or by toroidal compactification with a magnetic flux [95,96]. In Refs. [97][98][99][100][101][102], it was shown that the original model can be used to realize intersecting D-branes and Refs. [103][104][105] proposed matrix configurations that may correspond to phenomenologically viable low-energy effective theories. For related work, see also Refs. [106,107].

the Euclidean version
The Euclidean version of the IKKT matrix model, which we focus on in this paper, is obtained from the Lorentzian version by the Wick rotation The action is given by Eq. (2.1), where the contractions are now made with the metric δ µν (µ, ν = 1, . . . , 10), and the partition function is defined by  [22,23] despite the flat directions in the bosonic action S b . However, the Pfaffian Pf M = |Pf M|e iΓ is complex in general, which causes a severe sign problem in numerical simulations. The phase Γ fluctuates wildly for large matrices and plays a central role in the realization of the SSB of SO (10). In Ref. [24,25] it was shown that configurations with lower dimensions result in milder fluctuations of Γ, which points to the mechanism for favoring configurations less symmetric than SO (10).
Monte Carlo simulations of the Euclidean IKKT and related matrix models have a long history. Simplified versions can be defined by considering the reduction of the Ddimensional SYM theory to zero dimensions, which is possible for D = 3, 4, 6 and 10. The D = 3 model is ill-defined because the partition function is divergent [22,23]. The D = 4 model has a real non-negative fermion determinant in the effective action and Monte Carlo simulations showed that the SO(4) symmetry is not broken [15,17]. The D = 6 model has a complex fermion determinant and simulations are plagued by a severe sign problem as in the D = 10 case. Omitting the fermionic matrices [14] by simulating the bosonic model or omitting Γ [16,19,20] by simulating phase-quenched models, one finds no SSB.
Monte Carlo simulations including the complex phase were performed for the first time in Refs. [18-20, 33, 34]. These calculations used a reweighting-based method [18,33,34]. As an order parameter of the SSB, the eigenvalues λ µ of the "moment of inertia tensor" are defined with the ordering When the SSB of SO(D) to SO(d) with d < D is realized, the expectation values λ 1 = λ 2 = . . . = λ d are larger than the other λ µ (µ = d + 1, . . . , D). 3 The effect of the sign problem is reduced by simulating phase-quenched "microcanonical" systems with the constraints µ δ(λ µ − x µ ). Here the values of x µ are chosen appropriately assuming that the SSB of SO(D) to SO(d) with d < D is realized for each d. Varying the parameters x µ , one can sample in entropically highly suppressed regions of the configuration space in which the fluctuations of Γ are milder. The expectation values λ µ for each SO(d) symmetric vacuum are obtained by minimizing the free energy for the constrained system with respect to the parameters x µ using the saddle-point approximation, which is justified at large N.
The effect of the phase is factorized in the free energy, and it can be obtained by computing the average phase as a function of x µ . The results for λ µ obtained in this way for each SO(d) symmetric vacuum are found to be consistent with the results obtained by the GEM to be discussed below. However, comparison of the free energy for different vacua requires integration over x µ , which cannot be done accurately enough to draw a definite conclusion on the SSB pattern due to propagation of both systematic and statistical errors.
The issue of the SSB of SO(D) was also addressed analytically by using a systematic expansion called the GEM. Although the method involves only perturbative calculations, it can give nonperturbative information on the model to which it is applied [108], as we review in more detail in Section 4. In the context of matrix models, the GEM has been first applied to the BFSS matrix theory [109] and to simplified versions of the Euclidean IKKT model [110,111]. Then the SSB of rotational symmetry in the Euclidean IKKT model and related models has been investigated intensively [26][27][28][29][30][112][113][114][115][116][117]. One expands around a Gaussian action S 0 introduced by hand and containing many parameters, one for each quadratic term. In order to reduce the number of parameters, an "ansatz" that has an SO(d) symmetry is considered. The free energy and the expectation values of observables are calculated in an expansion around S 0 as functions of the parameters introduced. One can actually determine the region of the parameter space in which the free energy is independent of the parameters, and using these values of the parameters, one can obtain the free energy and the observables λ µ for each d.
The D = 6 model was studied in this way [29] and the SSB to SO(3) was found. The free energy and the spacetime extent were calculated up to the fifth order for the SO(d) ansatz with 2 ≤ d ≤ 5, and the free energy obtained for the SO(3) ansatz was found to be the minimum. The extended directions have an extent R 2 (d) = λ µ , µ = 1, . . . , d, which is equal to R 2 (3) ≈ 1.76 for d = 3. The shrunken directions have an extent r 2 = λ µ ≈ 0.223 (µ = d + 1, . . . , D), which turned out to be almost independent of d. Furthermore, the values R 2 (d) are such that they obey the constant volume property given by the relation where l is a length scale such that v ≡ l D gives the volume of spacetime. Its value was calculated and found to be l 2 ≈ 0.627. These values are consistent with the Monte Carlo simulations in Refs. [19,21]. A similar study was done also for the D = 10 model [30]. A systematic computation up to the third order was carried out for the SO(d) ansatzes with 2 ≤ d ≤ 7. The free energy for the SO(3) ansatz was found to be the minimum, suggesting also the SSB to SO (3). The values of the large and small extents of space R(d), r were calculated and found to have similar properties as in the D = 6 case. In particular, one obtains These results are consistent with the ones obtained by Monte Carlo simulations [20]. Note also that these values are obtained in the large-N limit and that they are finite. This should be contrasted with the results in the Lorentzian model, where space seems to expand indefinitely in the large-N limit.

Applying the CLM to the Euclidean IKKT model
As we reviewed in the previous section, the SSB of SO(10) symmetry in the Euclidean IKKT model is expected to occur due to the effect of the phase of the Pfaffian in (2.13), which causes the sign problem. The aim of the present work is to use the CLM to overcome this sign problem, and to understand the SSB pattern from first principles.
Here we apply the CLM to investigate the SSB of the SO(10) symmetry in the Euclidean IKKT model in a way similar to Refs. [21,53]. We discuss how to probe the SSB by using appropriate order parameters and taking appropriate limits. We also discuss important techniques used to avoid known problems in the CLM such as large excursions in the anti-Hermitian direction and the singular-drift problem [44][45][46][47][48]118]. These techniques include a deformation of the fermionic action, the adaptive stepsize and gauge cooling. By satisfying certain criteria [47,49], we can ensure that the CLM yields correct results for a large space of parameters. The discussion is brief and more details can be found in Ref. [21]. Our main results are presented in the last subsection.

the complex Langevin method
The model (2.13) we investigate can be written as where we define the effective action S eff = S b − log PfM, which is complex. In the CLM, we complexify the dynamical variables, which amounts to regarding A µ as general complex traceless matrices, and consider their fictitious time evolution governed by the complex Langevin equation, which is given as Here t is the fictitious time and η µ (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 is called the drift term, which is given as where A µ (t) is a general complex matrix solution of (3.2), t 0 is the thermalization time and T is large enough in order to obtain good statistics. Upon complexification of the matrices A µ (t), the observable O[A µ (t)] depends on general complex matrices. The analyticity of the function O[A µ ] plays a crucial role in the proof of the validity of (3.4) [44,45,47]. The numerical solution of (3.2) involves the discretization of the time t given as The square root in the last term comes from the chosen normalization of the η µ (t) so that

how to probe the SSB
In order to probe the SSB, we break the SO(10) symmetry explicitly by adding the terms to the action, where 0 < m 1 ≤ . . . ≤ m 10 , and take the ǫ → 0 limit after taking the large-N limit. As the order parameters, we consider [21,53] where no sum over µ is taken 4 . The λ µ are complex for a generic complexified configuration A µ (t), but after taking the average using Eq. (3.4), the λ µ become real. This is due to the symmetry of the drift term (3.3) under A µ −→ (A µ ) † for µ = 1, . . . , 9 and A 10 −→ −(A 10 ) † . Due to the choice of the ordering of m µ , we have λ 1 ≥ . . . ≥ λ 10 for finite N. If there is no SSB of the SO(10) rotational symmetry, all λ µ are equal in the N → ∞ and ǫ → 0 limits. If it turns out that some of them are different, we conclude that SSB occurs. µ ; t) of the general complex matrix 4 We avoid the use of the eigenvalues of T µν in Eq. (2.14) so that we do not enter into the subtleties involved in the holomorphicity of the observables. We have measured them, however, using an ordering based on their real part and the results are quantitatively identical to those obtained by using Eq. (3.7).

some techniques to make the CLM work
in the ensemble defined by P (A (R) , A (I) ; t) falls off exponentially or faster [47]. The above condition can be violated if the A µ (t) makes long excursions in the anti-Hermitian direction ("excursion problem"). In order to avoid it, we employ gauge cooling [50] in our simulations by minimizing the "Hermiticity norm" defined by at each step of (3.5). (See Ref. [21] for more details.) It was proven [47,51] that adding the gauge cooling procedure in the CLM does not affect the argument for its justification.
The stability of the evolution given by Eq. (3.5) is also controlled by the adaptive stepsize [119]. The typical magnitude of the drift term |∂S/∂(A µ ) ji |∆t ∼ u √ N ∆t in (3.5) has to be kept small compared with the typical |(A µ ) ij | (λ min /N) 1/2 , where λ min ≡ min µ |λ µ | for a given configuration. Therefore, at each step we adjust ∆t = min{∆t max , ǫ t √ λ min /(Nu) }, where ∆t max and ǫ t are tuned so that the results have negligible finite ∆t errors. Measurements are taken at fixed intervals of the Langevin time.
Another reason for the condition on the drift distribution to be violated is the singulardrift problem. This problem occurs because of the appearance of M −1 in Eq. (3.3) when the eigenvalues of M accumulate densely near zero. In order to avoid this problem, we deform the fermionic action by adding a term where m f is the deformation parameter. The addition of this term shifts the eigenvalue distribution of M away from the origin, enabling us to avoid the singular-drift problem. This technique was proposed originally in Ref. [53], where the singular-drift problem was indeed overcome in an SO(4) symmetric matrix model with a complex fermion determinant, and it was used successfully also in the D = 6 version of the IKKT matrix model [21]. Note that the above term (3.11) breaks the SO(10) symmetry down to SO(7) × SO(3) explicitly. Therefore, we investigate whether the remaining SO (7) symmetry breaks down to smaller subgroups as m f is varied and discuss what occurs at m f = 0.
Note that as m f → ∞, the fermionic degrees of freedom decouple and we obtain the dimensionally reduced Yang-Mills model without fermionic matrices ("bosonic model"), which is SO(10) symmetric [20]. Thus the deformation (3.11) can be regarded as an interpolation between the IKKT model and the bosonic model.
The presence of the singular-drift problem depends on the values of the parameters m f and ǫ. For large enough m f the problem is cured, but it reappears as m f is reduced for values of ǫ smaller than some value depending on m f . In our simulations, we choose the range of the simulated (m f , ǫ) for each N carefully so that the probability distribution of the magnitude of the drift (3.9) falls off faster than exponentially, ensuring that we do not have the singular-drift problem. In Appendix A, we show some results suggesting that the singular-drift problem vanishes for large enough N for given values of (m f , ǫ).

results of the CLM
In this subsection we present our main results obtained in the way discussed in the previous subsections. The deformation (3.11) breaks the SO(10) symmetry down to SO(7) × SO(3), and we examine if the SO(7) symmetry is broken down to a smaller group for m f = 3.0, 1.4, 1.0, 0.9, 0.7 to consider what happens for the undeformed IKKT model (3.1) corresponding to m f = 0. In order to probe the SSB, one has to take the N → ∞ limit first and then the ǫ → 0 limit.
In particular, we may confirm that the SO(3) symmetry remains unbroken by seeing that λ 1 = λ 2 and λ 3 agree in the N → ∞ and ǫ → 0 limits. On the other hand, this choice of m µ has a drawback that λ 6 and λ 7 are mixed up because of m 6 = m 7 , and hence one cannot distinguish SO(5) and SO (6) vacua. This does not cause any harm, however, as far as we find that λ 4 and λ 5 do not agree in the N → ∞ and ǫ → 0 limits, which implies that the SO(7) symmetry is broken to SO(4) or lower symmetries.
For given m µ in Eq. (3.6), the large-N limit is obtained by first computing the ratio 12) and then by making a large-N extrapolation The reason for investigating the ratio (3.12) instead of λ µ is that a large part of the ǫ dependence is canceled between the numerator and the denominator, which makes the ǫ → 0 extrapolation more reliable. The large-N extrapolation is performed by plotting ρ µ (m f , ǫ, N) against 1/N and making a quadratic fit with respect to 1/N. We find that the quadratic term is necessary in the fits especially for small (m f , ǫ). In Fig. 1 we show a typical case of such a fit. In Fig. 2 we plot the large-N extrapolated values ρ µ (m f , ǫ) as a function of ǫ for m f = 3.0, 1.4, 1.0, 0.9 and 0.7. In order to increase statistics, we average the ρ µ (m f , ǫ) for the µ that we expect to give equal values due to degeneracies in m µ . For m f = 3.0, we use m µ = (0.5, 0.5, 0.5, 1,2,4,8,8,8,8), which implies that the ρ µ (m f , ǫ) should be equal for µ = 1, 2, 3 and µ = 8, 9,10. For m f = 1.4, 1.0, 0.9, 0.7, we use m µ = (0.5, 0.5, 1,2,4,8,8,8,8,8), which implies that the ρ µ (m f , ǫ) should be equal for µ = 1, 2, µ = 6, 7 and µ = 8, 9, 10. The ǫ → 0 extrapolation is performed by fitting ρ µ (m f , ǫ) to a polynomial in ǫ.
From the extrapolated values ρ µ (m f ), we find that the SO(7) symmetry of the deformed model is not spontaneously broken at m f = 3.0, but it is actually broken to SO(4) for m f = 1.4 and to SO(3) for m f = 1.0, 0.9, 0.7. Thus as m f is decreased, the SO(7) symmetry seems to be spontaneously broken to smaller subgroups gradually in the same way as it was observed in the D = 6 case [21]. However, we consider that the symmetry is not going to be broken further down to SO (2) at smaller m f . This is based on the fact that the Pfaffian vanishes identically for strictly 2D configurations [24,25], which implies that the mechanism of SSB due to the phase of the Pfaffian no longer works there 6 . Hence our results are consistent with the results obtained by the GEM for the undeformed model, which show that the SO(3) vacuum has the smallest free energy.
When we make the fit for the ǫ → 0 extrapolation in Fig. 2, we excluded the data points at small ǫ except for m f = 3.0. These data points show a clear tendency towards the restoration of SO(10) rotational symmetry. We consider that this is due to the insufficient large-N extrapolation performed to obtain these points due to severe finite-N effects for small m f and small ǫ. As we discussed earlier, the parameter m f interpolates between the bosonic and supersymmetric models. In the bosonic model there is a strong attractive force between all the pairs of the N spacetime points defined by the eigenvalues of the A µ , which represents an O(N 2 ) effect [14], whereas in the supersymmetric model the attractive force is mostly canceled due to supersymmetry and it becomes an O(N) effect [2]. Also, a large value of ǫ reduces the spacetime volume according to Eq. (3.6). Therefore, for given N the density of spacetime points becomes small as one takes smaller values of m f and ǫ, and this is the reason why finite-N effects become severe in this region.
When finite-N effects are important, SSB is suppressed and configurations appear symmetric. The transition between the regions in which one sees true SSB effects and a (falsely) In the m f = 3.0 plot, the curves from top to bottom are (ρ 1 + ρ 2 + ρ 3 )/3, ρ 4 , ρ 5 , ρ 6 , ρ 7 and (ρ 8 + ρ 9 + ρ 10 )/3. For the other plots, the curves from top to bottom are (ρ 1 + ρ 2 )/2, ρ 3 , ρ 4 , ρ 5 , (ρ 6 + ρ 7 )/2 and (ρ 8 + ρ 9 + ρ 10 )/3. symmetric behavior can be seen as a clear smooth crossover in the three lower plots in Fig. 2 for m f = 0.7, 0.9 and 1.0. The crossover disappears for larger values of m f as expected. This is seen to happen already for m f = 3.0 in Fig. 2, whereas the crossover is very mild for m f = 1.4. The fitting region is chosen so that finite-N effects in the small-ǫ region do not affect the ǫ → 0 extrapolation.

Consistency check based on the GEM
In the previous section we discussed the pattern of the SSB of SO(10) as the deformation parameter m f in Eq. (3.11) is varied. As it has already been mentioned, m f interpolates between the IKKT matrix model (m f → 0), where SSB to SO(3) is expected to occur, and the bosonic model (m f → ∞), where there is no SSB. As the value of m f is decreased, the SSB pattern gradually changes from more symmetric vacua to less symmetric ones. However, we had to make extrapolations, first the large-N extrapolation and then the ǫ → 0 extrapolation, which introduce systematic errors. Therefore, it is important to verify the results using a completely different approach. In Section 2.2 we made a short introduction to the GEM and its application to the Euclidean IKKT matrix model. In particular, SO(10) was found to be broken down to SO(3) [30]. Here we apply the GEM to the IKKT matrix model deformed by the parameter m f . We perform a three-loop calculation using SO(d) symmetric ansatzes, where d = 6, 7, and calculate the free energy. We observe a trend that the free energy for the SO(6) vacuum becomes smaller than the SO(7) vacuum as we decrease m f . We also find that the extent of space obtained at m f = 3.0 agrees very well between the two methods.
It should be noted that the GEM does not depend on the large-N and small-ǫ extrapolations. The N → ∞ limit can be taken by simply using planar graphs, and the SSB can be readily probed by comparing the free energy obtained with various ansatzes for the Gaussian action with different symmetries. Furthermore, the systematic errors of the GEM come mainly from the truncation of the expansion and the determination of the parameters giving the physical solutions, which are completely different from the systematic errors of Monte Carlo simulations. Therefore the results of the two methods can be considered to be independent and their consistency provides more confidence in our conclusion.

applying the GEM to the deformed model
Here we review the GEM used in Ref. [30] to study the original IKKT matrix model and apply it to the deformed model with the fermionic mass-like term (3.11). The basic idea is to introduce a Gaussian action S 0 and to rewrite the original action S = S b + S f + ∆S f as where the first and second terms are regarded as the "classical action" and the "one-loop counterterm," respectively. This enables us to perform a perturbative expansion. The Gaussian action contains a large number of arbitrary parameters and the results depend on these parameters in general. On the other hand, physical quantities should not depend on these parameters. It turns out to be possible to find a range of these parameters in which observables are (almost) independent of their values. Let us consider the Gaussian action where M µ and A αβ are the parameters for the bosonic part and the fermionic part of the Gaussian action, respectively. We choose the ordering 0 < M 1 ≤ . . . ≤ M 10 so that the notation matches with (3.6). The complex 16 × 16 antisymmetric matrix A αβ in (4.2) can be expanded using the 10-dimensional gamma matrices Γ µ as where m µνρ is a totally antisymmetric 3-form. The effect of the fermion mass term (3.11) can be readily implemented by replacing the coefficient as m 8,9,10 → m 8,9,10 − m f . The partition function can be rewritten as where · 0 represents the expectation value with respect to the partition function for the Gaussian action Z 0 . One can expand the free energy F = − log Z perturbatively as where the subscript C in · C,0 means that only connected diagrams are summed over. Similarly, the expectation values of observables are given by In practice, we truncate the infinite series such as (4.4) and (4.5) at some finite order and evaluate each term using Feynman diagrams, where we restrict ourselves to planar diagrams since we are interested in the large-N limit.
As we already mentioned above, for a generic set of parameters M µ and A αβ , the results obtained in this way depend on their values. In order to find "physical" results that do not depend on these parameters, we search for stationary points of the free energy with respect to those parameters by solving the "self-consistency equations" The solutions to these equations are obtained numerically, and they are used as the probe of the plateaus for the values of the free energy in the space of the parameters M µ and A αβ . In order to consider the SO(d) symmetric vacuum, we impose the SO(d) symmetry on the Gaussian action by setting M 1 = . . . = M d and m µνρ = 0 unless µ, ν and ρ are different from each other with µ, ν, ρ ≥ d + 1. Thus we are left with (11 − d) parameters from M µ and 10−d C 3 parameters from m µνρ . However, it turns out hard to solve the self-consistency equation (4.6) with more than 5 parameters. In the previous work, we reduced the number of parameters by imposing some discrete symmetries Σ d to the shrunken directions as well, where SO(d) × Σ d are subgroups of SO (10). In the present work, we focus on two cases, in which we impose SO(d)×Z 3 with d = 6 and d = 7, where Z 3 represents a group of cyclic permutations of the 8th, 9th and 10th directions. Note that the imposed symmetries are subgroups of SO(7)×SO(3), which remains after adding the mass term (3.11).

results of the GEM
In Fig. 3 we plot the free energy calculated up to three loops for the solutions found with the SO(7) and SO(6) ansatzes against the fermion mass m f . We observe a clear tendency that the SO(6) symmetric vacuum is more favored as m f is decreased. However, the free energy for the two ansatzes tends to become degenerate as m f is increased. In this situation it is difficult to identify the critical point, given the accuracy of the GEM results.  (7) and SO (6) ansatzes are plotted against the fermion mass m f . We observe a clear tendency that the SO(6) symmetric vacuum is more favored as m f is decreased, whereas the free energy for the two ansatzes tends to be degenerate as m f is increased.
In Fig. 4, we plot the extent of space λ i (i = 1, . . . , 10) in each direction against m f for the SO(7) and SO (6) ansatzes. For the SO(7) ansatz (Left), we plot λ 1 = . . . = λ 7 , and λ 8 = λ 9 = λ 10 . We find that the two lines come close to each other, which is consistent with the fact that the deformed model becomes the bosonic model at m f = ∞, where the full SO(10) symmetry is expected to be restored. For the SO(6) ansatz (Right), we plot λ 1 = . . . = λ 6 , λ 7 and λ 8 = λ 9 = λ 10 . We find that the line in the middle corresponding to λ 7 goes up as m f increases and asymptotes to the line at the top corresponding to λ 1 = . . . = λ 6 . As a result, the plot looks almost identical to the plot for the SO (7) ansatz in the large-m f region. This clearly explains why the free energy for the two ansatz asymptotes to each other as m f increases. On the other hand, the line in the middle goes down as m f decreases and it seems to asymptote to the line at the bottom corresponding to λ 8 = λ 9 = λ 10 . This is consistent with the GEM results for the undeformed model, which suggest that the shrunken directions have a common extent. It is also natural from the viewpoint that the explicit breaking of SO(10) symmetry to SO(7) × SO(3) by the fermion mass term is removed as m f decreases.

Summary and discussions
In this paper we have applied the CLM to the Euclidean IKKT matrix model, which is conjectured to be a nonperturbative formulation of superstring theory in ten dimensions. This is the first time that a first principle study of this model produced clear results on the question of dynamical compactification of extra dimensions via SSB of the SO(10) rotational symmetry of the model. Monte Carlo simulations are plagued by a serious sign problem, which is overcome by applying the CLM. In order to avoid the singular-drift problem in the CLM, we deform the model by adding a mass-like term in the fermionic action parameterized by the parameter m f so that the conditions proposed in Ref. [47] can be met. By taking the large-N limit and then the ǫ → 0 limit, we have studied the SSB pattern as a function of the m f . At m f = 3.0 an SO(7) vacuum is found, which is the maximally symmetric vacuum of the deformed theory. At m f = 1.4 we find that our results are consistent with an SO(4) vacuum. At m f = 1.0, 0.9 and 0.7, the vacuum becomes SO(3) symmetric. Taking into account the argument that an SO(2) vacuum is not likely to be realized in this model [24,25], we conclude that our results are consistent with the ones obtained by using the GEM on the Euclidean IKKT matrix model [30], which predict SSB to an SO(3) vacuum.
As a consistency check, we performed independent calculations applying the GEM to the m f -deformed Euclidean IKKT matrix model. We did a three-loop calculation using SO (7) and SO (6) ansatzes and calculated the free energy. It turns out that the transition from an SO(7) symmetric vacuum to an SO(6) symmetric one occurs smoothly as m f decreases. Also the extent of space obtained at m f = 3.0 agrees very well between the two methods.
Our conclusion that the SO(10) rotational symmetry of the Euclidean IKKT model breaks down to SO(3) due to the phase of the Pfaffian is interesting, but it makes the model somewhat more difficult to interpret. Given the promising properties of the Lorentzian model reviewed in Section 2, we consider that the naive Wick rotation to the Euclidean model is not the right direction to pursue. On the other hand, the fact that the CLM enabled us to obtain a clear SSB pattern for the deformed model, which suffers from a severe sign problem, is encouraging. We hope that the CLM is equally useful in investigating the Lorentzian IKKT model, in particular in the presence of fermionic matrices, which are not included yet in Ref. [11].
Let us emphasize that the model has the amazing properties that spacetime, and possibly the matter content as well, are contained in the matrix degrees of freedom and that many interesting related questions can be answered dynamically. The surmounting evidence that the IKKT matrix model has nontrivial dynamics makes it a particularly promising candidate for a nonperturbative definition of superstring theory. By improving the algorithms that solve the sign problem and by acquiring more computational power that will allow us to study larger N, there is a great hope that we will find answers to many profound questions, in a similar way that it was done for QCD using the lattice gauge theory.
A The singular-drift problem vanishing at large N Rather surprisingly, we find that the singular-drift problem vanishes for large enough N for given values of (m f , ǫ). In Fig. 5 (Left), we plot the probability distribution p(u) of u defined in (3.9) for m f = 0.9 and ǫ = 0.16 with m µ = (0.5, 0.5, 1, 2, 4, 8,8,8,8,8) using various values of N within 16 ≤ N ≤ 96. We observe for N = 16 that the tail of the p(u) distribution is subexponential as a result of the singular-drift problem [47]. For N ≥ 32, the tail of the distribution is suppressed as an exponential or faster, indicating that the singular-drift problem is absent. This is confirmed by measuring the distribution p(ν) of the fourth root of the doubly degenerate non-negative eigenvalues d k of the positive definite Hermitian matrix D = MM † , where ν k = d when the ν k accumulate near the origin. The plot in Fig. 5 (Right) shows that the distribution of the ν k develops a small-ν finite cutoff ν c > 0 for N ≥ 32 proving that there is no singular-drift problem and that for N = 16, the ν k accumulate near the origin showing the appearance of the singular-drift problem.