The emergence of expanding space-time and intersecting D-branes from classical solutions in the Lorentzian type IIB matrix model

The type IIB matrix model is a promising candidate for a nonperturbative formulation of superstring theory. As such, it is expected to explain the origin of space--time and matter at the same time. This has been partially demonstrated by the previous Monte Carlo studies on the Lorentzian version of the model, which suggested the emergence of (3+1)-dimensional expanding space--time. Here we investigate the same model by solving numerically the classical equation of motion, which is expected to be valid at late times since the action becomes large due to the expansion of space. Many solutions are obtained by the gradient descent method starting from random matrix configurations, assuming a quasi-direct-product structure for the (3+1)-dimensions and the extra 6 dimensions. We find that these solutions generally admit the emergence of expanding space--time and a block-diagonal structure in the extra dimensions, the latter being important for the emergence of intersecting D-branes. For solutions corresponding to D-branes with appropriate dimensionality, the Dirac operator is shown to acquire a zero mode in the limit of infinite matrix size.


Introduction
Superstring theory has been investigated intensively as a unified theory including quantum gravity. However, in its perturbative formulation that has been established so far, there exist tremendously many stable vacua corresponding to various space-time dimensionality, gauge groups and matter contents. Also, it is known that the cosmic singularity at the beginning of the Universe cannot be resolved within perturbative superstring theory [1][2][3][4]. In order to address these issues, nonperturbative formulation is clearly needed, and the matrix models [5][6][7] have been proposed as its candidates.
In this paper, we focus on the type IIB matrix model [6], which is distinctive in that not only space but also time emerges dynamically from the matrix degrees of freedom. Indeed, it was shown by Monte Carlo simulation that (3+1)-dimensional expanding space-time appears from the Lorentzian version of the model [8]. The time scale that has been probed by the simulation is typically of the order of the Planck scale. This line of research has been continued in Refs. [9][10][11][12][13][14], whereas the Euclidean version has been investigated in Refs. [15,16] and references therein.
In order to probe the space-time at later times than the Planckian time, we need to increase the matrix size by many orders of magnitude, which is not feasible in numerical simulations.
On the other hand, it is expected that the classical approximation is valid at late times because the action becomes large due to the expansion of space. Various classical solutions that allow cosmological interpretations have been found [17][18][19][20], and a systematic method for constructing classical solutions was also developed [18]. (See also Refs. [21][22][23] for related papers in this direction.) However, the class of solutions that has been constructed so far looks too simple to describe our real world.
In fact, solving the classical equation of motion in the type IIB matrix model is quite different from ordinary classical dynamics, in which one solves a differential equation with some initial conditions. This is because the notion of time does not exist from the outset in this model, and it should emerge dynamically from the matrix degrees of freedom. Here we propose a method based on the gradient descent method, which enables us to generate infinitely many classical solutions starting from random matrix configurations. In particular, we assume a "quasi" direct-product structure for the (3+1) dimensions and the extra 6 dimensions, which is the most general structure compatible with the (3+1)-dimensional Lorentz symmetry [24]. Our ultimate goal is to determine the solution that describes our real world by identifying a specific solution that appears from the dominant configurations obtained by numerical simulations. We consider that some, if not all, of the qualitative features of such a solution should be captured by the solutions generated by our method.
In typical solutions, we find that the (3+1)-dimensional space-time exhibits an expanding behavior with a smooth structure. This is in sharp contrast to the situation with the previous Monte Carlo studies of the Lorentzian type IIB matrix model [13], where it was found that the 3-dimensional expanding space is described essentially by the Pauli matrices. This singular structure has been attributed to the approximation used to avoid the sign problem in the Monte Carlo simulation. Indeed a clear departure from the Pauli-matrix structure has been observed in Ref. [14], where the sign problem was treated correctly by the complex Langevin method for a simpler bosonic Lorentzian model. Based on this observation, it has been conjectured that a smooth (3+1)-dimensional expanding space-time should emerge dynamically from the Lorentzian type IIB matrix model in the large-N limit. Our results for the classical solutions support this conjecture.
Another issue we would like to address in this paper concerns how the Standard Model particles emerge from the type IIB matrix model at low energy. In fact, Refs. [25][26][27] discussed matrix configurations representing intersecting D-branes that can accommodate Standard Model fermions although it was not clear whether such configurations can be realized as classical solutions. Here we find that a slight extension of the simple direct-product structure, which we refer to as the "quasi" direct-product structure above, naturally realizes classical solutions with a block-diagonal structure in the extra dimensions, which corresponds to intersecting D-branes.
By choosing the dimensionality of the intersecting D-branes appropriately, we observe a trend that the Dirac operator acquires a zero mode in the large-N limit. We also confirm that the wave function of the zero mode is localized at a point, which is consistent with the picture of intersecting D-branes.
In this scenario, the Dirac zero modes appear in pairs with opposite chirality. As simple ways to obtain chiral fermions from classical solutions in the type IIB matrix model, one can either compactify the model on a noncommutative torus [28,29] or one can replace the matrices by operators [30] assuming N = ∞. However, it is still an open question whether the original type IIB matrix model can give rise to chiral fermions at low energy starting from finite N and taking the large-N limit.
The rest of this paper is organized as follows. In Sect. 2, we explain how we generate classical solutions of the type IIB matrix model assuming a quasi-direct-product structure for the (3+1) dimensions and the extra 6 dimensions. In Sect. 3, we investigate the space-time structure in the (3+1) dimensions and show that the expanding behavior with a smooth structure appears in typical solutions. In Sect. 4, we show that the block-diagonal structure in the extra dimensions that appears naturally in typical solutions can be viewed as a realization of intersecting D-branes. In particular, we show that the Dirac zero modes appear with localized wave functions when the dimensionality of the D-branes is chosen appropriately. Section 5 is devoted to a summary and discussions.

Generating classical solutions
In this section, after a brief review of the type IIB matrix model, we explain how we generate its classical solutions numerically.

The type IIB matrix model and its classical solutions
The action of the type IIB matrix model is given as [6] 2) where U ∈ SU(N ). The Euclidean version can be defined by the replacement A 0 = iA 10 , where A 10 is Hermitian, but here we stick to the Lorentzian version given above. In order to make the Lorentzian model well-defined, we impose the constraints corresponding to the IR cutoffs introduced in Ref. [8]. 1 At late times, we expect that the classical approximation is valid due to the expansion of space. We therefore solve the classical equation of motion of this model, which reads where ξ and ζ are the Lagrange multipliers corresponding to the constraints (2.5). Note here that if ξ = ζ = 0, classical solutions are always simultaneously diagonalizable, as is proved in Appendix A of Ref. [20]. Thus the existence of the constraints (2.5) is crucial in obtaining nontrivial solutions.
In this paper we search for solutions with a quasi-direct-product structure for the (3+1) dimensions and the extra 6 dimensions [24] given as where the N X × N X Hermitian matrices X µ represent the structure in the (3+1) dimensions, and the N Y × N Y Hermitian matrices Y a represent the structure in the extra 6 dimensions with This structure is actually compatible with the SO(3,1) symmetry to the direct-product structure. We will see that the matrix M = 1 plays an important role in making Y a block-diagonal, which will be interpreted later as dynamical generation of intersecting D-branes.

The algorithm
Since time does not exist a priori in the type IIB matrix model, solving the classical equation of motion (2.6) requires some method other than just solving differential equations unlike in ordinary classical dynamics. Here we define a "potential" By minimizing V using the gradient descent method starting from random configurations, we can generate classical solutions. Within our ansatz (2.7), the constraints (2.5) reduce to where Eq. (2.10) is chosen by using the redundancy under the rescaling X µ → sX µ and M → The initial configuration is prepared by generating Gaussian variables for each element of

7). Since the equation of motion (2.6) is invariant under rescaling
we make this rescaling at each step of the algorithm in such a way that the second equation of (2.5) is satisfied. The algorithm proceeds by updating X µ , Y a and M with the increment where the stepsize is taken to be small enough to make sure that V decreases at each step.
We repeat Eqs. (2.12) and (2.13) until V becomes as small as O(10 −4 ). The first equation of Eq. (2.10) is not imposed on configurations during this procedure, and we use it only to determine the value of κ after a classical solution is found. In this way, we can obtain infinitely many solutions with some κ depending on the initial random configuration and the initial values of ξ and ζ.

Typical solutions
In typical solutions generated by our algorithm, it turns out that the matrices M and Y a take the form up to the SU(N Y ) symmetry and N Y = n − + n 0 + n + . 2 This implies that Plugging Eq. (2.15) into the equation of motion (2.6) with the ansatz (2.7), we find that X µ and Y a decouple as In this work, we set the parameters ξ and ζ to be equal in order to respect the SO (3,1) symmetry and take them to be positive. Note that once we impose this condition on the initial 3 Space-time structure in the (3+1) dimensions Using the method described in the previous section, we obtain classical solutions of the type IIB matrix model within the ansatz (2.7). By generating the initial configurations randomly, we can obtain as many solutions as we wish satisfying (2.15)-(2.17). In this section, we focus on the matrices X µ in Eq. (2.7), which represent the structure of the (3+1)-dimensional space-time.

Band-diagonal structure
Using the SU(N ) symmetry (2.4) with U = g ⊗ 1, g ∈ SU(N X ), we can choose a basis in which X 0 is diagonal with its elements arranged in ascending order. In Fig. 1, we show the eigenvalues α p (p = 1, 2, . . . , N X ) of X 0 for a typical solution with N X = 64. We find that the eigenvalues are almost uniformly distributed, which suggests that we can define a continuous time in the In the basis that makes X 0 diagonal, the spatial matrices X i (i = 1, 2, 3) are not diagonal in general, but they actually turn out to be band-diagonal. In order to see it, we plot for the same solution in Fig. 2. We find that ∆ pq becomes very small for |p − q| ≥ n with some integer n, which is n ∼ 10 for the configuration in Fig. 2. The band-diagonal structure guarantees the locality of time, which enables us to extract the time evolution as we discuss in the next section.
The uniform distribution of α p and the band-diagonal structure in X i are common features of the solutions we obtain. They are actually shared also by the dominant configurations generated by previous Monte Carlo studies [9][10][11][12]. As we mentioned below (2.6), the classical solutions are strictly diagonal for ξ = ζ = 0. Therefore, the band-diagonal structure may be viewed as a deviation from the diagonal matrix due to finite ξ and ζ.

Extracting the time evolution
Next we discuss how we extract the time evolution from a given classical solution. The banddiagonal structure discussed in the previous section plays a crucial role here.
Let us define an n × n blockX i (t) by (X i (t)) rs = (X i ) k+r,k+s (r, s = 1, 2, . . . , n) , where k = 0, 1, . . . , N X − n and the argument t This block matrixX i (t) represents the space structure at time t, from which one can obtain the time evolution of the 3-dimensional space.
For instance, we define the extent of space at t by where the symbol "tr" represents the trace over an n×n matrix. In Fig. 3, we plot R 2 (t) for the typical solution discussed in Sect. 3.1. We find that R 2 (t) is roughly symmetric under t → −t reflecting the symmetry of the model under A 0 → −A 0 . In particular, an expanding behavior is seen at t > 0, which is consistent with the results of simulations [14] at late time.
Next we define a 3 × 3 real symmetric tensor whose eigenvalues λ i (t), which we order as λ 1 (t) ≤ λ 2 (t) ≤ λ 3 (t), represent how the 3dimensional space extends in each direction. Note that they are related to the extent of space In Fig. 4, we plot the eigenvalues λ i (t) for the typical solution discussed above. The three eigenvalues λ i (t) turn out to be roughly equal at each time. We have also confirmed that the three eigenvalues tend to come closer to each other as one increases the matrix size N X , suggesting that the SO(3) symmetry is recovered in the limit of infinite matrix size.  There are ten of them reflecting the chosen block size n = 10.
In order to probe the structure of the 3-dimensional space emerging from the classical solutions in more detail, we define an n × n matrix whose eigenvalues q r (t) represent the radial distribution of the points which describe the 3dimensional space. Note that the eigenvalues q r (t) are related to the extent of space R 2 (t) as R 2 (t) = 1 n n r=1 q r (t) . (3.8) In Fig. 5, we plot the eigenvalues of Q(t) for the same typical solution discussed above. We find that they are densely distributed at each time and that all the eigenvalues are growing with time in the t > 0 region. This is in sharp contrast to the results for the configurations obtained by the previous Monte Carlo studies of the Lorentzian type IIB matrix model [13], where only two of the eigenvalues grow with time, while the others remain small and constant.
As we mentioned in the Introduction, the latter behavior is interpreted as a consequence of the approximation used in the Monte Carlo studies to avoid the sign problem. Our results for the classical solutions suggest the emergence of smooth (3+1)-dimensional expanding space-time in the large-N limit.

The emergence of intersecting D-branes
In this section we discuss how chiral zero modes appear from the type IIB matrix model. A general discussion on the Dirac equation implies that it is necessary to have zero modes in the 6-dimensional Dirac operator. In fact, we show that the block-diagonal structure (2.14) of the typical solutions can be regarded as the emergence of intersecting D-branes. Indeed we find that zero modes appear in the limit of infinite matrix size if the dimensionality of the branes is chosen appropriately. where Ψ satisfies the Weyl condition

The Dirac equation
with the chirality operator Γ χ in 10 dimensions. Here we decompose the 10-dimensional gamma matrices Γ M as where ρ µ and γ a are the 4-dimensional and 6-dimensional gamma matrices, respectively. Note that the chirality operator Γ χ is also decomposed as where ρ χ and γ χ are the chirality operators in 4 and 6 dimensions, respectively.
Let us next decompose Eq. (4.1) into two terms as Note that the first term and the second term have opposite chirality in 4 dimensions as well as in 6 dimensions, which implies that each term has to vanish separately. Thus, in order to obtain chiral fermions in four dimensions, we need to have Dirac zero modes in 6 dimensions; i.e., Γ a [A a , Ψ] = 0 . (4.8) Let us now consider that A M is a classical solution with the quasi-direct-product structure (2.7). Since A a = 1 ⊗ Y a , the general solution to Eq. (4.8) can be obtained by decomposing Ψ as where the 4-dimensional and 6-dimensional gamma matrices act only on ψ (4d) and ψ (6d) , respectively. Thus, in order to satisfy Eq. For that purpose, it suffices to consider the solutions of (2.17) with a block-diagonal struc- where the block matrices Y Y , respectively. Corresponding to the block-diagonal structure (4.11) of Y a , we focus on the off-diagonal block ϕ in for which the eigenvalue problem reduces to to see in which case we obtain zero modes in the limit of infinite matrix size. Also, we look into the structure of the eigenvectors ϕ α , which we call the "wave functions" in what follows, and see whether it is localized as is expected from the picture of intersecting D-branes.
Since the Dirac operator under consideration anti-commutes with γ χ , the eigenvalues λ and −λ appear in pairs. If ϕ is an eigenvector corresponding to the eigenvalue λ, then γ χ ϕ is an eigenvector corresponding to the eigenvalue −λ. For the zero modes corresponding to λ = 0, one can choose them to have definite chirality γ χ ϕ = ±ϕ. In order to see the localization of the zero modes, we need to consider the right-handed and left-handed components separately, which can be extracted as (4.14) Furthermore, we should take into account that the classical solution of the form Eq. (4.11) has symmetry under Y ), which preserves the structure (4.11). Under this transformation, the wave functions of the corresponding Dirac operator transform as (4.17) We choose g (1) and g (2) in (4.17) in such a way that ϕ R1 is diagonal with positive elements in the descending order. This is nothing but the singular value decomposition of ϕ R1 , where the diagonal elements correspond to the singular values.
In what follows, we change the index a of Y a from a = 4, . . . , 9 to a = 1, . . . , 6. From the picture of the intersecting D-branes, we expect the emergence of Dirac zero modes when the two branes intersect at a point. This requires generically that the dimensionality of the branes adds up to six. If the sum of the dimensionality is less than six, the two branes do not intersect, and if the sum is more than six, the two branes intersect but not at a point. In order to specify the dimensionality of the brane, say to d, we set 6 − d components of Y a to zero. Under these conditions, we generate 3d and 4d solutions of (2.17) with ζ = 1 numerically, and use them as the block matrices in (4.11) in the following analysis. The adopted algorithm is the same as the one described in Sect. 2 except that only the Y a matrices are involved in our calculation here in order to reduce the computational cost.
The 2d solutions, on the other hand, can be constructed analytically by defining Z as   Eqs. (4.18) and (4.19) imply where L i (i = 1, 2, 3) are some representation of the SU(2) algebra. Here we restrict ourselves to the irreducible representation without loss of generality, since the Dirac spectrum for the reducible representation simply becomes the sum of those for the irreducible representations which the reducible representation decomposes to. The above construction suggests that the 2d brane is something like a "fuzzy disk", which can be obtained by projecting a fuzzy sphere on to a plane. More precisely, it should be regarded as two coinciding fuzzy disks corresponding to the two hemispheres of the fuzzy sphere.

The cases with intersection at a point
a , we set Y  In Fig. 7, we plot |(ϕ R1 ) pq | 2 and [(ϕ L5 ) pq | 2 , the component of the wave functions that has nonzero elements. We find that the wave functions are localized at the (1, 1) element. These results are consistent with the picture that the left-handed and right-handed zero modes appear from the intersecting D-branes.   wave function almost vanishes except for the spinor component α = 1. Similarly, we calculate |(ϕ Lα ) pq | 2 for each α, and find that the wave function almost vanishes except for the spinor component α = 5. In Fig. 9, we plot |(ϕ R1 ) pq | 2 and |(ϕ L5 ) pq | 2 , the component of the wave functions that has nonzero elements. We find that the wave functions are localized at the (1, 1) element although the localization is not as sharp as in the "3d-3d" case. These results are consistent with the picture that the left-handed and right-handed zero modes appear from the intersecting D-branes.

The cases without intersection
to the case with the "2d-3d" ansatz below. As in the "2d-4d" case discussed above, we have two-fold degeneracy in the eigenvalue spectrum. In Fig. 10, we plot the ratio µ 0 /µ 1 against 1/N (1) Y = 32, 48 and 64. We find that it does not converge to 0 in the N (1) Y → ∞ limit, which implies that we do not obtain zero modes in this case.

The cases with intersection but not at a point
Y = 32, 48 and 64. They do not converge to 0 in the N (1) Y → ∞ limit, which implies that we do not obtain zero modes in this case.
Next we consider the wave function corresponding to the lowest eigenvalue for one of the 16 cases with N (1) Y = 64. We calculate |(ϕ Rα ) pq | 2 for each α and find that the wave function almost vanishes except for the spinor component α = 1. Similarly, we calculate |(ϕ Lα ) pq | 2 for each α, and find that the wave function almost vanishes except for the spinor component α = 5.
In Fig. 13, we plot |(ϕ R1 ) pq | 2 and |(ϕ L5 ) pq | 2 , the component of the wave functions that have nonzero elements. We find that the wave functions are not localized but have a long tail along the diagonal line. These results are consistent with the picture that the intersection does not occur at a point. The zero modes do not appear, and the wave functions corresponding to the lowest eigenvalue do not localize. Similar behaviors are observed for the "4d-4d" ansatz.

Summary and discussions
In this paper, we have proposed a numerical method which enables us to solve the classical equation of motion of the type IIB matrix model. In particular, we impose a quasi-directproduct structure (2.7), which is the most general ansatz compatible with the SO(3,1) Lorentz invariance. Based on the solutions obtained in this way, we have investigated the space-time structure in the (3+1) dimensions and the structure in the extra 6 dimensions.
First, we have focused on the (3+1)-dimensional space-time structure, which is represented by X µ in Eq. (2.7). When X 0 is diagonalized, X i become band-diagonal, which enables us to extract the time evolution. We find that the 3-dimensional space expands with time, and moreover it is smooth unlike the singular structure observed by simulations with an approximation to avoid the sign problem [13]. Our results support the conjecture in Ref. [14] that the singular structure is caused by the approximation and that a smooth space-time should emerge in the N → ∞ limit if the sign problem is treated properly.
Next, we have focused on the structure in the extra dimensions. In fact, the existence of the matrix M = 1 in the quasi-direct-product structure (2.7) realizes naturally the appearance of intersecting D-branes represented by a block-diagonal structure in the extra dimensions.
In order to confirm this picture, we have investigated the eigenvalue spectrum of the Dirac operator for two intersecting D-branes with various dimensionality. In the "3d-3d" and "2d-4d" cases, our solutions give rise to Dirac zero modes in the limit of infinite matrix size. We have also found that the wave functions corresponding to the lowest eigenvalue are localized at a point, which is consistent with the picture that the zero modes appear from the intersection point. It should be emphasized that the zero modes were obtained for classical solutions of the type IIB matrix model unlike the previous studies [24,25,27,29], where the configurations that give rise to zero modes are constructed by hand. In the other cases in which the intersection either does not occur or occurs but not at a point, we do not obtain zero modes.
We consider that it is worth while to extend the present studies to larger matrix size. By obtaining X µ with larger N X , we can determine how the space expands with time for a longer time period. In particular, it would be interesting to see whether the expansion is exponential or obeys a power law possibly depending on time. It is also important to investigate whether the 3dimensional space becomes uniform and SO(3) symmetric. By obtaining Y a with larger N Y , we should be able to see the emergence of zero modes more explicitly. It would be also interesting to identify the zero modes in the fluctuation of the matrices Y a around the classical solution, which would correspond to the Higgs particles. Then we can calculate the Yukawa couplings from the overlap of the wave functions of the Dirac zero modes and the Higgs particles. Based on the information obtained in this way at the Planck scale, we may use the renormalization group to see whether the Standard Model appears at low energy. One of the most important issues to be addressed is whether one can obtain chiral fermions. We hope to report on the progress in these directions in the near future.