---
abstract: '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 decent method starting from random matrix configurations, assuming a quasi direct-product structure for the (3+1)-dimensions and the extra six dimensions. We find that these solutions generally admit the emergence of expanding space-time and intersecting D-branes, the latter being represented by a block-diagonal structure in the extra dimensions. 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.'
bibliography:
- 'ref\_spacetime\_chiralzeromode\_classicalsolution\_IIBMM.bib'
---

  -------------
    KEK-TH-2169
  -------------

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

Kohta H[atakeyama]{}$^{1,2)}$, Akira M[atsumoto]{}$^{3)}$, Jun N[ishimura]{}$^{3,4)}$,\
Asato T[suchiya]{}$^{1,2)}$ Atis Y[osprakob]{}$^{3)}$\

$^{1)}$ [*Department of Physics, Shizuoka University*]{}\
[*836 Ohya, Suruga-ku, Shizuoka 422-8529, Japan*]{}\
$^{2)}$ [*Graduate School of Science and Technology, Shizuoka University*]{}\
[*3-5-1 Johoku, Naka-ku, Hamamatsu 432-8011, Japan*]{}\
$^{3)}$ [*Department of Particle and Nuclear Physics, School of High Energy Accelerator Science,*]{}\
[*Graduate University for Advanced Studies (SOKENDAI)*]{}\
[*1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan*]{}\
$^{4)}$ [*KEK Theory center, High Energy Accelerator Organization*]{}\
[*1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan*]{}

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 group and matter contents. Also it is known that the cosmic singularity at the beginning of the Universe cannot be resolved within perturbative superstring theory [@Lawrence:2002aj; @Liu:2002kb; @Horowitz:2002mw; @Berkooz:2002je]. In order to address these issues, nonperturbative formulation is clearly needed, and the matrix models [@Banks:1996vh; @Ishibashi:1996xs; @Dijkgraaf:1997vv] have been proposed as its candidates.

In this paper, we focus on the type IIB matrix model [@Ishibashi:1996xs], 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 [@Kim:2011cr]. 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. [@Ito:2013ywa; @Ito:2015mxa; @Ito:2017rcr; @Azuma:2017dcb; @Aoki:2019tby; @Nishimura:2019qal], whereas the Euclidean version has been investigated in Refs. [@Anagnostopoulos:2013xga; @Anagnostopoulos:2017gos] 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 [@Kim:2011ts; @Kim:2012mw; @Chaney:2015ktw; @Steinacker:2017vqw], and a systematic method for constructing classical solutions was also developed [@Kim:2012mw]. (See also Refs. [@Steinacker:2019awe; @Steinacker:2019fcb; @Manolakos:2019tdf] 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 decent 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 six dimensions, which is the most general structure compatible with the (3+1)-dimensional Lorentz symmetry [@Nishimura:2013moa]. 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 [@Aoki:2019tby], 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. [@Nishimura:2019qal], 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. [@Chatzistavrakidis:2011gs; @Steinacker:2014fja; @Aoki:2014cya] 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 non-commutative torus [@Aoki:2008ik; @Aoki:2010gv] or one can replace the matrices by operators [@Honda:2019bdi] assuming $N=\infty$. 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 section \[sec:classicalsolution\], 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 six dimensions. In section \[Sec:analysis\_X\], 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 section \[sec: results\_chiralzeromode\], 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 \[conlusion\] is devoted to a summary and discussions.

Generating classical solutions {#sec:classicalsolution}
==============================

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 [@Ishibashi:1996xs] $$\begin{aligned}
S&=S_\mathrm{b} + S_\mathrm{f}\ , 
\label{S_tot}
\\
S_\mathrm{b} &= -\frac{1}{4g^2} 
{\mbox{Tr}}{\left( {{\left[ A^M, A^N \right]}} {{\left[ A_M, A_N \right]}} \right)} \ , 
\label{S_b}
\\
S_\mathrm{f} &= -\frac{1}{2g^2} 
{\mbox{Tr}}{\left( \Psi_{\alpha} (\mathcal{C} \Gamma^M)_{\alpha\beta} 
{{\left[ A_M, \Psi_{\beta} \right]}} \right)} \ ,
\label{S_f}\end{aligned}$$ where $A^M \ (M=0,\ldots,9)$ and $\Psi_{\alpha}\ (\alpha=1,\ldots,16)$ are $N \times N$ traceless Hermitian matrices, while $\Gamma^M$ and $\mathcal{C}$ are the gamma matrices and the charge conjugation matrix in ten dimensions after the Weyl projection. The indices $M$ and $N$ are raised and lowered by using the metric $\eta={\rm diag}(-1,1, \ldots , 1)$. The action is invariant under the SO(9,1) transformation, where $A^M$ and $\Psi_{\alpha}$ transform as a Lorentz vector and a Majorana-Weyl spinor, respectively. The action is also invariant under the SU$(N)$ transformation $$\label{eq: SUN_transformation}
A_M  \mapsto  U A_M U^\dagger \ , \quad
\Psi_\alpha  \mapsto U \Psi_\alpha U^\dagger \ ,$$ where $U \in \mathrm{SU}(N)$. The Euclidean version can be defined by the replacement $A_0 = i A_{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 $$\label{eq: IRcutoff}
\frac{1}{N}{\mbox{Tr}}(A_0)^2 = 
\kappa \ , \quad \frac{1}{N}{\mbox{Tr}}(A_I)^2 = 1 \quad (I=1,\ldots,9) \ ,$$ corresponding to the IR cutoffs[^1] introduced in Ref. [@Kim:2011cr].

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 $$\label{eq: EOM}
{{\left[ A^M, {{\left[ A_M, A_0 \right]}} \right]}}-\xi A_0 = 0 \ ,\quad
{{\left[ A^M, {{\left[ A_M, A_I \right]}} \right]}}-\zeta A_I = 0 \ ,$$ where $\xi$ and $\zeta$ are the Lagrange multipliers corresponding to the constraints . Note here that if $\xi=\zeta=0$, classical solutions are always simultaneously diagonalizable, as is proved in Appendix A of Ref. [@Steinacker:2017vqw]. Thus the existence of the constraints (\[eq: IRcutoff\]) 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 six dimensions [@Nishimura:2013moa] given as $$\label{eq: ansatz of A}
A_\mu = X_\mu \otimes M \quad (\mu=0,1,2,3) \ , \quad
A_a = \mathbf{1} \otimes Y_a \quad (a=4, \ldots, 9) \ ,$$ where the $N_{\rm X} \times N_{\rm X}$ Hermitian matrices $X_\mu$ represent the structure in the (3+1) dimensions, and the $N_{\rm Y} \times N_{\rm Y}$ Hermitian matrices $Y_a$ represent the structure in the extra six dimensions with $N=N_{\rm X} N_{\rm Y}$. This structure is actually compatible with the SO(3,1) symmetry $$\label{eq: Lorentzsymmetry}
U A_\mu U^\dagger = {O_\mu}^\nu A_\nu\ ,\quad
U A_a U^\dagger = A_a  \ ,$$ where $U= g \otimes \mathbf{1} \in \mathrm{SU}(N)$, $g \in \mathrm{SU}(N_{\rm X})$ and ${O_\mu}^\nu \in \mathrm{SO}(3,1)$ if $X_\mu$ satisfy $g X_\mu g^\dagger = {O_\mu}^\nu X_\nu$. Note that the $N_{\rm Y} \times N_{\rm Y}$ Hermitian matrix $M$ does not have to be $M = \mathbf{1}$, which would correspond to the direct-product structure. We will see that the matrix $M\neq \mathbf{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 (\[eq: EOM\]) requires some method other than just solving differential equations unlike in ordinary classical dynamics. Here we define a “potential” $$V = {\mbox{Tr}}{\left( {{\left[ A^M, {{\left[ A_M, A_0 \right]}} \right]}} - \xi A_0 \right)}^2 
+ {\mbox{Tr}}{\left( {{\left[ A^M, {{\left[ A_M, A_I \right]}} \right]}} - \zeta A_I \right)}^2 \ ,$$ which takes the minimum $V=0$ if and only if $A_M$ solves the classical equation of motion (\[eq: EOM\]). By minimizing $V$ using the gradient descent method starting from random configurations, we can generate classical solutions. Within our ansatz , the constraints reduce to $$\begin{aligned}
\label{eq:normalization}
\frac{1}{N_{\rm X}}\mathrm{Tr} {\left( X_0 \right)}^{2} = \kappa \ , 
\quad \frac{1}{N_{\rm Y}}\mathrm{Tr}M^{2} & =1 \ , \\
\frac{1}{N_{\rm X}}\sum_{i=1}^{3}
\mathrm{Tr}{\left( X_i \right)}^{2}+\frac{1}{N_{\rm Y}}
\sum_{a=4}^{9}\mathrm{Tr}{\left( Y_a \right)}^{2} & =1 \ ,\end{aligned}$$ where (\[eq:normalization\]) is chosen by using the redundancy under the rescaling $X_\mu \mapsto s X_\mu$ and $M \mapsto s^{-1} M$, which leaves $A_\mu = X_\mu \otimes M$ unaltered.

The initial configuration is prepared by generating Gaussian variables for each element of $X_\mu$, $Y_a$ and $M$ in (\[eq: ansatz of A\]). Since the equation of motion (\[eq: EOM\]) is invariant under rescaling $$A_M \to s A_M \ , \quad \xi \to s^2 \xi \ ,
\quad \zeta \to s^2 \zeta \ ,
\label{rescaleing_XYM}$$ we make this rescaling at each step of the algorithm in such a way that the second equation of is satisfied. The algorithm proceeds by updating $X_\mu$, $Y_a$ and $M$ with the increment $$\delta X_\mu = -\epsilon {\frac{{\partial}V}{{\partial}X^{*}_\mu}}\ ,\quad
\delta Y_a = -\epsilon {\frac{{\partial}V}{{\partial}Y^{*}_a}}\ ,\quad
\delta M = -\epsilon {\frac{{\partial}V}{{\partial}M^{*}}}\ ,
\label{updating_XYM}$$ where the stepsize $\epsilon$ is taken to be small enough to make sure that $V$ decreases at each step. We repeat (\[rescaleing\_XYM\]) and (\[updating\_XYM\]) until $V$ becomes as small as $O(10^{-4})$. The first equation of (\[eq:normalization\]) is not imposed on configurations during this procedure, and we use it only to determine the value of $\kappa$ after a classical solution is found. In this way, we can obtain infinitely many solutions with some $\kappa$ depending on the initial random configuration and the initial values of $\xi$ and $\zeta$.

typical solutions {#sec:typical-sol}
-----------------

In typical solutions generated by our algorithm, it turns out that the matrices $M$ and $Y_a$ take the form[^2] $$\label{eq: correspondence_MY}
M= \sqrt{\frac{N_{\rm Y}}{N_{\rm Y} - n_0 }}
\begin{pmatrix}
-\mathbf{1}_{n_{-}} &     & \\
    & \mathbf{0}_{n_0} & \\
    &      & \mathbf{1}_{n_{+}}
\end{pmatrix} \ , \quad
Y_a=
\begin{pmatrix}
Y_a^{(-)} &     & \\
    & Y_a^{(0)} & \\
    &      & Y_a^{(+)}
\end{pmatrix}$$ up to the SU($N_{\rm Y}$) symmetry and $N_{\rm Y}=n_{-} + n_0 + n_{+}$. This implies that $$\begin{aligned}
\frac{N_{\rm Y}-n_0}{N_{\rm Y}} M^3 =M \ , \quad \quad \quad
\label{eq: M} 
{{\left[ M, Y_a \right]}} =0 \ .\end{aligned}$$ Plugging (\[eq: M\]) into the equation of motion (\[eq: EOM\]) with the ansatz , we find that $X_\mu$ and $Y_a$ decouple as $$\begin{aligned}
{{\left[ X^\nu, {{\left[ X_\nu, X_0 \right]}} \right]}} - \xi X_0 &= 0  \ , \nonumber \\
{{\left[ X^\nu, {{\left[ X_\nu, X_i \right]}} \right]}} - \zeta X_i &= 0 \quad (i=1,2,3) \ ,
\label{eq: X_mu}
\\
\label{eq: Y_a}
{{\left[ Y^b, {{\left[ Y_b, Y_a \right]}} \right]}} - \zeta Y_a&= 0  \quad (a=4, \ldots , 9) \ . $$

In this work, we set the parameters $\xi$ and $\zeta$ 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 values of $\xi$ and $\zeta$, it is satisfied automatically at each step. The value of $\kappa$ turns out to be close to $1/3$ as far as we use the same Gaussian distribution for all the elements of $X_\mu$, $Y_a$ and $M$ in the initial configuration. A solution that belongs to the above class (\[eq: correspondence\_MY\]) is obtained without being trapped in a local minimum when the value of the parameter $\zeta$ lies within the range $$\label{eq:mass_region} 
0.1 \lesssim \zeta \lesssim 0.5 \ .$$

Space-time structure in the (3+1) dimensions {#Sec:analysis_X}
============================================

Using the method described in the previous section, we obtain classical solutions of the type IIB matrix model within the ansatz (\[eq: ansatz of A\]). By generating the initial configurations randomly, we can obtain as many solutions as we wish satisfying –. In this section, we focus on the matrices $X_{\mu}$ in (\[eq: ansatz of A\]), which represent the structure of the (3+1)-dimensional space-time.

band-diagonal structure {#SS:Config_X}
-----------------------

Using the SU($N$) symmetry with $U=g \otimes \mathbf{1}$, $g \in \mathrm{SU}(N_{\rm X})$, we can choose a basis in which $X_{0}$ is diagonal with its elements arranged in the ascending order. In Fig. \[fig:X0\], we show the eigenvalues $\alpha_{p}\ (p=1,2,\ldots,N_{\rm X})$ of $X_0$ for a typical solution with $N_{\rm X} = 64$. We find that the eigenvalues are almost uniformly distributed, which suggests that we can define a continuous time in the $N_{\rm X} \to \infty$ limit.

![The eigenvalues $\alpha_{p} \ (p=1,\ldots,N_{\rm X})$ of $X_{0}$ are plotted against the label $p$ in the ascending order for a typical solution with $N_{\rm X}=64$. []{data-label="fig:X0"}](X0.eps){width="10cm"}

![The quantity $\Delta_{pq}$ defined by Eq. (\[eq:ampX\]) is plotted against the labels $p$ and $q$ for the typical solution with $N_{\rm X}=64$. []{data-label="fig:ampX"}](M_2.eps){width="10cm"}

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 $$\label{eq:ampX}
\Delta_{pq} 
=\sum_{i=1}^{3}\left|(X_{i})_{pq}\right|^{2}$$ for the same solution in Fig. \[fig:ampX\]. We find that $\Delta_{pq}$ becomes very small for $|p-q|\geq n$ with some integer $n$, which is $n \sim 10$ for the configuration in Fig. \[fig:ampX\]. 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 $\alpha_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[@Ito:2013ywa; @Ito:2015mxa; @Ito:2017rcr; @Azuma:2017dcb]. As we mentioned below (\[eq: EOM\]), the classical solutions are strictly diagonal for $\xi=\zeta=0$. Therefore, the band-diagonal structure may be viewed as a deviation from the diagonal matrix due to finite $\xi$ and $\zeta$.

extracting the time evolution {#SS:Time_evolution}
-----------------------------

Next we discuss how we extract the time evolution from a given classical solution. The band-diagonal structure discussed in the previous section plays a crucial role here.

Let us define an $n\times n$ block $\bar{X}_i(t)$ by $$(\bar{X}_i (t))_{rs}
= (X_{i})_{k+r, k+s}
\quad (r, s=1,2,\ldots,n) \ ,$$ where $k=0,1,\ldots,N_{\rm X}-n$ and the argument $t$ represents the time defined by $$t
=\frac{1}{n}\sum_{r=1}^{n}\alpha_{k+r}  \ .$$ This block matrix $\bar{X}_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 $$R^{2}\left(t\right)
=\frac{1}{n}\sum_{i=1}^{3}
\mathrm{tr}\left(\bar{X}_{i} (t)\right)^2 \ ,$$ where the symbol “${\mbox{tr}}$” represents the trace over an $n \times n$ matrix. In Fig. \[fig:R2\], we plot $R^{2}(t)$ for the typical solution discussed in section \[SS:Config\_X\]. We find that $R^{2}(t)$ is roughly symmetric under $t\mapsto -t$ reflecting the symmetry of the model under $A_0 \mapsto - A_0$. In particular, an expanding behavior is seen at $t>0$, which is consistent with the results of simulations [@Nishimura:2019qal] at late time.

![The extent of space $R^{2}(t)$ is plotted against $t$ for the typical solution with $N_{\rm X}=64$ and the block size $n=10$. []{data-label="fig:R2"}](R2.eps){width="10cm"}

Next we define a $3\times3$ real symmetric tensor $$T_{ij}\left(t\right)=\frac{1}{n}\mathrm{tr}
\left(\bar{X}_{i}\left(t\right)\bar{X}_{j}\left(t\right)\right) \ ,$$ whose eigenvalues $\lambda_i(t)$, which we order as $\lambda_{1}(t) \leq \lambda_{2}(t) \leq \lambda_{3}(t)$, represent how the 3-dimensional space extends in each direction. Note that they are related to the extent of space $R^2(t)$ as $$R^{2}\left(t\right)=\sum_{i=1}^{3} \lambda_{i} (t)  \ .$$ In Fig. \[fig:T\], we plot the eigenvalues $\lambda_i(t)$ for the typical solution discussed above. The three eigenvalues $\lambda_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_{\rm X}$, suggesting that the SO(3) symmetry is recovered in the limit of infinite matrix size.

![The eigenvalues $\lambda_i (t)$ of $T_{ij}(t) \ (i,j=1,2,3)$ are plotted against $t$ for the typical solution with $N_{\rm X}=64$ and the block size $n=10$. []{data-label="fig:T"}](EV_Tij.eps){width="10cm"}

![The eigenvalues of $Q(t)$ are plotted against $t$ for the typical solution with $N_{\rm X}=64$. There are ten of them reflecting the chosen block size $n=10$. []{data-label="fig:Q"}](EV_Q.eps){width="10cm"}

In order to probe the structure of the 3-dimensional space emerging from the classical solutions in more detail, we define an $n \times n$ matrix $$Q\left(t\right)=\sum_{i=1}^{3} \bar{X}_{i}(t)^2  \ ,$$ whose eigenvalues $q_{r}(t)$ represent the radial distribution of the points which describe the 3-dimensional space. Note that the eigenvalues $q_{r}(t)$ are related to the extent of space $R^{2}(t)$ as $$R^{2} (t) 
=\frac{1}{n}\sum_{r=1}^{n} q_{r} (t) \ .$$ In Fig. \[fig:Q\], 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 [@Aoki:2019tby], 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 {#sec: results_chiralzeromode}
======================================

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 (\[eq: correspondence\_MY\]) 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.

the Dirac equation
------------------

From the fermionic action (\[S\_f\]) of the type IIB matrix model, one obtains the 10-dimensional Dirac equation $$\label{eq: 10d_Dirac}
\Gamma^M {{\left[ A_M, \Psi \right]}}=0 \ ,$$ where $\Psi$ satisfies the Weyl condition $$\label{eq: chirality_in_10d}
\Gamma_\chi \Psi = \Psi$$ with the chirality operator $\Gamma_\chi$ in 10 dimensions. Here we decompose the 10-dimensional gamma matrices $\Gamma^M$ as $$\Gamma^\mu= \rho^\mu \otimes \mathbf{1} \ , \quad
\Gamma^a= \rho_\chi \otimes \gamma^a \ ,$$ where $\rho^\mu$ and $\gamma^a$ are the 4-dimensional and 6-dimensional gamma matrices, respectively. Note that the chirality operator $\Gamma_\chi$ is also decomposed as $$\Gamma_\chi = \rho_\chi \otimes \gamma_\chi \ ,$$ where $\rho_\chi$ and $\gamma_\chi$ are the chirality operators in 4 and 6 dimensions, respectively.

Suppose $\Psi$ is chiral in four dimensions meaning $$\label{eq: chirality_in_4d_Psi}
(\rho_\chi \otimes \mathbf{1}) \Psi = \pm \Psi \ .$$ Due to the chirality (\[eq: chirality\_in\_10d\]) in 10 dimensions, Eq. (\[eq: chirality\_in\_4d\_Psi\]) implies $$\label{eq: chirality_in_6d_Psi}
{\left( \mathbf{1} \otimes \gamma_\chi \right)} \Psi = \pm \Psi \ ;$$ i.e., $\Psi$ is chiral also in six dimensions.

Let us next decompose Eq. (\[eq: 10d\_Dirac\]) into two terms as $$\Gamma^\mu {{\left[ A_\mu, \Psi \right]}}+ \Gamma^a {{\left[ A_a, \Psi \right]}}=0 \ .
\label{eq:ac-decompose}$$ 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 six dimensions; i.e., $$\Gamma^a {{\left[ A_a, \Psi \right]}}=0 \ .
\label{Dirac-eq-6d}$$

Let us now consider that $A_M$ is a classical solution with the quasi direct-product structure (\[eq: ansatz of A\]). Since $A_a = \mathbf{1} \otimes Y_a$, the general solution to (\[Dirac-eq-6d\]) can be obtained by decomposing $\Psi$ as $$\Psi = \psi^\mathrm{(4d)} \otimes \psi^\mathrm{(6d)} \ ,$$ where the 4-dimensional and 6-dimensional gamma matrices act only on $\psi^\mathrm{(4d)}$ and $\psi^\mathrm{(6d)}$, respectively. Thus, in order to satisfy (\[Dirac-eq-6d\]), we only need to require $$\label{eq: chirality_in_6d_psi}
\gamma^a {{\left[ Y_a, \psi_a^\mathrm{(6d)} \right]}}=0 \ .$$

Dirac zero modes from intersecting D-branes
-------------------------------------------

When the configuration of $Y_a$ takes the block-diagonal form (\[eq: correspondence\_MY\]), the eigenvalue problem associated with the Dirac operator on the left-hand side of (\[eq: chirality\_in\_6d\_psi\]) reduces to that for each block matrix of $\psi^\mathrm{(6d)}$. Note that the diagonal blocks of $Y_a$ can be interpreted as D-branes, whereas an off-diagonal block of $\psi^\mathrm{(6d)}$ can be interpreted as the open string ending on the corresponding two D-branes. When the two D-branes intersect with each other at a point, the length of the open string can become zero as it comes close to the intersecting point, and therefore a chiral zero mode appears there. In what follows, we confirm that this indeed happens for the typical solutions.

For that purpose, it suffices to consider the solutions of with a block-diagonal structure $$\begin{aligned}
\label{eq: form_Ya}
Y_a &=
\begin{pmatrix}
Y_a^{(1)} & \\
 & Y_a^{(2)}
\end{pmatrix} \ ,\end{aligned}$$ where the block matrices $Y_a^{(1)}$ and $Y_a^{(2)}$ have to satisfy with the same $\zeta$, which we set to $\zeta=1$ without loss of generality. We denote the size of the block matrices $Y_a^{(1)}$ and $Y_a^{(2)}$ as $N_{\rm Y}^{(1)}$ and $N_{\rm Y}^{(2)}$, respectively.

Corresponding to the block-diagonal structure (\[eq: form\_Ya\]) of $Y_a$, we focus on the off-diagonal block $\varphi$ in $$\begin{aligned}
\psi^\mathrm{(6d)}&=
\begin{pmatrix}
 & \varphi \\
 & 
\end{pmatrix} \ ,
\label{eq: off-diagonal_block_psi}\end{aligned}$$ for which the eigenvalue problem reduces to $$\label{eq: Diraceq_6d}
\gamma^a_{\alpha \beta} {\left[ Y_a^{(1)} \varphi_\beta
- \varphi_\beta Y_a^{(2)} \right]} =\lambda \varphi_\alpha \ ,$$ where $\varphi_\alpha \ (\alpha=1,\ldots,8)$ are the eigenvectors. The eigenvalue $\lambda$ corresponds to the mass in the (3+1) dimensions as one can see from Eq. (\[eq:ac-decompose\]).

We obtain numerically the eigenvalue spectrum of the Dirac operator for various pairs of classical solutions $Y_a^{(1)}$ and $Y_a^{(2)}$ 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 $\varphi_\alpha$, 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 $\gamma_\chi$, the eigenvalues $\lambda$ and $- \lambda$ appear in pairs. If $\varphi$ is an eigenvector corresponding to the eigenvalue $\lambda$, then $\gamma_{\chi}\varphi$ is an eigenvector corresponding to the eigenvalue $- \lambda$. For the zero modes corresponding to $\lambda=0$, one can choose them to have definite chirality $\gamma_\chi \varphi = \pm \varphi$. 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 $$\varphi_{\mathrm{R} \alpha} = 
\frac{1+ \gamma_\chi}{2} \varphi_\alpha \ , \quad
\varphi_{\mathrm{L} \alpha} = 
\frac{1- \gamma_\chi}{2} \varphi_\alpha \ .
\label{eq: unitary_transf_psi}$$

Furthermore, we should take into account that the classical solution of the form (\[eq: form\_Ya\]) has symmetry under $$\label{eq: UEX}
Y_a  \mapsto  g \, Y_a \, g^\dagger \ ,$$ where the $g \in \mathrm{SU}(N_{\rm Y})$ is block-diagonal $$g =
\begin{pmatrix}
g^{(1)} & 0 \\
0 & g^{(2)}
\end{pmatrix}$$ with $g^{(1)} \in \mathrm{SU}(N_{\rm Y}^{(1)})$ and $g^{(2)} \in \mathrm{SU}(N_{\rm Y}^{(2)})$, which preserves the structure . Under this transformation, the wave functions of the corresponding Dirac operator transform as $$\label{eq: varphi_LR_prime}
\varphi_\alpha \mapsto 
g^{(1)} \varphi_\alpha g^{(2)\dagger}  \ .
$$ We choose $g^{(1)}$ and $g^{(2)}$ in in such a way that $\varphi_{\mathrm{R}1}$ is diagonal with positive elements in the descending order. This is nothing but the singular value decomposition of $\varphi_{\mathrm{R}1}$, where the diagonal elements correspond to the singular values.

In what follows, we change the index $a$ of $Y_a$ from $a=4,\ldots,9$ to $a=1,\ldots,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.

We generate 3d and 4d solutions numerically using the algorithm described in section \[sec:classicalsolution\]. The 2d solutions, on the other hand, can be constructed analytically by defining $Z$ as $$\label{eq: Z}
{{\left[ Y_1, Y_2 \right]}}=i Z \ ,$$ so that we can reduce with $\zeta=1$ to $$\label{eq: Y1Z_Y2Z}
{{\left[ Y_2, Z \right]}} = i  Y_1  \ , \quad
{{\left[ Z, Y_1 \right]}}= i  Y_2 \ .$$ Eqs.  and imply $$Y_1 = L_1 \ , \quad Y_2 = L_2 \ , \quad Z = L_3 \ ,
\label{2d-solution-generator}$$ 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 onto 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
--------------------------------------

#### “3d-3d” ansatz {#subsec: 3d-3d}

For the first brane $Y_a^{(1)}$, we set $Y_4^{(1)} = Y_5^{(1)} = Y_6^{(1)} = 0$, and for the second brane $Y_a^{(2)}$, we set $Y_1^{(2)} = Y_2^{(2)} = Y_3^{(2)} = 0$. We use four 3d solutions for each of them. Thus we consider $4\times 4=16$ backgrounds with the matrix size $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $40$, $48$, $56$ and $64$, on which we solve the eigenvalue problem (\[eq: Diraceq\_6d\]).

![The ratio $\mu_0/\mu_1$ with the “3d-3d” ansatz is plotted against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $40$, $48$, $56$ and $64$. The dashed line represents a fit to the quadratic function $a+b/N_{\rm Y}^{(1)}+c/(N_{\rm Y}^{(1)})^2$ of $1/N_{\rm Y}^{(1)}$ with $a=-0.04(7)$, $b=34(6)$ and $c=-5(1) \times 10^2$. []{data-label="ratio_EV_3d"}](ratio_3d-3d.eps){width="10.0cm"}

![The wave functions $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $|(\varphi_{\mathrm{L}5})_{pq}|^2$ $(p=1,\ldots,N_{\rm Y}^{(1)},  \ q=1,\ldots,N_{\rm Y}^{(2)})$ with the “3d-3d” ansatz are plotted for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. []{data-label="wavefn_3d-3d"}](wavefn_3d-3d_2.eps){width="16.5cm"}

We obtain the smallest absolute value $|\lambda|$ of the eigenvalues for each of the 16 cases and denote the average as $\mu_0$. Similarly, we denote as $\mu_1$ the average of the second smallest absolute value $|\lambda|$. Since the normalization of the solutions is somewhat arbitrary for different matrix size, we take the ratio $\mu_0/\mu_1$ and plot it against $1/N_{\rm Y}^{(1)}$ in Fig. \[ratio\_EV\_3d\]. The fit to the quadratic function of $1/N_{\rm Y}^{(1)}$ suggests that zero modes appear in the $N_{\rm Y}^{(1)} \to \infty$ limit.

Next we consider the wave function corresponding to the lowest eigenvalue for one of the 16 cases with $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. We calculate $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ for each $\alpha$ and find that the wave function almost vanishes except for the spinor component $\alpha=1$. Similarly, we calculate $|(\varphi_{\mathrm{L} \alpha})_{pq}|^2$ for each $\alpha$, and find that the wave function almost vanishes except for the spinor component $\alpha=5$. In Fig. \[wavefn\_3d-3d\], we plot $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $[(\varphi_{\mathrm{L}5})_{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.

#### “2d-4d” ansatz {#subsec: 2d-4d}

For the first brane $Y_a^{(1)}$, we set $Y_3^{(1)} =Y_4^{(1)} = Y_5^{(1)} = Y_6^{(1)} = 0$ and use the analytic solution corresponding to the irreducible representation of the SU(2) algebra with $N_{\rm Y}^{(1)}=24$, $28$, $32$ and $36$. For the second brane $Y_a^{(2)}$, we set $Y_1^{(2)} = Y_2^{(2)} = 0$ and use sixteen 4d solutions obtained[^3] with $N_{\rm Y}^{(2)}=(N_{\rm Y}^{(1)})^2/16$. Thus we consider $16$ backgrounds, on which we solve the eigenvalue problem (\[eq: Diraceq\_6d\]).

In the present case, the eigenvalues turn out to have two-fold degeneracy, which may be understood from the fact that the 2d brane is actually something like two coinciding fuzzy disks. (See Eq. (\[2d-solution-generator\]) and the lines below.) In Fig. \[ratio\_EV\_2d-4d\_NY1=\[24-28-32-36\]\], we plot the ratio $\mu_0/\mu_1$ against $1/N_{\rm Y}^{(1)}$. The fit to the quadratic function of $1/N_{\rm Y}^{(1)}$ suggests that zero modes appear in the $N_{\rm Y}^{(1)} \to \infty$ limit.

![The ratio $\mu_0/\mu_1$ with the “2d-4d” ansatz is plotted against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=24$, $28$, $32$ and $36$ with $N_{\rm Y}^{(2)}=(N_{\rm Y}^{(1)})^2/16$. The dashed line is a fit to $a+ b/N_{\rm Y}^{(1)}+ c/(N_{\rm Y}^{(1)})^2$ with $a=-0.01(16)$, $b=13(9)$ and $c=-1.3(1.3)\times 10^2$.[]{data-label="ratio_EV_2d-4d_NY1=[24-28-32-36]"}](ratio_2d-4d.eps){width="10.0cm"}

Next we consider the wave function corresponding to one of the 2 lowest eigenvalues[^4] for one of the 16 cases with $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. We calculate $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ for each $\alpha$ and find that the wave function almost vanishes except for the spinor component $\alpha=1$. Similarly, we calculate $|(\varphi_{\mathrm{L}\alpha})_{pq}|^2$ for each $\alpha$, and find that the wave function almost vanishes except for the spinor component $\alpha=5$. In Fig. \[wavefn\_2d-4d\], we plot $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $|(\varphi_{\mathrm{L}5})_{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 wave functions $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $|(\varphi_{\mathrm{L}5})_{pq}|^2$ $(p=1,\ldots,N_{\rm Y}^{(1)},\ q=1,\ldots,N_{\rm Y}^{(2)})$ with the “2d-4d” ansatz are plotted for $N_{\rm Y}^{(1)}=36$ and $N_{\rm Y}^{(2)}=81$. []{data-label="wavefn_2d-4d"}](wavefn_2d-4d.eps){width="16.5cm"}

the cases without intersection
------------------------------

#### “2d-3d” ansatz {#subsec: 2d-3d}

For the first brane $Y_a^{(1)}$, we set $Y_3^{(1)} =Y_4^{(1)} = Y_5^{(1)} = Y_6^{(1)} = 0$ and use the analytic solution corresponding to the irreducible representation of the SU(2) algebra with $N_{\rm Y}^{(1)}=32$, $48$ and $64$. For the second brane $Y_a^{(2)}$, we set $Y_1^{(2)} = Y_2^{(2)} =Y_3^{(2)} = 0$ and use sixteen 3d solutions obtained with $N_{\rm Y}^{(2)}=N_{\rm Y}^{(1)}$. Thus we consider $16$ backgrounds, on which we solve the eigenvalue problem (\[eq: Diraceq\_6d\]).

As in the “2d-4d” case discussed above, we have two-fold degeneracy in the eigenvalue spectrum. In Fig. \[ratio\_EV\_2d-3d\_NY1=\[32-48-64\]\], we plot the ratio $\mu_0/\mu_1$ against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $48$ and $64$. We find that it does not converge to $0$ in the $N_{\rm Y}^{(1)} \to \infty$ limit, which implies that we do not obtain zero modes in this case.

![The ratio $\mu_0/\mu_1$ with the “2d-3d” ansatz is plotted against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $48$ and $64$.[]{data-label="ratio_EV_2d-3d_NY1=[32-48-64]"}](ratio_2d-3d.eps){width="10.0cm"}

![The spinor components $\alpha = 1,2,5,6$ of the wave functions $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ and $|(\varphi_{\mathrm{L}\alpha})_{pq}|^2$ $(p=1,\ldots,N_{\rm Y}^{(1)},\ q=1,\ldots,N_{\rm Y}^{(2)})$ with the “2d-3d” ansatz are plotted for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. []{data-label="wavefn_2d-3d"}](wavefn_2d-3d_2.eps){width="16.5cm"}

Next we consider the wave function corresponding to one of the 2 lowest eigenvalues for one of the 16 cases with $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. We calculate $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ for each $\alpha$ and find that the wave function almost vanishes except for the spinor components $\alpha=1,2,5,6$. Similarly, we calculate $|(\varphi_{\mathrm{L}\alpha})_{pq}|^2$ for each $\alpha$, and find that the wave function almost vanishes except for the spinor components $\alpha=1,2,5,6$. In Fig. \[wavefn\_2d-3d\], we plot $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ and $|(\varphi_{\mathrm{L}\alpha})_{pq}|^2$ for $\alpha=1,2,5,6$, the components of the wave functions that have nonzero elements. We find that $\varphi_{\mathrm{R}\alpha}$ is localized at the $(1,1)$ element for $\alpha=1,2$ and at the $(1,2)$ element for $\alpha=5,6$, while $\varphi_{\mathrm{L}\alpha}$ is localized at the $(1,2)$ element for $\alpha=1,2$ and at the $(1,1)$ element for $\alpha=5,6$. These results are consistent with the picture that the 2d brane and the 3d brane do not intersect, but the points on each brane separated by the minimum distance are uniquely determined.

the cases with intersection but not at a point
----------------------------------------------

#### “3d-4d” ansatz {#subsec: 3d-4d}

For the first brane $Y_a^{(1)}$, we set $Y_4^{(1)} = Y_5^{(1)} = Y_6^{(1)} = 0$ and use four 3d solutions, while for the second brane $Y_a^{(2)}$, we set $Y_1^{(2)} = Y_2^{(2)} = 0$ and use four 4d solutions. Thus we consider $4\times 4=16$ backgrounds with the matrix size $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $48$ and $64$, on which we solve the eigenvalue problem (\[eq: Diraceq\_6d\]).

In Fig. \[ratio\_EV\_3d-4d\_NY1=\[32-48-64\]\], we plot the ratio $\mu_0/\mu_1$ against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $48$ and $64$. They do not converge to $0$ in the $N_{\rm Y}^{(1)} \to \infty$ limit, which implies that we do not obtain zero modes in this case.

![The ratio $\mu_0/\mu_1$ with the “3d-4d” ansatz is plotted against $1/N_{\rm Y}^{(1)}$ for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=32$, $48$ and $64$.[]{data-label="ratio_EV_3d-4d_NY1=[32-48-64]"}](ratio_3d-4d.eps){width="10.0cm"}

![The wave functions $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $|(\varphi_{\mathrm{L}5})_{pq}|^2$ $(p=1,\ldots,N_{\rm Y}^{(1)},\ q=1,\ldots,N_{\rm Y}^{(2)})$ with the “3d-4d” ansatz are plotted for $N_{\rm Y}^{(1)}=N_{\rm Y}^{(2)}=64$. []{data-label="wavefn_3d-4d"}](wavefn_3d-4d_2.eps){width="16.5cm"}

Next we consider the wave function corresponding to the lowest eigenvalue for one of the 16 cases with $N_{\rm Y}^{(1)}=64$. We calculate $|(\varphi_{\mathrm{R}\alpha})_{pq}|^2$ for each $\alpha$ and find that the wave function almost vanishes except for the spinor component $\alpha=1$. Similarly, we calculate $|(\varphi_{\mathrm{L}\alpha})_{pq}|^2$ for each $\alpha$, and find that the wave function almost vanishes except for the spinor component $\alpha=5$. In Fig. \[wavefn\_3d-4d\], we plot $|(\varphi_{\mathrm{R}1})_{pq}|^2$ and $|(\varphi_{\mathrm{L}5})_{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 {#conlusion}
=======================

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 direct-product structure , 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 six dimensions.

First, we have focused on the (3+1)-dimensional space-time structure, which is represented by $X_\mu$ in . 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 [@Aoki:2019tby]. Our results support the conjecture in Ref. [@Nishimura:2019qal] that the singular structure is caused by the approximation and that a smooth space-time should emerge in the $N \to \infty$ 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 \neq {\bf 1}$ in the quasi direct-product structure 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 [@Aoki:2010gv; @Chatzistavrakidis:2011gs; @Nishimura:2013moa; @Aoki:2014cya], 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_\mu$ with larger $N_{\rm 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 3-dimensional space becomes uniform and SO(3) symmetric. By obtaining $Y_a$ with larger $N_{\rm 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.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors would like to thank H. Kawai and H. Steinacker for valuable discussions. Computation was carried out on PC clusters at KEK and XC40 at YITP in Kyoto University. A.T. was supported in part by Grant-in-Aid for Scientific Research (No. 18K03614) from Japan Society for the Promotion of Science. K.H. is supported by Grant-in-Aid for JSPS Fellows (No. 19J10002).

[^1]: Strictly speaking, the IR cutoffs introduced in Ref. [@Kim:2011cr] take the form of inequalities unlike (\[eq: IRcutoff\]). This is not a big difference, however, since the inequalities are actually saturated in numerical simulations for entropic reasons.

[^2]: The coefficient $\frac{N_{\rm Y}}{N_{\rm Y} - n_0}$ for $M$ is simply due to the chosen normalization $\frac{1}{N_{\rm Y}}\mathrm{Tr}M^{2} =1$.

[^3]: The chosen matrix size $N_{\rm Y}^{(2)}$ for the 4d brane is motivated from the fact that the degrees of freedom on a lattice with a linear extent $L$ grow as $L^2$ and $L^4$ for 2d and 4d cases, respectively. The factor of $1/16$ in $N_{\rm Y}^{(2)}=(N_{\rm Y}^{(1)})^2/16$ is introduced to avoid having too large $N_{\rm Y}^{(2)}$ to perform explicit calculations.

[^4]: The situation with the other lowest eigenvalue is qualitatively the same. The same comment applies also to the case with the “2d-3d” ansatz below.
