Power-law expansion of the Universe from the bosonic Lorentzian type IIB matrix model

Recent studies on the Lorentzian version of the type IIB matrix model show that (3+1)D expanding universe emerges dynamically from (9+1)D space-time predicted by superstring theory. Here we study a bosonic matrix model obtained by omitting the fermionic matrices. With the adopted simplification and the usage of a large-scale parallel computer, we are able to perform Monte Carlo calculations with matrix size up to N = 512, which is twenty times larger than that used previously for the studies of the original model. When the matrix size is larger than some critical value Nc ≃ 110, we find that (3+1)D expanding universe emerges dynamically with a clear large-N scaling property. Furthermore, the observed increase of the spatial extent with time t at sufficiently late times is consistent with a power-law behavior t1/2, which is reminiscent of the expanding behavior of the Friedmann-Robertson-Walker universe in the radiation dominated era. We discuss possible implications of this result on the original supersymmetric model including fermionic matrices.


Introduction
Understanding how our universe began is a fascinating subject in theoretical physics. This is, of course, not so easy since one has to deal with quantum gravity near cosmic singularity, which appears in general relativity. As a promising way to describe quantum gravity, string theory has been studied for many years. However, as far as one applies perturbative string theory around a time dependent background, cosmic singularity cannot be resolved in general [1][2][3][4]. It is therefore necessary to study string theory in a nonperturbative background-independent fashion. As nonperturbative formulations of superstring/M theory, matrix models that can be obtained formally by dimensionally reducing 10D N = 1 super Yang-Mills theory to d = 0 [5], d = 1 [6] and d = 2 [7] were proposed. As a closely related direction, ref. [8,9] proposes a conformal field theory, which is holographically dual to inflationary models.
The type IIB matrix model [5] is one of these proposals corresponding to the d = 0 case above. A peculiar feature of this model is that both space and time emerge dynamically from the matrix degrees of freedom. In this context, the idea of emergent gravity has been pursued [10][11][12][13][14][15][16][17] in gauge theories on non-commutative space that appear from the type IIB matrix model for a particular class of backgrounds [18][19][20][21]. There are also other proposals for the description of curved space-time in matrix models [22][23][24].
Until recently, the type IIB matrix model was studied after the "Wick rotation" [25][26][27][28][29][30][31][32][33][34][35][36][37][38] since the partition function of the Euclidean model obtained in this way is shown to be finite [39,40]. However, the Euclidean model is clearly not suitable for studying the realtime dynamics since the time coordinate is treated as purely imaginary due to the "Wick JHEP11(2015)070 rotation". Moreover, it is known that the Wick rotation is more subtle in quantum gravity than in quantum field theories at the nonperturbative level (See, for instance, refs. [41,42]). In fact, according to a recent study of the Euclidean model using the Gaussian expansion method, the emergent space-time seems to have only three dimensions [43].
On the other hand, the Lorentzian version of the type IIB matrix model has been studied for the first time in ref. [44]. Unlike the Euclidean version, one has to introduce infrared cutoffs in the temporal and spatial directions to make the partition function finite. Despite this subtlety, the time-evolution of the "universe" was extracted from the matrix configurations generated by Monte Carlo simulation, and a large-N scaling behavior was observed with the matrix size N ≤ 16. Quite remarkably, the SO(9) rotational symmetry of the 9D space turned out to be spontaneously broken down to SO(3) at some critical time, after which only three out of the nine spatial directions start to expand.
In order to see what happens at later times in this model, one needs to increase the matrix size, which makes the Monte Carlo simulation more and more time-consuming. As an alternative approach, one can use the classical approximation [45][46][47][48][49] to investigate possible behaviors at late times since each term in the action becomes large as the expansion of the "universe" proceeds. A general prescription to construct solutions to the classical equations of motion was given in ref. [46]. One can actually construct classical solutions corresponding to an expanding (3+1)D universe, which naturally solve the cosmological constant problem [46]. As a closely related progress, it was found that matrix configurations with intersecting fuzzy spheres in the extra dimensions can accommodate the standard model fermions [50][51][52][53][54][55][56].
In fact, it is known that the classical equations of motion of the matrix model have infinitely many solutions [46]. Therefore, in order to determine which classical solution is actually realized at late times, we need to study the time-evolution of the "universe" at least for a sufficiently long time by performing Monte Carlo simulation. To that end, we previously studied a simplified model that describes the early time behaviors of the original model [57]. With the matrix size N ≤ 64, we observed a clear exponentially expanding behavior, which is reminiscent of the inflation. Monte Carlo studies of the original model with N = 24 [58] yielded results consistent with this observation.
In this paper we study a bosonic model, which can be obtained by simply omitting the fermionic matrices. This simplification and the usage of a large-scale parallel computer enable us to perform Monte Carlo simulation with the matrix size up to N = 512. Unlike the original model, the eigenvalue distribution of the temporal matrix turns out to have a finite extent without introducing a cutoff in the temporal direction. The extent of the eigenvalue distribution is independent of N for N < N c 110, but it increases with N for N ≥ N c . We find that the properties of the model changes drastically at the critical N = N c . For N < N c , the dominant matrix configurations do not allow extraction of a well-defined time-evolution. For N ≥ N c , on the other hand, we can extract a meaningful time-evolution, which shows that the SO(9) rotational symmetry is broken spontaneously down to SO(3) symmetry at some point in time similarly to the original model. The large-N scaling behavior is clearly observed. The expanding behavior of the spatial extent can be fitted by an exponential function only for a rather short period, and JHEP11(2015)070 after that the expansion seems to slow down. In fact, we find that the expanding behavior at late times is consistent with a power-law t 1/2 , which is the expanding behavior of the Friedmann-Robertson-Walker universe in the radiation dominated era. From these results, we speculate that the exponential expansion of the space in the original model suggested in the previous works actually terminates at some point in time and turns into a power law similarly to the bosonic model. This would imply that the number of e-foldings is determined dynamically in the Lorentzian type IIB matrix model.
The rest of this paper is organized as follows. In section 2 we briefly review some important properties of the Lorentzian type IIB matrix model. In section 3 we define the bosonic model we study in this paper, and discuss the existence of the critical N , at which the properties of the model change drastically. In sections 4 and 5 we discuss in detail the properties of the model below and above N c , respectively. Section 6 is devoted to a summary and discussions. In appendix A we give the details of our simulation. In appendix B we present our results for the (5+1)D version of the bosonic model.

Brief review of the Lorentzian type IIB matrix model
The action of the Lorentzian type IIB matrix model is given by [5] 2) where the bosonic N × N matrices A µ (µ = 0, . . . , 9) and the fermionic matrices Ψ α (α = 1, . . . , 16) are both traceless and Hermitian. Γ µ are 10D gamma-matrices after the Weyl projection and C is the charge conjugation matrix. The "coupling constant" g is merely a scale parameter in this model since it can be absorbed by rescaling A µ and Ψ appropriately. The indices µ and ν are contracted using the Lorentzian metric η µν = diag (−1, 1, . . . , 1). The Euclidean version can be obtained by making the "Wick rotation" A 0 = iA 10 , where A 10 is supposed to be Hermitian. The partition function for the Lorentzian version is proposed in ref. [44] as with the action (2.1). The "i" in front of the action is motivated from the fact that the string world-sheet metric should also have a Lorentzian signature. By integrating out the fermionic matrices, we obtain the Pfaffian which is real unlike in the Euclidean case [38]. Note also that the bosonic action (2.2) can be written as

JHEP11(2015)070
where we have introduced the Hermitian matrices F µν = i [A µ , A ν ]. Since the two terms in the last expression have opposite signs, S b is not positive semi-definite, and it is not bounded from below.
In order to make the partition function (2.4) finite, one needs to introduce infrared cutoffs in both the temporal and spatial directions, for instance, as [44] It was found that the two parameters κ and Λ can be removed in the large-N limit in such a way that physical quantities scale [44]. Therefore the resulting theory can be characterized by a single scale parameter.
In the present work, it will be important to understand the reason why we need to introduce the cutoff (2.7) in the temporal direction. Note first that one can use the SU (N ) symmetry of the model to bring the temporal matrix A 0 into the diagonal form where α 1 < · · · < α N . (2.9) By "fixing the gauge" in this way, we can rewrite the partition function (2.4) as where ∆(α) is the van der Monde determinant. The factor ∆(α) 2 in (2.10) appears from the Fadeev-Popov procedure for the gauge fixing, and it acts as a repulsive potential between the eigenvalues α i of A 0 . Here we consider a situation in which the eigenvalues of A 0 are well separated from each other. Then the action S = S b + S f can be expanded as 12) omitting the subleading terms for large |α a − α b |. Integrating out A i at one loop neglecting the zero modes corresponding to diagonal elements, we obtain ∆(α) −18 . On the other hand, integrating out Ψ α at one loop similarly, we obtain ∆(α) 16 . Thus we find that the potential between α i is canceled exactly at the one-loop level. This is actually a consequence of supersymmetry [5] of the model (2.4). Owing to this property, the eigenvalue distribution of A 0 would extend to infinity even for finite N if the cutoff (2.7) were absent.

JHEP11(2015)070
In fact, after some manipulation and rescaling of A µ , we can rewrite the partition function (2.4) as [44] where θ (x) is the Heaviside step function. This form allows us to perform Monte Carlo simulation of the Lorentzian model without the sign problem unlike the Euclidean model. 1 This is surprising considering that the original model (2.4) has a severe sign problem (even without fermions) since the bosonic part of the action (2.2) appears in the integrand of (2.4) as a pure phase. The basic idea in reformulating the model as (2.14) is to integrate out first the fluctuations A µ → ρA µ associated with the overall scale factor ρ. Since the bosonic action S b is a homogeneous polynomial in A µ of the fourth order, it behaves as S b → ρ 4 S b under the scale transformation. Note also that the Pfaffian as well as the measure dA simply gives a power of ρ under the scale transformation, which allows us to perform the integration over ρ first. For sufficiently large Λ, the result of ρ integration vanishes due to violent fluctuations of the phase factor e iρ 4 S b unless S b = 0, and hence we obtain the first delta function in (2.14). In order to perform the remaining integral, one needs to impose a constraint in order not to doubly integrate over the scale factor. This is realized by imposing a constraint 1 N Tr (A i ) 2 = 1, which is compatible with the cutoff (2.8). This gives the second delta function in (2.14). See appendix A of ref. [57] for a more detailed derivation of (2.14).
It turns out that one can extract a time-evolution from configurations generated by simulating (2.14). A crucial observation is that the spatial matrices A i have a band-diagonal structure in the SU(N ) basis in which A 0 has the diagonal form (2.9). More precisely, there exists some integer n such that the elements of spatial matrices (A i ) ab for |a − b| > n are much smaller than those for |a − b| ≤ n. Based on this observation, we may naturally consider n × n matrices as representing the state of the universe at time t, where I, J = 1, . . . , n and ν = 0, 1, . . . , N − n. The time t in (2.15) is defined by corresponding to the n × n matricesĀ i . For example, we can define the extent of space at time t as

JHEP11(2015)070
where the symbol tr represents a trace over the n × n block. We also define the "moment of inertia tensor" which is a 9 × 9 real symmetric matrix. The eigenvalues of T ij (t), which we denote by represent the spatial extent in each of the nine directions at time t. Note that the expectation values λ i (t) tend to be equal in the large-N limit if the SO(9) symmetry is not spontaneously broken. This is the case at early times of the time-evolution. After a critical time t c , on the other hand, we find that three largest eigenvalues λ i (t) (i = 1, 2, 3) become significantly larger than the others, which implies that the SO(9) symmetry is spontaneously broken down to SO(3). It would be interesting to study a long time-evolution of the model and see how the expansion of space proceeds. This requires very large matrices, which makes the simulation unfeasible. In the previous work [57], we studied a simplified model, in which the Pfaffian is replaced by the one-loop contribution ∆(α) 16 mentioned above. This replacement is expected to be valid at early times, where the expansion of space has not proceeded much and the leading term in (2.13) is indeed dominant. According to the argument below (2.13), the potential between the eigenvalues of A 0 is canceled at one loop and hence the cutoff (2.7) in the temporal direction is needed in this simplified model as well as in the original model. On the other hand, this simplified model can be simulated with much less effort than the original model. 2 In ref. [57] the (5+1)D version of the simplified model was studied with the matrix size N ≤ 64, and the SO(5) symmetry was found to be broken spontaneously down to SO(3) at some point in time similarly to the original model. Moreover, the expanding behavior of the 3D space turned out to be exponential, 3 and no tendencies of slowing down were observed within the scaling region. Analogous behaviors were also confirmed for the (9+1)D version of the simplified model. In the original model, on the other hand, the subleading term in the fermionic action (2.13) becomes important at late times as the expansion proceeds, and hence it can affect the expanding behavior.

The bosonic Lorentzian type IIB matrix model
In this paper we study a bosonic model in which the fermionic matrices are simply omitted. The partition function is given by (3.1)

JHEP11(2015)070
In section 2 we reviewed an argument for the necessity of the temporal cutoff in the original model and the simplified model for early times. In the present case of the bosonic model (3.1), the same argument suggests that there is no need to introduce the temporal cutoff (2.7); the eigenvalues of the temporal matrix A 0 attract one another at the one-loop level in the absence of supersymmetry. Note, however, that we still need a cutoff (2.8) in the spatial direction since the bosonic action (2.2) is not positive semi-definite unlike in the Euclidean model, and therefore it does not suppress configurations with large A i . Corresponding to (2.14), we can study the bosonic Lorentzian type IIB matrix model by simulating which requires computational efforts comparable to the simplified model for early times reviewed in the previous section. We have used a large-scale parallel computer to simulate the model (3.3) with the matrix size up to N = 512, which enables us to investigate a long time-evolution. See appendix A for the details of the simulation. We measure the quantity 1 N Tr (A 0 ) 2 , which represents the extent of the eigenvalue distribution of A 0 . As we mentioned above, this quantity turns out to be finite in the model (3.3) although we do not introduce a cutoff in the temporal direction such as (2.7). In figure 1 (Left) we plot the results against N . At small N , it is almost independent of N . However, for N ≥ N c = 112, it begins to increase linearly with N . Figure 1 (right) shows the N dependence of the expectation values λ i (t) of the nine eigenvalues of T ij (t) evaluated at t = t peak , where R 2 (t) becomes maximum. 4 For small N , there is no significant gap between the nine eigenvalues, whereas for N ≥ N c , we observe a big gap between λ 3 (t peak ) and λ 4 (t peak ) . We will see in section 5 that the SO(9) symmetry is broken down to SO(3) after a critical time similarly to the original Lorentzian type IIB matrix model.

Properties of the bosonic model for N < N c
In this section we discuss the properties of the bosonic model for N < N c . In order to extract the time-evolution, we need to determine the block size n to be used in eq. (2.15). In figure 2 (Left) we plot the magnitude of the off-diagonal elements of A i against the time separation α a − α b for N = 110. The origin in the horizontal axis corresponds to the diagonal elements. We observe a nice scaling behavior for all the matrix elements. However, the magnitude falls off rather smoothly as one goes in the off-diagonal direction, which means that the dominant matrix configurations do not have a band-diagonal structure.
In this situation, we cannot naturally define the block matrices (2.15) representing the state at each time and hence the notion of time-evolution becomes obscure. Let us JHEP11(2015)070  nevertheless try to extract the "time-evolution" using n = 14 as the block size, which is the value obtained for N = N c = 112 in the way described in the next section. In figure 2 (Right) we plot the expectation values λ i (t) for N = 110. It turns out that there is only little t-dependence, and there is no clear gap between the eigenvalues for all t. The situation for smaller N is similar to the N = 110 case. In figure 3 we plot the extent of space R 2 (t) as a function of t for N = 64, 96 and 110 obtained with the same block size n = 14. The dependence on N turns out to be modest. 2), which is not spontaneously broken. However, the reason why the time evolution turns out to be an expansion and a contraction of the "universe" is not obvious. We can only say that it appears as a consequence of the path integral (3.2).
In order to study the large-N scaling property, we perform simulation for N = 256, 384, 512 as well. In figure 5 (Top-Left) we zoom up the region near the origin in figure 4 (Left). From this figure, we determine the block size for N = 128 to be n = 20. Similarly, from the other figures in figure 5, we determine the block size for N = 256, N = 384 and 5 The fact that the spatial dimensionality after the spontaneous symmetry breaking turned out to be the same as in the original model is understandable from the view point of the mechanism suggested in ref. [44], which involves only the boson part of the action.  512 to be n = 24, n = 28 and 32, respectively. The fact that the block size n increases with N is natural since the block matrices representing the state of the universe at fixed time should have infinitely many degrees of freedom in the large-N limit. The definition of the critical time t c is ambiguous at finite N . See figure 6 (Left), where we plot the expectation values λ i (t) of the eigenvalues of T ij (t) against t for N = 512. Note first that the appearance of a gap between λ 3 (t) and λ 4 (t) signals the spontaneous symmetry breaking of SO(9) to SO(3). Let us therefore define the separation d j (t) = λ j (t) − λ j+1 (t) . Then we find that the symmetric phase can be characterized by d 1 (t) > d 2 (t) > · · · > d 8 (t), while in the broken phase we find d 2 (t) < d 3 (t). Therefore we may define the critical time t c by the largest value of t such that d 1 (t) > d 2 (t) > · · · > d 8 (t) holds for t ≤ t . For instance, the critical time t c obtained in this way for N = 512 from figure 6 (Left) is t c = −0.76559 (7). Similarly we obtain t c = −0.76930 (7) for N = 384. Applying the same procedure to smaller N , we find that the large-N scaling behavior in figure 6 (Right) becomes less clear due to finite N effects. We absorb these finite N effects by adjusting the value of t c slightly 6 from the one determined from the above procedure. As is proposed in the original Lorentzian type IIB matrix model [44], we use the extent of space R(t c ) at the critical time to fix the scale. Explicit values of R(t c ) are given in table 1 together with the block size n and the critical time t c for each N .
In figure 6 (Right) the extent of space R 2 (t) is plotted against t. The large-N scaling behavior is observed 7 by shifting the time coordinate so that the critical time comes to the origin and by plotting dimensionful quantities in units of R(t c ). The observed large-N scaling shows that the theory one obtains in the large-N limit is characterized by one scale parameter R(t c ) and it does not contain any dimensionless parameters.
It turns out that the behavior of R 2 (t) at t > t c can be fitted to an exponential function only for a finite range, and it seems to slow down at later times in contrast to the JHEP11(2015)070 exponential behavior observed in ref. [57] in a simplified model for early times. We can actually fit our N = 512 data within 1.85 ≤ (t − t c )/R(t c ) ≤ 2.5 to a straight line, which corresponds to the power-law expansion Whether this behavior can be confirmed in the scaling region at larger N remains to be seen. In appendix B we present the results for the (5+1)D version of the bosonic type IIB matrix model. While we observe qualitatively the same behaviors, there are also some interesting quantitative differences. In order to understand the observed large-N scaling further, we investigate how the continuum limit and the infinite volume limit in the temporal direction are achieved in the large-N limit. Here we restrict ourselves to N ≥ 256 since N = 128 is too close to the critical value N c = 112. As the "lattice spacing" in the temporal direction, we consider the separation of data points in figure 6 (Right) in the horizontal direction. This quantity is actually t-dependent, and it can be defined more explicitly as δt R(tc) , where δt is the difference of (2.16) between adjacent blocks. In figure 7 (Left) we plot this t-dependent "lattice spacing", choosing the horizontal axis to be the same as in figure 6 (Right). We find clear tendency that the "lattice spacing" at the same point on the horizontal axis decreases as N increases. As the "volume" in the temporal direction, we define Using this quantity, we can also define an average "lattice spacing" ε = ∆/ν, where ν is the number of data points within the region [t c , t peak ]. The values of ε and ∆ obtained for each N are given in table 1. We find that the average "lattice spacing" ε decreases and the "volume" increases as N becomes large. In figure 7 (Right) we plot ε and ∆ against N in the log scale. The straight lines represent fits to the power-law behaviors, although the behaviors may be subject to a qualitative change at larger N . In particular, it is an interesting dynamical question whether ∆ → ∞ or ∆ → const. in the large-N limit. In the former case, the expansion of space continues forever, since the time t peak at the peak cannot be reached within a finite time. In the latter case, on the other hand, the space stops expanding in a finite time and starts to shrink. By addressing this issue in the original model, one can, in principle, predict the fate of our Universe.

Summary and discussions
In this paper we have studied a bosonic model, which can be obtained from the Lorentzian type IIB matrix model by omitting the fermionic matrices. Due to the attractive potential between the eigenvalues of A 0 arising from integrating out the spatial matrices at one loop, the eigenvalue distribution of A 0 has a finite extent even if one does not introduce the temporal cutoff (2.7). In the original supersymmetric model, this attractive potential is canceled by the repulsive potential arising from integrating out the fermionic matrices at one loop and from the van der Monde determinant. The simplified model for early times  Table 1. The block size n, the critical time t c and the extent of space R(t c ) at the critical time, which are used to make the plot in figure 6 (Right). We also present the explicit values of the average "lattice spacing" ε and the "volume" ∆ in the temporal direction, which are plotted in figure 7 (Right). studied in ref. [57] has the same feature since the Pfaffian is replaced by the one-loop contribution. Note that this simplified model as well as the original supersymmetric model allows the extraction of a sensible time-evolution. We may naturally speculate that this is possible only when the eigenvalues of A 0 can fluctuate to large values as in these models. Indeed for N < N c , we observe that the extent of the eigenvalue distribution of A 0 is almost independent of N , and the dominant matrix configurations do not allow the extraction of a sensible time-evolution. On the other hand, for N ≥ N c , we find that the extent of the eigenvalue distribution of A 0 grows linearly in N , and a sensible timeevolution can be extracted. We find that the SO(9) symmetry is spontaneously broken down to SO(3) at some critical time similarly to the original model [44]. The quantity such as R(t) shows a clear large-N scaling behavior. The growth of R(t) at late times turns out to be consistent with R(t) ∝ t 1/2 , which is reminiscent of the behavior of the Friedmann-Robertson-Walker universe in the radiation dominated era.
Let us recall that in the simplified model for early times, the growth of R(t) was observed to be exponential [57]. In that model, only the first term in (2.13) was used to represent the effect of fermionic matrices. We consider that the exponential expansion occurs also in the original supersymmetric model at early times as is suggested by direct

JHEP11(2015)070
Monte Carlo studies up to N ≤ 24 [58]. At late times, however, the subleading term in (2.13) becomes important due to the expansion of space, and that would affect the expanding behavior. Note that the repulsive potential for the eigenvalues of A 0 is obtained from integrating out the fermionic matrices without the subleading term. Therefore one of the effects of the subleading term would be to make the repulsive potential less effective. Considering that the bosonic model mimics such a situation, we speculate that the exponential expansion in the original model changes into a power-law expansion at some point in time, where the subleading term in (2.13) becomes important. According to this scenario, the number of e-foldings is determined dynamically in the Lorentzian type IIB matrix model. It would be interesting to confirm the transition directly by simulating the original supersymmetric model. An attempt in doing this with a systematic approximation is in progress.

Acknowledgments
We would like to thank K. Anagnostopoulos

A Details of Monte Carlo simulation
In this appendix we give the details on how we perform Monte Carlo simulation of the bosonic model (3.3).
First the delta functions in (3.3) are replaced by Gaussian potentials as where the coefficients γ C and γ L are taken large enough to fix each observable to the specified value. The values used in actual simulation are given in table 2.
Another important issue we have to take care of is the spontaneous breaking of the shift symmetry A 0 → A 0 + α1. For instance, let us consider calculating the expectation value R 2 (t) defined in (2.17). The peak of this quantity measured for each configuration fluctuates considerably. This reflects the ambiguity in choosing the origin of the time coordinate, and we should fix it before taking the ensemble average. Here we fix it by JHEP11(2015)070 introducing a potential where the values of the coefficient γ sym used in our simulation are given in table 2.
To summarize, the model we simulate is given by The simulation of the model (A.5) can be performed by using the Hybrid Monte Carlo (HMC) method. First we rewrite the model by introducing auxiliary variables p a and (X i ) ab (a, b = 1, . . . , N ) with the action Here p a are real variables, whereas X i are traceless Hermitian matrices. We update all the variables in the model (A.6) in the following way. First we regard p a as the conjugate momenta of α a and X i as the conjugate momenta of A i . Then we regard S HMC in (A.6) as the Hamiltonian H and solve the classical equations of motion obtained as the Hamilton equations for some fictitious time τ . This part of the algorithm is called the Molecular Dynamics. In order to solve the Hamilton equations (A.7) numerically, we discretize them using the so-called leap-frog discretization, which maintains reversibility with respect to τ . Starting from the previous configuration at τ = 0, we obtain a new configuration at τ = τ f by solving (A.7) with the step size ∆τ so that τ f = N τ · ∆τ , where N τ is the number of steps. We accept the new configuration with the probability min (1, exp (−∆S HMC )), where ∆S HMC ≡ S HMC (τ f ) − S HMC (0), following the idea of the Metropolis algorithm to satisfy the detailed balance. The important point here is that S HMC is nothing but the Hamiltonian H, which is preserved in the classical dynamics if the equations (A.7) are solved exactly. In fact, ∆S HMC becomes non-zero due to the discretization, but it is guaranteed to be a small quantity of the order of (∆τ ) 2 . By choosing sufficiently small ∆τ , the acceptance JHEP11(2015)070  Table 2. The values of the parameters γ C , γ L and γ sym in (A.5) used in our simulation. We also give the values of the parameters in the HMC algorithm: the number of steps N τ in the Molecular Dynamics and its step size ∆τ . In the last column, we give the number of "trajectories", which represents how many times we solve the Molecular Dynamics after thermalization to achieve the statistics of our data.
rate can be kept reasonably high, which enables the system to move around efficiently in the configuration space. Note also that the auxiliary variables p a and (X i ) ab appear only as the Gaussian terms in (A.6). Therefore, we can update them independently by using normalized Gaussian random numbers. This procedure of refreshing the conjugate momenta should be done each time we start a Molecular Dynamics procedure.
To summarize, the HMC algorithm as applied to our system can be described as follows.
1. Generate initial configurations of p a (0) and X i (0) with Gaussian distribution ∝ e In the HMC algorithm, there are two parameters 8 ∆τ and τ f . In the present work we choose them as in table 2.

B Results for the (5+1)D version of the bosonic model
In this appendix we present our results 9 for a bosonic model that can be obtained by omitting fermionic matrices in the (5+1)D version of the type IIB matrix model. The latter model is obtained formally by dimensional reducing the 6D N = 1 super Yang-Mills theory to a point, and it consists of six bosonic matrices A µ (µ = 1, . . . , 6) and four fermionic matrices Ψ α (α = 1, . . . , 4) representing four components of a 6D Weyl spinor. The form of the bosonic part of the action is the same as that of the original type IIB matrix model, which is given in (2.2). 8 These parameters can be optimized as follows. For fixed τ f , it is optimal to choose ∆τ so that ∆τ × (acceptance rate) is maximized. Then τ f can be optimized to minimize the auto-correlation time in units of one step in the Molecular Dynamics. 9 Preliminary results shown in figure 9 (Right) are published in the proceedings [58].    Table 4. The block size n, the critical time t c , the extent of space R(t c ) at the critical time, which are used in the (5+1)D model to make the plot in figure 9 (Right). We also present the explicit values of the average "lattice spacing" ε and the "volume" ∆ in the temporal direction, which are plotted in figure 10 (Right).
In figure 8 (Left), we plot the extent 1 N Tr (A 0 ) 2 of the eigenvalue distribution of A 0 against N . In figure 8 (Right), we plot the expectation values λ i (t) of the five eigenvalues of T ij (t) at t = t peak against N . While the qualitative behaviors are the same as in the (9+1)D case shown in figure 1, we find that the critical N c is smaller and the slope of the linearly increasing behavior of 1 N Tr (A 0 ) 2 for N ≥ N c is larger. We can understand this difference by considering the attractive potential between the eigenvalues of A 0 discussed below eq. (2.12). For general spatial dimensionality d, one obtains a factor ∆ −2d (α) from integrating out the spatial matrices A i at one loop. This factor acts as an attractive potential between the eigenvalues of A 0 , and it is stronger for larger d. Below we discuss the properties of the (5+1)D model for N ≥ N c . (The parameters used in the simulation are listed in table 3.) We have determined the block size to be n = 8, 10, 12 for N = 64, 96, 128, respectively, from the fall-off of the off-diagonal elements of A i as is done for the (9+1)D case in section 5. In figure 9 (Left) we plot the expectation values λ i (t) of the five eigenvalues of T ij (t) for N = 128, which shows that the SO(5) symmetry is broken spontaneously down to SO(3) after a critical time. From this kind of figures, we can determine the critical time t c for each N as described 10 in section 5.

JHEP11(2015)070
In figure 9 (Right) we show the large-N scaling behavior of the extent of space R 2 (t). Explicit values of R(t c ), together with the block size n and the critical time t c , which are used to make this plot, are given in table 4. The power-law expansion (5.1) is observed at late times similarly to the (9+1)D model.
In figure 10 (Left) we plot the t-dependent "lattice spacing", which shows how the continuum limit is achieved as N increases. The average "lattice spacing" ε and the "volume" ∆ in the temporal direction are given in table 4. In figure 10 (Right) we plot them in the log scale. The straight lines represent fits to the power-law behaviors.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.