On the structure of the emergent 3d expanding space in the Lorentzian type IIB matrix model

TL;DR

This paper investigates the structure of emergent 3D expanding space in the Lorentzian type IIB matrix model, revealing that the expansion is linked to Pauli matrices and discussing implications of simulation approximations.

Contribution

It analyzes the space-time structure in the Lorentzian type IIB matrix model and its simplified versions, highlighting the impact of simulation approximations on the emergent geometry.

Findings

Expanding space is described by Pauli matrices in the simplified models.

Simulation approximation replaces $e^{iS_b}$ with $e^{eta S_b}$, affecting results.

Discussion on how to recover the original model's space-time structure.

Abstract

The emergence of (3+1)-dimensional expanding space-time in the Lorentzian type IIB matrix model is an intriguing phenomenon which was observed in Monte Carlo studies of this model. In particular, this may be taken as a support to the conjecture that the model is a nonperturbative formulation of superstring theory in (9+1) dimensions. In this paper we investigate the space-time structure of the matrices generated by simulating this model and its simplified versions, and find that the expanding part of the space is described essentially by the Pauli matrices. We argue that this is due to an approximation used in the simulation to avoid the sign problem, which actually amounts to replacing $e^{iS_{\rm b}}$ by $e^{\beta S_{\rm b}}$ ($\beta>0$) in the partition function, where $S_{\rm b}$ is the bosonic part of the action. We also discuss the possibility of obtaining a regular space-time with the (3+1)-dimensional expanding behavior in the original model with the correct $e^{iS_{\rm b}}$ factor.

Tables

Table 1. Table 1: The parameter sets ( N 𝑁 N , C 𝐶 C , κ 𝜅 \kappa ) used for the simulation of the VDM model are listed. We also present the block size n 𝑛 n , the “volume” Δ Δ \Delta and the “lattice spacing” ϵ italic-ϵ \epsilon determined from the data for each parameter set.

$N$	$C$	$\kappa$	$n$	$\Delta$	$ϵ$
64	8.81	0.14	24	1.0990(16)	0.0550(1)
96	0	2.00	14	1.3811(41)	0.1151(3)
64	0	2.00	10	1.2726(63)	0.1591(8)
64	0	4.00	7	1.3762(87)	0.2752(17)

Equations

$$S$$

$$S_{\mathrm{b}}$$

$$S_{\mathrm{f}}$$

$$Z=\int dAd\Psi \text{mbox}{e}^{iS\left[A,\Psi \right]}=\int dA\mathrm{Pf}\mathcal{M}(A)\text{mbox}{e}^{i{S}_{\mathrm{b}}},$$

$$S_{\mathrm{b}}=\frac{1}{4}\mathrm{Tr}\left(F_{\mu \nu }{F}^{\mu \nu }\right)=\frac{1}{4}\left\{-2\mathrm{Tr}{\left(F_{0i}\right)}^2+\mathrm{Tr}{\left(F_{ij}\right)}^2\right\},$$

$$\frac{1}{N}\mathrm{Tr}{\left\{{\left(A_0\right)}^2\right\}}^p$$

$$\frac{1}{N}\mathrm{Tr}{\left\{{\left(A_i\right)}^2\right\}}^p$$

$$Z$$

$$A_0=\mathrm{diag}({\alpha }_1,\dots ,{\alpha }_N),\text{mbox}where{\alpha }_1<\cdots <{\alpha }_N.$$

$${\left({\bar{A}}_i\right)}_{IJ}(t)\equiv {\left(A_i\right)}_{\nu +I,\nu +J}$$

$$t\equiv \frac{1}{n}\sum _{I=1}^n{\alpha }_{\nu +I},$$

$$R^2(t)=⟨\sum _{i=1}^9\frac{1}{n}\mathrm{tr}{\left({\bar{A}}_i(t)\right)}^2⟩,$$

$$T_{ij}(t)=\frac{1}{n}\mathrm{tr}\left({\bar{A}}_i(t){\bar{A}}_j(t)\right),$$

$${\lambda }_1(t)>{\lambda }_2(t)>\cdots >{\lambda }_9(t).$$

$$R^2(t)=⟨\mathrm{tr}T⟩=\sum _{i=1}^9⟨{\lambda }_i(t)⟩.$$

$$\mathrm{Pf}\mathcal{M}(A)⟹\{\begin{array}{ll}\Delta (\alpha {)}^{16}& \text{for the VDM model},\\ 1& \text{for the bosonic model},\end{array}$$

$$R(t)\sim \text{mbox}{e}^{\Lambda t},$$

$$R(t)\sim t^{1/2}$$

$$Q(t)\equiv \sum _{i=1}^9({\bar{A}}_i(t){)}^2,$$

$$q_1(t)<\cdots <q_n(t).$$

$$R^2(t)=⟨\frac{1}{n}\mathrm{tr}Q(t)⟩=⟨\frac{1}{n}\sum _{k=1}^n{q}_k(t)⟩.$$

$${\bar{A}}^{(i)}(t)=\sum _{j=1}^9{v}_j^{(i)}(t){\bar{A}}_j(t)$$

$$a_1^{(i)}(t)<\cdots <a_n^{(i)}(t).$$

$$\lambda_{i}(t)=\frac{1}{n}\sum_{k=1}^{n}\Big{(}a^{(i)}_{k}(t)\Big{)}^{2}\ .$$

$$\Delta \equiv \frac{t_{\mathrm{peak}}-t_{\mathrm{c}}}{R\left(t_{\mathrm{c}}\right)},ϵ\equiv \frac{\Delta }{\nu },$$

$$\bar{A}^{(i)}(t)=\left(\begin{array}[]{cc}\ast&\ast\\ \ast&X_{i}(t)\end{array}\right)$$

$$[X_i,X_j]=ic{ϵ}_{ijk}{X}_k$$

$$S(c)\equiv \mathrm{tr}(i{ϵ}_{ijk}[X_i,X_j]+2c{X}_k{)}^2,$$

$$\tilde{c}=-\frac{i{ϵ}_{ijk}\mathrm{tr}(X_k[X_i,X_j])}{2\mathrm{tr}(X_l^2)}.$$

$$Z$$

$$Z=$$

$$f(S_{\mathrm{b}})\equiv {\int }_0^{L^{2p}}d(r^p)r^{9(N^2-1)-1}\text{mbox}{e}^{r^2(i{S}_{\mathrm{b}}-ϵ|S_{\mathrm{b}}|)},$$

Full text

KEK-TH-2110

On the structure of the emergent 3d expanding space

in the Lorentzian type IIB matrix model

Toshihiro Aokia111 E-mail address : toshaoki@post.kek.jp, Mitsuaki Hirasawaa222 E-mail address : mitsuaki@post.kek.jp, Yuta Itob333 E-mail address : yito@post.kek.jp,

Jun Nishimuraab444 E-mail address : jnishi@post.kek.jp and Asato Tsuchiyac555 E-mail address : tsuchiya.asato@shizuoka.ac.jp

*a**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*

*b**KEK Theory Center, High Energy Accelerator Research Organization,

1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan*

*c**Department of Physics, Shizuoka University,

836 Ohya, Suruga-ku, Shizuoka 422-8529, Japan*

abstract

The emergence of (3+1)-dimensional expanding space-time in the Lorentzian type IIB matrix model is an intriguing phenomenon which was observed in Monte Carlo studies of this model. In particular, this may be taken as a support to the conjecture that the model is a nonperturbative formulation of superstring theory in (9+1) dimensions. In this paper we investigate the space-time structure of the matrices generated by simulating this model and its simplified versions, and find that the expanding part of the space is described essentially by the Pauli matrices. We argue that this is due to an approximation used in the simulation to avoid the sign problem, which actually amounts to replacing $\mbox{e}^{iS_{\rm b}}$ by $\mbox{e}^{\beta S_{\rm b}}$ ($\beta>0$) in the partition function, where $S_{\rm b}$ is the bosonic part of the action. We also discuss the possibility of obtaining a regular space-time with the (3+1)-dimensional expanding behavior in the original model with the correct $\mbox{e}^{iS_{\rm b}}$ factor.

1 Introduction

Superstring theory is a promising candidate for a unified theory that includes quantum gravity consistently. One of the striking consequences of this theory is that the space-time should have ten dimensions. Therefore, in order to make the theory compatible with our (3+1)d world, the extra six dimensions have to be compactified somehow. Depending on the structure of these compact extra dimensions, one can obtain various quantum field theories in the (3+1)d space-time at low energy. This issue has been investigated extensively at the perturbative level including D-branes configurations as a background accounting for certain nonperturbative effects, and it led to tremendously many consistent vacua, the situation which is called the landscape. However, it remains to be seen what really happens if one formulates the theory in a fully nonperturbative manner as one does in the case of quantum field theory using the lattice formulation.

The type IIB matrix model [1] was proposed as such a nonperturbative formulation of superstring theory. Formally, the action of the model can be obtained by dimensionally reducing the action of 10d $\mathcal{N}=1$ SYM theory to 0d, and it actually has maximal $\mathcal{N}=2$ supersymmetry in 10d. The space-time does not exist a priori, and it is represented by the eigenvalue distribution of the ten bosonic matrices $A_{\mu}$ $(\mu=0,\cdots,9)$. This is manifested by the fact that translations in the supersymmetry algebra turn out to be realized by the shifts $A_{\mu}\mapsto A_{\mu}+\alpha_{\mu}{\bf 1}$ in the ten bosonic matrices. The model, therefore, has the potential to clarify a possible nonperturbative mechanism for dynamical compactification in superstring theory. The Euclidean version of the model was investigated from this viewpoint, and the spontaneous breaking of SO(10) rotational symmetry was suggested by various approaches [2, 3, 4, 5, 6, 7]. However, latest calculations based on the Gaussian expansion method suggested that the SO(10) symmetry is broken down to SO(3) instead of SO(4) [5].

This provided a strong motivation to consider the Lorentzian version of the model. Monte Carlo simulation was performed in ref. [8], and the results turned out to be intriguing. In the SU($N$) basis which diagonalizes the temporal matrix $A_{0}$, the spatial matrices $A_{i}$ $(i=1,\cdots,9)$ have a band-diagonal structure, which enabled the extraction of the real-time evolution. In this way, it was found that only three out of nine directions start to expand at some critical time, which implies that the model predicts the emergence of a (3+1)d expanding space-time from superstring theory in (9+1)d. The expanding behavior for a longer time was investigated by simulating simplified models. The obtained results suggested a scenario for the full model that the expansion is exponential at early times [9], which is reminiscent of the inflation, and that it turns into a power law [10] at later times, which is reminiscent of the Friedmann-Robertson-Walker universe in the radiation dominated era. See also refs. [11, 12, 13] for closely related work.

Due to the expansion of space, it is expected that the dominant configurations can be well approximated at late times by some classical solution of the Lorentzian type IIB matrix model. Indeed several types of classical solutions representing expanding space-time have been constructed [14, 15, 16, 17, 18, 19, 20, 21]. Also, matrix configurations with various structures in the extra dimensions are considered to realize chiral fermions in the (3+1)d space-time. Earlier attempts used slightly modified models obtained, for instance, by orbifolding [22, 23] or by toroidal compactification with magnetic fluxes [24]. More recently [25, 26], it was shown that the original model can be used to realize the idea of intersecting D-branes [27], which led to the proposal of matrix configurations that can give rise to phenomenologically viable low-energy effective theories [28, 29, 30].

In this paper we investigate the space-time structure of the matrix configurations generated by Monte Carlo simulation of the Lorentzian type IIB matrix model and the simplified models. In particular, we calculate the eigenvalues of the submatrices of the spatial matrices $A_{i}$ corresponding to each time slice and find that only two eigenvalues grow in magnitude after the critical time. A more detailed analysis shows that the expanding 3d space is described essentially by the Pauli matrices. Namely the space is actually more like a fuzzy sphere although it has been called “3d” in the sense that it has three extended directions. While we keep on using the word 3d in this sense in what follows, we refer to the space with the Pauli-matrix structure as a “singular 3d space”. We observe that the situation remains unaltered even at late times or in the continuum limit for the simplified models, and it is shared by the original model with maximal supersymmetry as well. This raises the important question of whether this model can generate a 3d space with continuum geometry, which we refer to as a “regular 3d space”.

In fact, Monte Carlo simulation of the Lorentzian type IIB matrix model is not straightforward due to the phase factor $\mbox{e}^{iS_{\rm b}}$ in the partition function, where $S_{\rm b}$ is the bosonic part of the action. The importance sampling is not applicable as it is and one has to face the sign problem if one uses reweighting for this factor. In this work as well as in the previous studies, this problem is avoided by integrating out the scale factor of the bosonic matrices $A_{\mu}$ first and using certain approximation. Here we point out a subtlety in this approximation, and argue that it actually amounts to replacing the phase factor $\mbox{e}^{iS_{\rm b}}$ by a positive weight $\mbox{e}^{\beta S_{\rm b}}$ ($\beta>0$). This new interpretation of the simulations naturally explains not only the emergence of the band-diagonal structure in the spatial matrices $A_{i}$, which is crucial in extracting the real-time evolution, but also the (3+1)d expanding behavior with the Pauli-matrix structure. We also discuss the possibility of obtaining a regular space-time in the original model with the phase factor $\mbox{e}^{iS_{\rm b}}$ without spoiling the (3+1)d expanding behavior. Some results supporting this possibility are reported in a separate paper [31], where the sign problem is overcome by using the complex Langevin method.

The rest of this paper is organized as follows. In section 2 we briefly review the Lorentzian type IIB matrix model. In section 3 we discuss the space-time structure of the matrix configurations obtained by simulation, and show that they are essentially described by the Pauli matrices. In section 4 we provide theoretical understanding of the obtained results, and discuss the possibility of obtaining a regular space-time with the (3+1)d expanding behavior if the sign problem is treated correctly. Section 5 is devoted to a summary and discussions.

2 Brief review of the Lorentzian type IIB matrix model

In this section, we define the Lorentzian type IIB matrix model and its simplified versions, and review some results obtained by Monte Carlo simulations.

2.1 Definition of the Lorentzian type IIB matrix model

The action of the type IIB matrix model is given as [1]

[TABLE]

where $A_{\mu}$ ($\mu=0,1,\cdots,9$) and $\Psi_{\alpha}$ ($\alpha=1,\cdots,16$) are bosonic and fermionic $N\times N$ traceless Hermitian matrices. The indices $\mu$ and $\nu$ are contracted with the Lorentzian metric $\eta_{\mu\nu}=\mathrm{diag}(-1,1,\ldots,1)$. The $16\times 16$ matrices $\Gamma^{\mu}$ and $\mathcal{C}$ are the 10-dimensional gamma matrices and the charge conjugation matrix, respectively, obtained after the Weyl projection. The action (2.1) has a manifest SO(9,1) Lorentz symmetry, under which $A_{\mu}$ and $\Psi_{\alpha}$ transform as a Lorentz vector and a Majorana-Weyl spinor, respectively.

The partition function of the Lorentzian type IIB matrix model is defined as [8]

[TABLE]

where the “$i$” in front of the action is motivated from the fact that the string worldsheet metric has a Lorentzian signature. Note that the bosonic action $S_{\mathrm{b}}$ can be written as

[TABLE]

where we have introduced the Hermitian matrices $F_{\mu\nu}=i\left[A_{\mu},A_{\nu}\right]$. Hence $S_{\mathrm{b}}$ is not positive semi-definite unlike in the Euclidean case. Note also that, unlike in the Euclidean version [32, 33], the matrix integral in (2.4) is divergent because $\mbox{e}^{iS_{\mathrm{b}}}$ is a pure phase factor and the Pfaffian $\mathrm{Pf}\mathcal{M}(A)$ obtained by integrating out the fermionic matrices is a polynomial in $A_{\mu}$.

In order to make the partition function (2.4) finite, we need to introduce the IR cutoffs both in the temporal and spatial directions, for instance, as

[TABLE]

The power $p$ is a parameter, which can be used to test how much the obtained results depend on the way the IR cutoff is introduced [34]. While $p=1$ would be a natural choice, it was proposed that $p$ should be chosen to be a slightly larger value in order to make the results almost independent of $p$. Too large values of $p$ lead to pathological behaviors, however.

The Pfaffian $\mathrm{Pf}\mathcal{M}(A)$ in (2.4) is real in the Lorentzian version unlike in the Euclidean version, where it becomes complex due to the replacement $A_{0}=iA_{10}$. However, the phase factor $\mbox{e}^{iS_{\mathrm{b}}}$ causes the sign problem when one tries to investigate the Lorentzian model by Monte Carlo methods. Here, we avoid this problem666Strictly speaking, the model (2.8) is not completely free of sign-problem because the Pfaffian is real but not positive semi-definite. However, configurations with positive Pfaffian dominates the path integral (2.8) at large $N$, and therefore one can safely replace the Pfaffian by its absolute value in the simulation. following previous work [8, 9, 10] by rewriting the partition function (2.4) as

[TABLE]

where $\theta(x)$ is the Heaviside step function. This can be obtained by integrating out the overall scale factor of the bosonic matrices $A_{\mu}$ first and using certain approximation as discussed in section 4. The parameter $C$ should be set to zero according to the “derivation”, but we generalize the model by choosing $C\neq 0$, which allows us to obtain results for larger matrices in the original $C=0$ model by using smaller matrices [35, 9]. See Appendix B of ref. [9] for the details of the Monte Carlo simulation of the model (2.8).

2.2 SSB of rotational SO(9) symmetry

Next we discuss how one can extract the time-evolution from a given matrix configuration generated by Monte Carlo simulation [8]. Since the eigenvalues of the temporal matrix $A_{0}$ represents time, we work in an SU($N$) basis which diagonalizes $A_{0}$ as

[TABLE]

In this basis, the spatial matrices $A_{i}$ turn out to have an approximate band-diagonal structure. By this, we mean that there exists777In practice, the integer $n$ can be determined by observing the scaling behavior for $\sum_{i}|(A_{i})_{ab}|^{2}$ with $(a+b)/2$ fixed to different values corresponding to different time slices. See section 5 of ref. [10] for the details. some integer $n$ such that the elements of the spatial matrices $\left(A_{i}\right)_{ab}$ for $|a-b|>n$ are much smaller than those for $|a-b|<n$. Thanks to this structure, we can naturally consider the $n\times n$ submatrices $\bar{A}_{i}$

[TABLE]

representing the state at time $t$ defined by

[TABLE]

where $I,J=1,\ldots,n$ and $\nu=0,1,\ldots,N-n$. For example, we can define the extent of the 9d space at time $t$ using $\bar{A}_{i}(t)$ as

[TABLE]

where the symbol “tr” represents a trace over the $n\times n$ submatrix. We can also define the “moment of inertia tensor”

[TABLE]

which is a $9\times 9$ real symmetric tensor. The eigenvalues of $T_{ij}(t)$ represent the spatial extent in each of the nine directions at time $t$, and we denote them by $\lambda_{i}(t)$ with the ordering

[TABLE]

Note that $R^{2}(t)$ and $\lambda_{i}(t)$ are related as

[TABLE]

The expectation values $\langle\lambda_{i}(t)\rangle$ can be used as the order parameters for the spontaneous breaking of the rotational SO(9) symmetry of the model. If the nine eigenvalues do not approach a common value in the large-$N$ limit, we conclude that the SO(9) symmetry is spontaneously broken. From the Monte Carlo simulations of the model (2.8), it was found [8] that the three eigenvalues $\langle\lambda_{i}(t)\rangle$ ($i=1,2,3$) start to grow with $t$ after a critical time $t_{\mathrm{c}}$, which implies that the SO(9) symmetry is spontaneously broken down to SO(3) for $t>t_{\rm c}$. (See refs. [9, 10] for a precise definition of the critical time $t_{\rm c}$, which we use in this work.)

2.3 Expanding behaviors in the simplified models

It is interesting to investigate how the 3d space expands with time. For that, one clearly needs to increase the matrix size, which is very time-consuming due to the existence of the Pfaffian in (2.8). This led to the proposal of the simplified models, the VDM model [9] and the bosonic model [10], which amounts to replacing the Pfaffian as

[TABLE]

where $\Delta(\alpha)\equiv\prod_{a>b}^{N}(\alpha_{a}-\alpha_{b})$ is the van der Monde (VDM) determinant. This replacement reduces the computational cost from $O(N^{5})$ to $O(N^{3})$, which enables simulations with considerably large matrix size. These two models are expected to describe the qualitative behaviors of the original model at early times and at late times, respectively.

In both these models, the spontaneous breaking of the SO(9) rotational symmetry to SO(3) was observed after some critical time as in the original model, and the rate of expansion at late times was investigated. In the VDM model, the extent of space $R(t)$ defined in (2.12) exhibits an exponential growth [9]

[TABLE]

which is reminiscent of inflation888This behavior was observed also in the original model [36] although the matrix size used was not large enough to confirm the long-time behavior. , and this behavior does not seem to change with increasing $t$. In the bosonic model, on the other hand, the exponential expansion observed at early times changes into a power-law expansion [10]

[TABLE]

at later times, which is reminiscent of the Friedmann-Robertson-Walker Universe at the radiation dominated era. Based on these results, it has been speculated that the extent of space $R(t)$ in the original model shows an exponential growth at early times and a power-law expansion at later times. If true, it implies that the e-folding or the duration of the cosmic inflation may be determined dynamically in the original model.

3 Space-time structure of the matrix configurations

In this section, we investigate the space-time structure of the matrix configurations generated by the Monte Carlo simulation of the model (2.8) and the simplified models (2.16).

3.1 Results for the bosonic model

In this subsection, we consider the bosonic model, which is a simplified model for the late time behaviors. Let us first look at the basic quantities such as the extent of space $R^{2}(t)$ and the eigenvalues $\langle\lambda_{i}(t)\rangle$ of $T_{ij}(t)$. In Fig. 1 we plot the extent of space $R^{2}(t)/R^{2}(t_{\mathrm{c}})$ (Top-Left) and the normalized eigenvalues $\langle\lambda_{i}(t)\rangle/R^{2}(t_{\mathrm{c}})$ of $T_{ij}(t)$ (Top-Right) against $(t-t_{\mathrm{c}})/R(t_{\mathrm{c}})$ for $N=256$, $C=100$, $\kappa=1.0$ with the block size $n=18$ in (2.12). Here and for all the other plots in Fig. 1, we only present the results in the $t<0$ region since the results are symmetric999This does not mean that the Big Crunch occurs in this model because the time difference between the symmetric point $t=0$ and the critical time $t=t_{\rm c}$ seems to diverge in physical units in an appropriate large-$N$ limit. See section 3.3. under the time reflection $t\mapsto-t$. The power $p$ in the IR cutoff (2.6) and (2.7) is chosen to be $p=1.5$, which is found to be large enough to make the results almost independent of $p$ (See Appendix A.). Let us recall that $R^{2}(t)$ is related to $\langle\lambda_{i}(t)\rangle$ through (2.15). While the extent of space $R^{2}(t)/R^{2}(t_{\mathrm{c}})$ grows with $t$ for $t>t_{\rm c}$, it is only three out of nine eigenvalues of $T_{ij}(t)$ that grow with $t$, which suggests that the rotational SO(9) symmetry is broken spontaneously to SO(3). These results are analogous to the previous results obtained for $p=1$ [10].

The simplest way to probe the space-time structure is to define an $n\times n$ matrix

[TABLE]

which is invariant under SO(9) rotations. Let us denote its eigenvalues as $q_{k}(t)$ ($k=1,\cdots,n$) with the ordering

[TABLE]

These eigenvalues tell us how the space spreads in the radial direction at each time $t$.

In Fig. 1 (Middle-Left), we plot the eigenvalues $q_{k}(t)/R^{2}(t_{\mathrm{c}})$ against $(t-t_{\rm c})/R(t_{\rm c})$. We find that the two largest eigenvalues grow with $t$, but not the others. Let us note that the eigenvalues of $Q(t)$ are related to the extent of space $R^{2}(t)$ as

[TABLE]

This implies that the time-dependence of $R^{2}(t)$ seen in the Top-Left panel is caused only by the two largest eigenvalues of $Q(t)$.

Let us next discuss the space-time structure in the three extended directions and the six shrunken directions separately. Since we are dealing with spontaneous symmetry breaking, we need to choose the frame properly in order to distinguish these directions. Suppose $v_{j}^{(i)}(t)$ ($j=1,\cdots,9$) are the normalized eigenvectors of the “moment of inertia tensor” (2.13) corresponding to the eigenvalues $\lambda_{i}(t)$ with the ordering (2.14). Then, we can define the $n\times n$ matrix corresponding to the spatial direction with the extent $\lambda_{i}$ as

[TABLE]

and its eigenvalues $a^{(i)}_{k}(t)$ ($k=1,\cdots,n$) with the ordering

[TABLE]

In Fig. 1 (Middle-Right), we plot the eigenvalues $a^{(1)}_{k}(t)/R(t_{\mathrm{c}})$ against $(t-t_{\mathrm{c}})/R(t_{\mathrm{c}})$. We find that only two eigenvalues $a_{1}^{(1)}(t)$ and $a_{n}^{(1)}(t)$ grow in magnitude with time $t$, and all the others remain close to zero. Similar behaviors are seen also for the eigenvalues $a^{(2)}_{k}(t)$ and $a^{(3)}_{k}(t)$ obtained for the other extended directions. In Fig. 1 (Bottom-Left), we zoom up the same plot to make visible the eigenvalues close to zero. In Fig. 1 (Bottom-Right), we plot the eigenvalues $a^{(4)}_{k}(t)/R(t_{\mathrm{c}})$ against $(t-t_{\mathrm{c}})/R(t_{\mathrm{c}})$. We find that all the eigenvalues remain close to zero. Similar behaviors are seen also for the eigenvalues $a^{(5)}_{k}(t),\cdots,a^{(9)}_{k}(t)$ obtained for the other shrunken directions. Comparing the two plots at the bottom of Fig. 1, we notice that the eigenvalue distribution of $\bar{A}^{(i)}$ is almost identical for the extended directions and the shrunken directions except for the two eigenvalues with large magnitude.

Similarly to (3.3), the eigenvalues of $\bar{A}^{(i)}(t)$ are related to the extent of space $\lambda_{i}(t)$ in the $i$th direction as

[TABLE]

Our observation implies that the spontaneous symmetry breaking of the SO(9) rotational symmetry seen in the Top-Right panel is caused only by the two eigenvalues of $\bar{A}^{(i)}(t)$ with large magnitude.

3.2 Including fermionic contributions

In order to seek for the possibility to obtain a regular space-time, we repeat the analysis in the previous subsection in the case of the original model (2.8) including fermionic contributions. Since the cost of Monte Carlo simulations increases from O($N^{3}$) to O($N^{5}$), here we restrict ourselves to a rather small matrix size $N=16$.

In Fig. 2 we plot the same quantities as in Fig. 1 for the original model with $N=16$, $C=3.91$, $\kappa=0.38$ and the block size $n=6$. The power $p$ in the IR cutoff (2.6) and (2.7) is chosen to be $p=1.6$, which is found to be large enough to make the results almost independent of $p$ (See Appendix A.). These results are qualitatively the same as those obtained for the bosonic model. While the fermionic matrices are expected to play an important role in the properties of the model such as the expanding behavior, they do not seem to affect the singular space-time structure.

3.3 Taking the continuum limit

As yet another possibility to obtain a regular space-time, let us consider taking the continuum limit. Here we use the VDM model, which is a simplified model for the early time behaviors. In Fig. 3 (Top-Left), we plot the extent of space $R^{2}(t)/R^{2}(t_{\mathrm{c}})$ against time $(t-t_{\mathrm{c}})/R(t_{\mathrm{c}})$ for various $N$, $C$ and $\kappa$ with the block size $n$ listed in table 1. The power $p$ in the IR cutoff (2.6) and (2.7) is chosen as $p=1.4$ following ref. [34]. From this plot, we observe a clear scaling behavior for $(t-t_{\rm c})/R(t_{\rm c})\lesssim 0.40$.

In Fig. 3 (Top-Right), we plot the normalized eigenvalues $\langle\lambda_{i}(t)\rangle/R^{2}(t_{\mathrm{c}})$ of $T_{ij}(t)$ for the VDM model with $N=96$, $C=0$ and $\kappa=2$. Similar behaviors are obtained for the other parameter sets. We find that three out of nine eigenvalues of $T_{ij}(t)$ grow with time, which suggests that the rotational SO(9) symmetry is broken spontaneously to SO(3) for $t>t_{\rm c}$. These results are similar to those obtained in refs. [9, 34].

In order to discuss the continuum limit, let us define the “volume” and the “lattice spacing” in the temporal direction as [9]

[TABLE]

where $t_{\mathrm{peak}}$ represents the position of the peak in $R^{2}(t)$ and $\nu$ is the number of data points of $R^{2}(t)$ contained within $t_{\rm c}<t\leq t_{\mathrm{peak}}$. Roughly speaking, the lattice spacing $\epsilon$ represents the average horizontal spacing between the adjacent data points of $R^{2}(t)/R^{2}(t_{\mathrm{c}})$. In table 1, we present the volume $\Delta$ and the lattice spacing $\epsilon$ obtained for each parameter set $(N,C,\kappa)$ used in Fig. 3. The deviation from the scaling behavior for $(t-t_{\rm c})/R(t_{\rm c})>0.40$ seen in Fig. 3 can be understood either as the finite volume effects or as the finite lattice spacing effects depending on the parameter set.

In what follows, we focus on the point $(t-t_{\rm c})/R(t_{\rm c})\sim 0.40$, at which the results for $R^{2}(t)/R^{2}(t_{\rm c})$ with the four parameter sets agree with each other. In Fig. 3 (Middle-Left), we plot the normalized eigenvalues $\langle q_{k}(t)\rangle/R^{2}(t_{\rm c})$ $(k=1,\cdots,n)$ of $Q(t)$ against their label $(k-1)/(n-1)$ for the four parameter sets. This reveals a clear scaling behavior except for the two largest eigenvalues, which grow as the lattice spacing $\epsilon$ decreases. This scaling behavior is consistent with the scaling of the ratio $R^{2}(t)/R^{2}(t_{\rm c})$ in the continuum limit [9, 10] seen in the Top-Left panel considering the relation (3.3). Note, however, that the time dependence of $R^{2}(t)/R^{2}(t_{\rm c})$ is caused by the two largest eigenvalues of $Q(t)$ as we have seen in the previous subsections. Therefore, the scaling of $R^{2}(t)/R^{2}(t_{\rm c})$ implies that the two largest eigenvalues of $Q(t)$ should grow linearly in $n$ in the continuum limit. This is confirmed numerically in Fig. 4 (Left) assuming the presence of $1/n$ corrections.

Let us next consider the space-time structure in the extended directions and the shrunken directions separately. In Fig. 3 (Middle-Right), we plot the eigenvalues of $\bar{A}^{(1)}(t)/R(t_{\rm c})$ obtained at $(t-t_{\rm c})/R(t_{\rm c})\approx 0.40$ against the label $(k-1)/(n-1)$. Here again we observe a clear scaling behavior except for the ones at both ends of the spectrum. Similar behaviors are obtained for the other extended directions. According to the same argument as in the previous paragraph, we can deduce that the normalized eigenvalues at both ends of the spectrum grow in magnitude as O($\sqrt{n}$) in the continuum limit, which is confirmed in Fig. 4 (Right) assuming the presence of $1/n$ corrections.

In Fig. 3 (Bottom), we plot the eigenvalues of $\bar{A}^{(4)}(t)/R(t_{\rm c})$ obtained at $(t-t_{\rm c})/R(t_{\rm c})\approx 0.40$ against the label $(k-1)/(n-1)$. We observe a clear scaling behavior here as well. In fact, the eigenvalues are almost the same as those for the extended directions except for the ones at both ends. Similar behaviors are obtained for the other shrunken directions.

Thus we find in the VDM model that the singular space-time structure becomes even more pronounced in the continuum limit instead of getting milder. It is surprising that the two eigenvalues of $\bar{A}^{(i)}(t)/R(t_{\mathrm{c}})$ ($i=1,2,3$ ) actually diverges in the continuum limit although the extent of space defined by $R^{2}(t)/R^{2}(t_{\rm c})$ remains finite. It is these two eigenvalues that cause the spontaneous breaking of the SO(9) rotational symmetry and the expansion of space. All the other eigenvalues of $\bar{A}^{(i)}(t)/R(t_{\mathrm{c}})$ remain finite and contribute only to the time-independent SO(9) symmetric part of the “moment of inertia tensor” $T_{ij}(t)$.

3.4 The Pauli-matrix structure

In this subsection, we provide deeper understanding of the singular space-time structure observed in the previous subsections. Let us work in the SU($n$) basis which diagonalizes $Q(t)$ at each time $t$ with the ordering (3.2), and consider the $2\times 2$ submatrix $X_{i}(t)$ in the bottom-right corner of

[TABLE]

for the extended directions $i=1,2,3$. Here we use the VDM model with the parameter sets given in table 1 and take the continuum limit focusing on the time $(t-t_{\rm c})/R(t_{\rm c})\approx 0.40$ as we did in section 3.3.

We show below that the three matrices $X_{i}$ in (3.8) tend to satisfy the SU(2) Lie algebra

[TABLE]

for some real constant $c$ in the continuum limit. In order to determine the optimal value of $c$, we consider a quantity

[TABLE]

which represents the violation of the relation (3.9). The value of $c$ that minimizes $S(c)$ can be readily obtained as

[TABLE]

Using $c=\tilde{c}$ as the optimal value for each configuration, we investigate to what extent the relation (3.9) is satisfied.

In Fig. 5, we show a scatter plot for the real part (Left) and the imaginary part (Right) of each side of (3.9). The quantities on both sides are normalized by ${\rm tr}(X_{l}^{2})$ so that they become invariant under the scale transformation $X_{i}\mapsto{\rm const.}X_{i}$. We observe that the data points tend to converge to the line $y=x$ as one goes from the top to the bottom corresponding to decreasing the lattice spacing $\epsilon$ (See table 1.). This shows that the $2\times 2$ matrices $X_{i}$ ($i=1,2,3$) tend to satisfy (3.9) in the continuum limit.

Thus we conclude that the singular space-time structure observed for the matrix configurations generated by simulations is essentially associated with the Pauli matrices. The Pauli matrices may be regarded as the simplest matrix configuration that has SO(3) symmetry in the sense that their SO(3) rotation can be absorbed by an appropriate SU($N$) transformation. Given the situation characterized by the two large eigenvalues of $Q(t)$, the appearance of the Pauli-matrix structure may not be that surprising.

4 The new interpretation of the simulation

In this section, we attribute the observed Pauli-matrix structure to the approximation involved in deriving the partition function (2.8), which was used in Monte Carlo simulation. We point out a subtlety in the approximation, and argue that the approximation amounts to replacing $\mbox{e}^{iS_{\mathrm{b}}}$ by $\mbox{e}^{\beta S_{\mathrm{b}}}$ in the original partition function (2.4). This new interpretation of the simulation provides us with a natural understanding of the (3+1)d expanding behavior with the Pauli-matrix structure discussed in section 3. We also speculate on a possible scenario for the original model with the correct $\mbox{e}^{iS_{\mathrm{b}}}$ factor.

4.1 The “derivation” of

the partition function (2.8)

Let us first review how one can obtain the partition function (2.8) used in Monte Carlo simulation from the original partition function (2.4). (This was done in Appendix A of ref. [9] for $p=1$, but here we generalize it to arbitrary $p$.)

Note that the integrand of the partition function (2.4) involves a phase factor $\mbox{e}^{iS_{\rm b}}$. As is commonly done in integrating oscillating functions, we introduce the convergence factor $\mbox{e}^{-\epsilon|S_{\rm b}|}$ and take the $\epsilon\rightarrow 0$ limit after the integration.

The partition function can then be rewritten as

[TABLE]

where $\kappa$ and $L$ are the cutoff parameters introduced in (2.6) and (2.7), respectively. Rescaling the variables $A_{\mu}\mapsto r^{1/2}A_{\mu}$ in the integrand, we get

[TABLE]

Here we have defined the function $f(S_{\rm b})$ by

[TABLE]

which is a complex-valued function with the property $f(-S_{\rm b})=f(S_{\rm b})^{*}$.

For $|S_{\rm b}|\ll\frac{1}{L^{4}}$, the function can be well approximated by

[TABLE]

For $|S_{\rm b}|\gtrsim\frac{1}{L^{4}}$, on the other hand, the phase of the integrand in (4.3) starts to oscillate violently in the region $r\gtrsim 1/\sqrt{|S_{\rm b}|}$, and hence the integral decreases rapidly in magnitude for increasing $|S_{\rm b}|$. In particular, the asymptotic behavior of $f(S_{\rm b})$ for $S_{\rm b}\gg\frac{1}{L^{4}}$ can be estimated as

[TABLE]

by deforming the integration contour in (4.3).

Recalling eq. (2.5), the condition $|S_{\rm b}|\ll\frac{1}{L^{4}}$ for (4.4) can be rewritten as

[TABLE]

Therefore, assuming that the right-hand side $\frac{4}{NL^{4}}$ of (4.6) becomes small at large $N$, we may make a replacement

[TABLE]

up to a normalization constant. For the bosonic model and the VDM model, one simply has to replace the Pfaffian in (4.1) and (4.2) as (2.16).

4.2 Subtlety in the derivation and the new interpretation

The only step in the derivation that may go wrong is the replacement (4.7). The subtlety in this replacement can be seen as follows. Note that the phase factor $\mbox{e}^{iS_{\rm b}}$ in the partition function (2.4) favors configurations at which the bosonic action $S_{\rm b}$ is stationary. On the other hand, the above approximation essentially replaces the phase factor $\mbox{e}^{iS_{\rm b}}$ by the delta function $\delta(S_{\rm b})$, which amounts to picking up configurations at which $S_{\rm b}$ is stationary only under rescaling $A_{\mu}\mapsto{\rm const.}A_{\mu}$. While it is true that $|f(S_{\rm b})|$ is sharply peaked at $S_{\rm b}=0$, the function $f(S_{\rm b})$ is actually a complex-valued function, whose phase rotates violently around $S_{\rm b}=0$. This effect of the phase should be responsible for favoring the configurations at which $S_{\rm b}$ is stationary. The approximation ignores this effect completely, and hence it cannot be justified.

If the model (2.8) is not equivalent to the original model (2.4), what kind of model does it actually correspond to? Here we point out that the constraint on $S_{\rm b}$ that appears in (2.8) may be regarded as the constraint one uses in defining a microcanonical ensemble. From this viewpoint, we consider that the model (2.8) is actually equivalent to the corresponding canonical ensemble with the Boltzmann weight $\mbox{e}^{\beta S_{\rm b}}$. The real parameter $\beta$ depends on the parameter $C$ in the constraint101010This connection also provides clear justification of the renormalization-group-like method [9, 35], which amounts to tuning the parameter $C$ in order to obtain the late-time behaviors with smaller matrix size.. As we will see below, we consider that the model (2.8) corresponds essentially to replacing $\mbox{e}^{iS_{\rm b}}$ by $\mbox{e}^{\beta S_{\rm b}}$ with $\beta>0$.

For $\beta>0$, the first term in (2.5) that appears in $\mbox{e}^{\beta S_{\rm b}}$ favors configurations in which $A_{0}$ and $A_{i}$ commute. This means that the spatial matrices $A_{i}$ tend to become diagonal in the SU($N$) basis which diagonalizes $A_{0}$. On the other hand, the second term in (2.5) favors configurations in which the noncommutativity among the spatial matrices $A_{i}$ is large. The band-diagonal structure, which plays a crucial role in extracting the real-time evolution as in section 2.2, can be understood as a consequence of the balance of these two effects.

We can also understand the reason for the (3+1)d expanding behavior with the Pauli-matrix structure. Here we assume that the first term in (2.5) is not important except in realizing the band-diagonal structure and focus on the effect of the second term in (2.5), which favors large $\mathrm{Tr}\left(F_{ij}\right)^{2}$, where $F_{ij}=i\left[A_{i},A_{j}\right]$. We also have to take into account the constraint $\frac{1}{N}\mathrm{Tr}\left\{\left(A_{i}\right)^{2}\right\}^{p}=1$, where we set $p=1$ in what follows.

Simplifying the band-diagonal structure of the spatial matrices $A_{i}$ ($i=1,\cdots,9$), we consider the block-diagonal structure given as

[TABLE]

where $n$ is the common block size and $B$ is the number of blocks satisfying $N=nB$. Within this ansatz, we would like to maximize $\mbox{Tr\,}(F_{ij})^{2}$ under the constraint $\frac{1}{N}\mbox{Tr\,}(A_{i})^{2}=1$. Note that we have

[TABLE]

where we have defined $\bar{F}^{(b)}_{ij}=i[\bar{A}^{(b)}_{i},\bar{A}^{(b)}_{j}]$ for each block $b$.

Let us solve the maximization problem in two steps. First we fix

[TABLE]

and maximize $\mbox{Tr\,}(F_{ij})^{2}$ under this constraint. Following the discussion given in ref. [8], the solution to this first maximization problem can be written in terms of the Pauli matrices $\sigma_{i}$ as

[TABLE]

for $i=1,2,3$ and $\bar{A}^{(b)}_{i}=0$ otherwise, up to the symmetries of the problem such as the SO(9) rotational symmetry and the SU($n$) symmetry within each block. The value of $\mbox{Tr\,}(F_{ij})^{2}$ for (4.13) is given as

[TABLE]

As the second step of the maximization, we maximize (4.14) under the constraint (4.12). The maximum is given when all but one of the $r_{b}$’s are zero.

In reality, one should also take into account the entropic factor due to quantum fluctuations, which is expected to favor certain distribution of $r_{b}$. Due to the time-reversal symmetry $A_{0}\mapsto-A_{0}$ of the model, the most natural distribution would be that $r_{b}$ is large around $t=0$ and decreases with $|t|$. Thus we can understand the appearance of the (3+1)d expanding behavior with the Pauli-matrix structure.

4.3 A possible scenario for the original model

In the previous subsections, we have argued that the model (2.8) used for Monte Carlo simulation actually corresponds to a model with $\mbox{e}^{\beta S_{\rm b}}$ instead of $\mbox{e}^{iS_{\rm b}}$ in (2.4). This new interpretation explains naturally the (3+1)d expanding behavior with the Pauli-matrix structure. The crucial question then is what happens for the model with the correct $\mbox{e}^{iS_{\rm b}}$ factor. It is not easy to answer this question due to the sign problem, which occurs because $\mbox{e}^{iS_{\rm b}}$ is a pure phase factor and one cannot regard the integrand of the partition function (2.4) as the probability distribution. Here we speculate on a possible scenario based on the results obtained so far.

For that purpose, let us consider a generalized model with a factor $\mbox{e}^{\beta(\cos\theta+i\sin\theta)S_{\rm b}}$ ($0\leq\theta\leq\pi/2$), which interpolates the two models. At $\theta=0$, we obtain the model with the positive definite factor $\mbox{e}^{\beta S_{\rm b}}$ we have been studying, whereas at $\theta=\pi/2$, we obtain the model with $\mbox{e}^{i\beta S_{\rm b}}$ we are aiming at. The scale parameter $\beta$ can be absorbed, if one wishes, by the redefinition $A_{\mu}\mapsto\beta^{-1/4}A_{\mu}$ and the replacement $L\mapsto\beta^{1/4}L$ in (2.7).

As far as $\theta<\pi/2$, the real part of the coefficient of $S_{\rm b}$ is positive. Therefore, certain effects favoring the band-diagonal structure and the Pauli-matrix structure in $A_{i}$ are at work. Note also that the classical equation of motion is common to all values of $\theta$. In fact, the classical equation of motion becomes valid at late times if the expansion of space occurs because each term in the bosonic action becomes large [14, 15]. Therefore, if some classical solution dominates for $\theta=0$, the same solution may well dominate also for other $\theta$ less than some value $\theta_{0}$. From this argument, we speculate that the models with $0\leq\theta\leq\theta_{0}$ are qualitatively the same.

As one approaches $\theta=\pi/2$, the real part of the coefficient of $S_{\rm b}$ becomes small, and different classical solutions may dominate. Note that the matrix configurations with the Pauli-matrix structure are obtained essentially by maximizing $S_{\rm b}$, but the classical solutions that can be obtained by extremizing $S_{\rm b}$ instead of maximizing it should have more variety. Indeed we have generated numerically many classical solutions that have (3+1)d expanding behavior and find for all of them that the matrix $Q(t)$ defined in (3.1) has a smooth eigenvalue distribution [37]. This is understandable since the configurations with the Pauli-matrix structure are actually disfavored entropically. Recall, for instance, that only two eigenvalues of the matrix $Q(t)$ are large, meaning that the entropy for such configurations must be small. It should be mentioned, however, that from the above classical analysis alone, one cannot single out the (3+1)d expanding space-time because there are also other solutions with different dimensionality. Whether the (3+1)d expanding behavior remains even for $\theta\sim\pi/2$ is therefore a highly nontrivial question.

5 Summary and discussions

In this paper we have investigated the space-time structure of the matrix configurations obtained in Monte Carlo studies of the Lorentzian type IIB matrix model and the simplified models. In these models, the time-evolution can be extracted from the matrix configurations by working in the SU($N$) basis which diagonalizes the temporal matrix $A_{0}$. The $n\times n$ spatial submatrices $\bar{A}_{i}(t)$ ($i=1,\cdots,9$) at each time $t$ show that only three out of nine directions expand after some critical time suggesting the SSB of rotational SO(9) symmetry to SO(3). By calculating the eigenvalues of $\bar{A}_{i}(t)$ at each $t$, however, we have found that only two of them increase in magnitude with $t$ in the extended directions, while the rest are independent of $t$ and SO(9) symmetric. This implies that the SSB is caused only by the two eigenvalues. In the continuum limit, the magnitude of the two eigenvalues diverges in physical units and the spatial matrices $\bar{A}_{i}(t)$ approach a configuration which is essentially described by the Pauli matrices.

We have attributed this problem to the approximation used in Monte Carlo simulation to avoid the sign problem, which actually amounts to replacing $\mbox{e}^{iS_{\mathrm{b}}}$ by $\mbox{e}^{\beta S_{\mathrm{b}}}$ in the partition function (2.4) of the Lorentzian type IIB matrix model. This new interpretation of the Monte Carlo simulation enables us to understand the interesting aspects of the obtained results such as the band-diagonal structure of the spatial matrices $A_{i}$ as well as the appearance of the (3+1)d expanding behavior with the Pauli-matrix structure.

In order to discuss what happens in the original model, we have considered a model with a factor $\mbox{e}^{\beta(\cos\theta+i\sin\theta)S_{\rm b}}$, which interpolates the model we have been studying ($\theta=0$) and the model we are aiming at ($\theta=\pi/2$). Using some arguments based on the classical equation of motion, which is common to all $\theta$, we have speculated that it is possible to obtain a regular space-time structure with the (3+1)d expanding behavior by approaching $\theta=\pi/2$ in the large-$N$ limit. The crucial point is that the Pauli-matrix structure is obtained by maximizing the action at the expense of reducing the entropy. By approaching $\theta=\pi/2$, one may obtain classical solutions which only extremize the action that have larger entropy due to a smooth eigenvalue distribution of the matrix $Q(t)$. The existence of such classical solutions with the (3+1)d expanding behavior has been confirmed numerically [37]. Whether such classical solutions appear from the full quantum theory by approaching $\theta=\pi/2$ remains to be seen.

Monte Carlo simulation of the interpolating model for $\theta\neq 0$ is difficult since the complex weight $\mbox{e}^{\beta(\cos\theta+i\sin\theta)S_{\rm b}}$ causes the sign problem. As a promising approach to overcome this problem, we may use the complex Langevin method [38, 39], which has attracted much attention recently [40, 41, 42, 43, 44, 45, 46]. It was successful also in investigating the SSB of rotational symmetry in the 6d Euclidean type IIB matrix model [7]. Preliminary results [31] for the bosonic Lorentzian model suggest that by approaching $\theta=\pi/2$, one obtains clear deviations from the Pauli-matrix structure without losing the (3+1)d expanding behavior. We hope to see whether a regular (3+1)d expanding space-time emerges or not in the near future.

Acknowledgements

The authors would like to thank K.N. Anagnostopoulos, T. Azuma, K. Hatakeyama, S. Iso, H. Kawai, A. Matsumoto, H. Steinacker and A. Yosprakob for valuable discussions. This research used computational resources of the K computer of the HPCI system provided by the AICS through the HPCI System Research Project “Quest for the ultimate laws and the birth of the universe” (Project ID : hp170229, hp180178). Computation was carried out also on PC clusters at KEK. The supercomputer FX10 at University of Tokyo has been used in developing our code for parallel computing. A. T. was supported in part by Grant-in-Aid for Scientific Research (No. 18K03614) from Japan Society for the Promotion of Science.

Appendix A The determination of the parameter $p$

In this appendix, we explain how we determine the parameter $p$ in the IR cutoff (2.6) and (2.7). While a naive choice would be $p=1$, it was proposed in ref. [34] that one should choose a slightly larger value so that the results become almost independent of $p$. There it was found in the VDM model that the results for the extent of space $R^{2}(t)$ become independent of $p$ when $p$ is larger111111For the values of $p$ in this region, it was also observed [34] from the analysis of the Schwinger-Dyson equations that the effect of the IR cutoff decreases as one takes the infinite volume limit. than $p_{\rm c}=1.2\sim 1.3$. Based on this observation, we used $p=1.4$ when we simulate the VDM model in section 3.3.

Here we repeat the same analysis in the case of the bosonic model and the original model. In Fig. 6, we plot the extent of space $R^{2}(t)/R^{2}(t_{\rm c})$ against time $(t-t_{\rm c})/R(t_{\rm c})$ for the bosonic model (Left) and the original model (Right), respectively, with various values of $p$. For all values of $p$, we find that only three directions start to expand at some critical time $t_{\rm c}$. In the bosonic model, the results scale for $p=1.3,1.4,1.5$ except for the data around the peak of $R^{2}(t)$. Similar scaling behavior is observed for the original model for $p=1.4,1.5,1.6$. Based on these results, we use $p=1.5$ for the bosonic model and $p=1.6$ for the original model in sections 3.1 and 3.2, respectively.