Complex Langevin analysis of the space-time structure in the Lorentzian type IIB matrix model

The Lorentzian type IIB matrix model has been studied as a promising candidate for a nonperturbative formulation of superstring theory. In particular, the emergence of (3+1)D expanding space-time was observed by Monte Carlo studies of this model. It has been found recently, however, that the matrix configurations generated by the simulation is singular in that the submatrices representing the expanding 3D space have only two large eigenvalues associated with the Pauli matrices. This problem has been attributed to the approximation used to avoid the sign problem in simulating the model. Here we investigate the model using the complex Langevin method to overcome the sign problem instead of using the approximation. Our results indicate a clear departure from the Pauli-matrix structure, while the (3+1)D expanding behavior is kept intact.


Introduction
Nonperturbative studies often provide totally new perspectives in quantum field theories. Confinement of quarks, for instance, has been vividly demonstrated by the strong coupling expansion in the lattice gauge theory [1]. We consider the emergence of (3+1)D expanding space-time in the Lorentzian version of the type IIB matrix model is another case of this sort [2]. This model was conjectured to be a nonperturbative formulation of superstring theory [3] analogous to the lattice gauge theory in QCD. The model has ten bosonic N × N Hermitian matrices, which represent ten-dimensional space-time in the large-N limit. One of the most interesting features is that the eigenvalue distribution of the ten bosonic matrices can collapse to a lower-dimensional manifold, which may be interpreted as the actual spacetime dynamically generated in this model. If this really happens, it implies that the (9+1)D Lorentz symmetry of the model is spontaneously broken.
Monte Carlo studies of the type IIB matrix model are extremely hard, however, due to the so-called sign problem caused by the complex weight in the partition function. In the Eulicidean version, it comes from the Pfaffian that is obtained by integrating out fermionic matrices, while in the Lorentzian version, it comes from the phase factor e iS b with the bosonic action S b . If we treat the phase of the complex weight by reweighting, huge cancellation among configurations with different phases occurs, which makes the calculation impractical. Recently the complex Langevin method (CLM) [4,5] has been attracting much attention as a promising approach to this problem [6][7][8][9][10][11][12]. In particular, it has been applied successfully to the Euclidean version of the 6D type IIB matrix model [13], and the spontaneous breaking of the rotational SO(6) symmetry to SO(3) suggested by the Gaussian expansion method [14] has been confirmed.
In ref. [2] and our subsequent work [15][16][17][18] on the Lorentzian type IIB matrix model and its simplified models, the sign problem was avoided by integrating out the scale factor of the bosonic matrices by hand, which yields a function of the bosonic action S b sharply peaked at the origin. Approximating this function by a sharply peaked Gaussian function, we can perform Monte Carlo simulations without the sign problem. The emergence of (3+1)D expanding space-time was obtained in this way [2]. The expanding behavior for a longer time was investigated by simulating the simplified models. The obtained results suggested a scenario for the full model that the expansion is exponential at early times [15], which is reminiscent of the inflation, and that it turns into a power law [16] at later times, which is reminiscent of the Friedmann-Robertson-Walker universe in the radiation dominated era. See also refs. [19][20][21] for closely related work.
It has been found recently, however, that the matrix configurations generated by the simulation is singular in that the submatrices representing the expanding 3D space have only two large eigenvalues associated with the Pauli matrices [22]. This problem has been attributed to the aforementioned approximation used to avoid the sign problem since the function obtained after integrating out the scale factor is actually complex-valued, and the effect of the phase is not taken into account. It was realized that the approximation actually amounts to replacing the phase factor e iS b by a positive definite weight e cS b with some constant c > 0. This new interpretation of the simulation provides clear understanding of the observed Pauli-matrix structure and the (3+1)D expanding behavior. It has also been argued that a regular space-time may be obtained if the phase factor e iS b is used correctly. This is a very nontrivial issue, however, since losing the Pauli-matrix structure may also imply losing the (3+1)D expanding behavior at the same time.
In this paper we address this issue by using the CLM to solve the sign problem instead of using the aforementioned approximation. Note that the Lorentzian type IIB matrix model needs to be regularized in some way or another because the phase factor e iS b in the partition function cannot suppress the contribution from the bosonic matrices with arbitrary large elements. 3 Here we use the infrared cutoffs on both the spatial and temporal matrices analogous to the ones used in the previous work [2]. We also find it useful to introduce two deformation parameters (s, k), which correspond to the Wick rotations on the worldsheet and in the target space, respectively. These parameters enable us to interpolate between the Lorentzian version (s, k) = (0, 0) and the Euclidean version (s, k) = (1, 1) .
First we focus on (s, k) = (−1, 0) in the deformation parameter space, where we do not have the sign problem. In fact, this case corresponds to the approximate model investigated in our previous work. We observe the emergence of (3+1)D expanding space-time with the Pauli-matrix structure. Then we tune the worldsheet deformation parameter s close to that for the Lorentzian model (s = 0) keeping the target space deformation parameter k in such a way that the space-time noncommutativity is minimized. There, we find it possible to obtain a smoother space-time structure without losing the (3+1)D expanding behavior. The deviation from the Pauli-matrix structure was not seen for the matrix size N ≤ 64 within the parameter region that can be explored by the CLM, and it becomes more prominent as we increase N from 128 to 192. We consider that the two deformation parameters s and k should be tuned eventually to (s, k) = (0, 0) in the large-N limit. Whether a smooth classical space-time picture appears in that limit at sufficiently late time is an important open question, which can be answered along the line of this research.
The rest of this paper is organized as follows. In section 2 we define the Lorentzian 3 This situation is in sharp contrast to the Euclidean version [23,24], in which the phase factor e iS b is is a real non-negative quantity. type IIB matrix model and introduce the infrared cutoffs as well as the two deformation parameters s and k. In section 3 we discuss how we apply the CLM to the Lorentzian type IIB matrix model. In section 4 we focus on (s, k) = (−1, 0) in the deformation parameter space, which corresponds to the approximate model investigated in the previous work. Indeed we observe the emergence of (3+1)D expanding space-time with the Pauli-matrix structure. In section 5 we show our results for the worldsheet deformation parameter s close to that for the Lorentzian model (s = 0) with the target space deformation parameter k chosen in such a way that the space-time noncommutativity is minimized. We observe a clear departure from the Pauli-matrix structure, while the (3+1)D expanding behavior is still being observed. Section 6 is devoted to a summary and discussions.

Definition of the Lorentzian type IIB matrix model
The action of the Lorentzian type IIB matrix model is given by [3] 2) where the bosonic variables A µ (µ = 0, . . . , 9) and the fermionic variables Ψ α (α = 1, . . . , 16) are N × N Hermitian matrices. Γ µ are 10D gamma-matrices after the Weyl projection and C is the charge conjugation matrix. The "coupling constant" g is merely a scale parameter in this model since it can be absorbed by rescaling A µ and Ψ appropriately. The indices µ and ν are contracted using the Lorentzian metric η µν = diag (−1, 1, . . . , 1). The Euclidean version can be obtained by making a "Wick rotation" A 0 = iA 10 , where A 10 is Hermitian. The partition function for the Lorentzian version is proposed in ref. [2] as with the action (2.1). The "i" in front of the action is motivated from the fact that the string worldsheet metric should also have a Lorentzian signature. By integrating out the fermionic matrices, we obtain the Pfaffian which is real unlike in the Euclidean case [25]. Note also that the bosonic action (2.2) can be written as where we have introduced the Hermitian matrices Since the two terms in the last expression have opposite signs, S b is not positive semi-definite, and it is not bounded from below. In order to make the partition function (2.4) finite, we need to introduce infrared cutoffs in both the temporal and spatial directions, for instance, as We can use the SU (N) symmetry of the model to bring the temporal matrix A 0 into the diagonal form By "fixing the gauge" in this way, we can rewrite the partition function (2.4) as where ∆(α) is the van der Monde determinant. The factor ∆(α) 2 in (2.10) appears from the Fadeev-Popov procedure for the gauge fixing, and it acts as a repulsive potential between the eigenvalues α i of A 0 . We can extract a time-evolution from configurations generated by simulating (2.10). A crucial observation is that the spatial matrices A i have a band-diagonal structure in the SU(N) basis in which A 0 has the diagonal form (2.9). More precisely, there exists some integer n such that the elements of spatial matrices (A i ) ab for |a − b| > n are much smaller than those for |a − b| ≤ n. Based on this observation, we may naturally consider n × n submatrices of A i defined as where I, J = 1, . . . , n, ν = 0, 1, . . . , N − n, and t is defined by We interpret theĀ i (t) as representing the state of the universe at time t. UsingĀ i (t), we can define, for example, the extent of space at time t as where the symbol tr represents a trace over the n × n submatrix. We also define the "moment of inertia tensor" which is a 9 × 9 real symmetric matrix. The eigenvalues of T ij (t), which we denote by λ i (t) with the order represent the spatial extent in each of the nine directions at time t. Note that the expectation values λ i (t) tend to be equal in the large-N limit if the SO(9) symmetry is not spontaneously broken. This is the case at early times of the time-evolution. After a critical time t c , on the other hand, it was found [2] that the three largest eigenvalues λ i (t) (i = 1, 2, 3) become significantly larger than the rest, which implies that the SO(9) symmetry is spontaneously broken down to SO(3).
Here we introduce two deformation parameters s and k, which correspond to Wick rotations on the worldsheet and in the target space, respectively. Let us introduceS = −iS b so that the factor e iS b in the partition function (2.10) is rewritten as e −S . We introduce the first parameter s (−1 ≤ s ≤ 1) corresponding to the Wick rotation on the worldsheet as The second parameter k (0 ≤ k ≤ 1) corresponding to the Wick rotation in the target space can be introduced by the replacement A 0 → e −ikπ/2 A 0 . The action (2.17) becomesS Note that the coefficient of the first term in (2.18) can be made real non-negative by choosing the parameters so that ie isπ/2 e −ikπ = 1, which implies k = (1+s)/2. For this choice, the bosonic action is most effective in minimizing the noncommutativity between the spatial matrices A i and the temporal matrix A 0 . For 0 ≤ k < s/2, on the other hand, the real part of the coefficient becomes negative, which favors maximum noncommutativity between A i and A 0 . As a result, the eigenvalues of A 0 lump up into two clusters separated from each other, and we cannot obtain a continuous time. The Lorentzian model (s, k) = (0, 0) lies on the boundary of this unphysical region. In this work, we keep away from this region by restricting ourselves to the cases satisfying k = (1 + s)/2.
Taking into account the infrared cutoffs (2.7) and (2.8), we arrive at the partition function where θ(x) is the Heaviside step function andS is given by (2.18). By rescaling A µ → LA µ and β → L −4 A µ , we can set L = 1 without loss of generality.

Applying the CLM to the Lorentzian model
We apply the CLM to the model (2.19). From now on, we omit the Pfaffian and consider the 6D version, which consists of A 0 and A i (i = 1, · · · , 5), for simplicity. The first step of the CLM is to complexify the real variables. As for the spatial matrices A i , we simply treat them as general complex matrices instead of Hermitian matrices. As for the temporal matrix A 0 , which is diagonalized as (2.9), we have to take into account the ordering of the eigenvalues. For that purpose, we make the change of variables as so that the ordering is implemented automatically, and then complexify τ a (a = 1, · · · , N − 1). We have chosen to set α 1 = 0 using the shift symmetry A 0 → A 0 + const.1 of the action. In order to respect this symmetry, we decide to impose the cutoff like (2.7) only on the traceless partÃ 0 = A 0 − 1 N Tr A 0 in this work. The Heaviside function in (2.19) is difficult to treat in the CLM as it is. Here we mimic its effect by introducing the potential where the power p is set to p = 4 in this work, and the coefficients γ s and γ t are chosen to be large enough to make 1 N Tr (A i ) 2 and 1 N Tr (Ã 0 ) 2 fluctuate around some constants. 4 The effective action then reads where the last term comes from the Jacobian associated with the change of variables (3.1). The complex Langevin equation is given by The expectation values of observables can be calculated by defining them holomorphically for complexified τ a and A i and taking an average using the configurations generated by solving the discretized version of (3.4) for sufficiently long time. In order for this method to work, the probability distribution of the drift terms, namely the first terms on the righthand side of (3.4), has to fall off exponentially [11]. We have checked that this criterion is indeed satisfied for all the values of parameters used in this paper.

Emergence of (3+1)D expanding behavior
In this section we consider (s, k) = (−1, 0) in the parameter space. The action is given by which is real, and the CLM reduces to the ordinary Langevin method. The first term in (4.1) tries to minimize the space-time noncommutativity, which has the effects of making the spatial matrices close to diagonal in the basis (2.9). On the other hand, the second term favors maximal noncommutativity among spatial matrices. Figure 1 shows our results 5 for N = 128, κ = 0.13, β = 2. The block size for (2.12) is chosen to be n = 16. In the Top panel, we plot the extent of space R 2 (t) defined by (2.14) against t. The result is symmetric under the reflection t − t p → −(t − t p ), where t p represents the time at which R 2 (t) is peaked, due to the symmetry of the model under A 0 → −Ã 0 .
Next we discuss the SSB of SO(5) symmetry by considering the moment of inertia tensor (2.15). In the Bottom-Left panel, we plot the eigenvalues λ i (t) of T ij (t), which shows that only three out of five eigenvalues become large in the time region around t = t p . This suggests that the rotational SO(5) symmetry of the 6D bosonic model is broken down to SO (3) in that time region. These results are qualitatively the same as what has been obtained in ref. [16], which is consistent with the speculation [22] that the previous simulations correspond to the parameter choice (s, k) = (−1, 0).
As is known from the previous work [2], the time difference between the peak (t = t p ) and the critical time at which the SSB occurs increases in physical units as we take the large-N limit. Therefore, the reflection symmetry with respect to t does not necessarily imply that the Big Crunch occurs in the finite future.
The mechanism of this SSB can be understood as follows [22]. Since the first term in (4.1) favors A i close to diagonal, we may consider the submatricesĀ i (t) as the effective degrees of freedom. The infrared cutoff (2.8) fixes Tr {Ā i (t)} 2 to some constant, and the second term in (4.1) favors maximal noncommutativity betweenĀ i (t). According to the argument in ref. [2], this leads toĀ i (t) ∝ σ i ⊕ 0 n−2 for i = 1, 2, 3 andĀ i (t) = 0 n for i ≥ 4 up to SO(5) rotations, where σ i are the Pauli matrices. In order to confirm this mechanism, we calculate the matrix and plot the four largest eigenvalues of Q(t) in Fig. 1 (Bottom-Right). Indeed we find that only two of them are large, while the rest are very small in the time region in which the SSB occurs.

Departure from the Pauli-matrix structure
In this section we tune the worldsheet deformation parameter s to some values near s = 0, which is the target value for the Lorentzian model, keeping the target-space deformation parameter k to be k = (1 + s)/2, which minimizes the space-time noncommutativity. The action reads The only difference from (4.1) is the second term with the coefficient e −i π 2 (1−s) whose real part changes its sign at s = 0. This implies, in particular, that for s > 0 the second term starts to minimize the noncommutativity among the spatial matrices. Therefore, we may anticipate a drastic change of the behavior around s = 0. In fact, for the values of s below what is reported below, we do not see any qualitative difference from the results obtained at s = −1.  Unlike the (s, k) = (−1, 0) case, the action becomes complex for s > −1 in general. Therefore, the quantity such as R 2 (t) defined in (2.14) is not guaranteed to be real positive. In the Top panel, we plot the real and imaginary parts of R 2 (t). We find that R 2 (t) is dominated by the real part near the peak.
Let us also take a look at the "Hermiticity norm" forĀ i (t) defined by using a configuration generated by the simulation. The result is plotted in the Top panel as well. Note that h(t) = 0 implies that the matricesĀ i (t) are all Hermitian, while h(t) = 1 implies that they are all anti-Hermitian. We find that h(t) is small and hence theĀ i (t) are close to Hermitian near the peak, which is consistent with our observation that R 2 (t) is dominated by the real part in this region. This property supports our previous speculation [ 26,27] that some classical solution, which is typically represented by a real configuration, dominates the path integral in the time region near the peak due to the expansion of space.
In the Bottom panels, we plot the same quantities 6 as in Fig. 1 (Bottom). From the left panel, we observe that (3+1)D expanding behavior persists even at slightly positive s, while the right panel reveals a clear departure from the Pauli-matrix structure.
We perform a similar analysis with the matrix size N increased from N = 128 to N = 192. Figure 3 shows our results for (s, k) = (0.0118, 0.5059), N = 192, κ = 0.0044, β = 64, n = 24. In the Top panel, we plot the real and imaginary parts of R 2 (t) and the Hermiticity norm h(t) defined by (5.2). As we have seen in Fig. 2 (Top), the spatial matricesĀ i (t) are close to Hermitian near the peak of R 2 (t), which suggests that the behavior in this region is semi-classical. In the Bottom panels, we plot the same quantities as in Fig. 2 (Bottom). We find that the departure from the Pauli-matrix structure is even more pronounced 7 , while the (3+1)D expanding behavior is kept intact.
We have also investigated the model with N < 128. For N = 32, 64, while the results at (s, k) = (−1, 0) are similar to those for N = 128, the departure from the Pauli-matrix structure does not show up at all even at s ∼ 0. As we increase s further with k = (1+s)/2, the Hermiticity of the configurations is completely lost, and the criterion for justifying the CLM is found to be violated. Thus the use of large values of N seems to be crucial in investigating the model near the target values (s, k) = (0, 0).

Summary and discussions
The Lorentzian type IIB matrix model is a promising candidate for a nonperturbative formulation of superstring theory. Monte Carlo studies of this model is extremely hard due to the sign problem caused by the phase factor e iS b in the partition function. Previous work avoided this problem by integrating out the scale factor of the bosonic matrices and using an approximation. However, it was noticed recently that this approximation actually amounts to replacing e iS b by e cS b for some c > 0. This suggests the importance of studying the model without such an approximation.
In this paper we have investigated the space-time structure based on the complex Langevin simulation of the (5+1)D bosonic version of the model with the deformation parameters s and k corresponding to the Wick rotations on the worldsheet and in the target space, respectively, The original model corresponds to (s, k) = (0, 0), whereas our previous simulations were speculated to correspond to the (s, k) = (−1, 0) case [22]. Our results for (s, k) = (−1, 0) indeed reproduced the (3+1)D expanding behavior with the Pauli-matrix structure as expected. Then we tuned the parameter s towards the region s ∼ 0 restricting ourselves to k = (1 + s)/2 in order to stabilize our simulation. The results indeed showed a clear departure from the Pauli-matrix structure, while the (3+1)D expanding behavior is kept intact. The spatial matrices turn out to be close to Hermitian near the peak of the spatial extent R 2 (t) even for s ∼ 0, which confirms our expectation [26,27] that some classical solution dominates at late times.
The appearance of the Pauli-matrix structure for (s, k) = (−1, 0) is due to the Tr (F ij ) 2 term in the action, which tries to make the spatial matrices A i maximally noncommutative. The situation changes drastically around s = 0, where the coefficient of the Tr (F ij ) 2 term becomes pure imaginary. On the other hand, there are infinitely many classical solutions [26,27], which have (3+1)D expanding behavior without the Pauli-matrix structure. (See also refs. [28][29][30][31][32][33] for related work.) We therefore consider it possible that the space-time structure becomes smooth without losing the (3+1)D expanding behavior in the large-N limit. The departure from the Pauli-matrix structure observed at s ∼ 0 supports this possibility. Some future directions are in order. The most important thing to do is to repeat the same analysis with increased matrix size N. In particular, we need to confirm the appearance of a smooth space-time at (s, k) ∼ (0, 0). While this issue may not depend much on the effects of the fermionic matrices, it would be certainly desirable to include them eventually. Unfortunately, this is not straightforward since the complex Langevin method may suffer from the singular-drift problem due to the near-zero eigenvalues of the Dirac operator. The deformation technique [12] used successfully in studying the Euclidean version [13] is worth trying, though. We consider that the dominance of classical solutions at late times [26,27] supported by our results is important because it enables us to understand possible late-time behaviors of this model by solving classical equations of motion. For instance, we may try to find classical solutions [34][35][36][37], which can accommodate Standard Model particles as excitations around them. Work in this direction is ongoing [38].
We hope that the simulation method as well as the obtained results discussed in this paper is useful in understanding the dynamics of the Lorentzian matrix model further.