A renormalization group method for studying the early universe in the Lorentzian IIB matrix model

We propose a new method for studying the early universe in the Lorentzian version of the IIB matrix model, which is considered to be a nonperturbative formulation of superstring theory. This method is based on the idea of renormalization group, and it enables us to study the time-evolution of the universe for much longer time than in the previous work, which showed that the SO(9) rotational symmetry is spontaneously broken down to SO(3) after a"critical time". We demonstrate how this method works in a simplified model, which is expected to capture the behaviors of the original model when the space is not so large. In particular, we present clear evidence that the three-dimensional space expands exponentially after the critical time in this simplified model.


Introduction
Understanding how our universe began is one of the most fundamental and fascinating themes in theoretical physics. For instance, there are good reasons to believe that our universe underwent a rapid expansion called inflation before the Big Bang. While there are various models which describe inflation phenomenologically, we still do not have a description based on first-principle calculations in a complete quantum gravity theory like superstring theory. The most crucial problem with superstring theory is that it is defined only perturbatively around consistent backgrounds, and within such perturbative formulations, it is known that the cosmic singularity is not resolved generally [1][2][3][4]. In order to overcome this problem, one really needs to use a fully nonperturbative formulation. In fact, there exist concrete proposals for such a formulation using supersymmetric matrix models [5][6][7]. These models can be obtained formally 1 by dimensionally reducing tendimensional N = 1 super Yang-Mills theory to d = 0, 1, 2 dimensions, respectively. Based on these proposals, various issues of the early universe have been discussed [8][9][10][11][12][13][14][15][16][17]. As a closely related direction, ref. [18] proposes a conformal field theory, which is holographically dual to inflationary models. (See also ref. [19] and references therein.) The IIB matrix model [5] is one of the matrix model proposals corresponding to the d = 0 case mentioned above, in which not only space but also time should emerge dynamically from the matrix degrees of freedom. This aspect of the IIB matrix model has also been discussed intensively as "emergent gravity" [20][21][22][23][24][25][26][27] in noncommutative field theories, which appear in the IIB matrix model for a particular class of classical backgrounds [28][29][30][31]. (See also ref. [32] for a different proposal for emergence of curved space-time in the matrix model.) Until quite recently, the IIB matrix model was studied after making a Wick rotation since the partition function of the Euclidean matrix model obtained in this way was shown to be finite [33,34]. In fact, a lot of efforts have been devoted to identifying the matrix configurations that dominate the partition function using various methods [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49]. However, the Euclidean matrix model is clearly not applicable to cosmology since it does not provide the real-time dynamics. Moreover, it is known that the Wick rotation is more subtle in quantum gravity than in quantum field theory at the nonperturbative level (See, for instance, refs. [50,51]). Indeed a recent study based on the Gaussian expansion method suggests that the space-time obtained dynamically in the Euclidean matrix model does not seem to correspond to our four-dimensional space-time [49].
Motivated by all these problems with the Euclidean IIB matrix model, three of the authors (S.-W.K., J.N. and A.T.) studied the Lorentzian version of the IIB matrix model by Monte Carlo simulation for the first time [52]. Unlike the Euclidean version, one has to introduce infrared cutoffs in both spatial and temporal directions in order to make the partition function finite. However, it was found that these two cutoffs can be removed in the large-N limit in such a way that physical quantities scale. The eigenvalue distribution of the matrix representing the time extends in that limit, and the dominant matrix configurations have a very nontrivial structure, which enables us to naturally extract the time-evolution. Quite surprisingly, it was found that a phase transition occurs at some point in time, and after that, only three out of nine spatial directions start to expand. This phase transition can be interpreted as the birth of our 3-dimensional universe in superstring theory. It should be emphasized that the results seem to suggest that the space-time dimensionality is determined uniquely by the nonperturbative dynamics of superstrings unlike in perturbative string theory, in which consistent backgrounds can have various space-time dimensionality.
As another important property of the Lorentzian IIB matrix model, it is expected that the classical approximation becomes valid at late times [53,54]. The reason for this is that each term in the action has large contribution from the degrees of freedom at late times due to the expansion of the universe. One can actually construct a simple solution representing an expanding (3+1)-dimensional universe, which naturally solves the cosmological constant problem [54]. It has also been argued that local field theory emerges from low-lying fluctuation modes around a solution representing a commutative spacetime [55]. In fact, the classical equations of motion of the matrix model have infinitely many solutions, which is reminiscent of the so-called landscape in superstring theory. Unlike the situation with the landscape, however, there is a definite criterion to pick up a particular solution describing the late-time behaviors since we have a well-defined partition function.
Clearly it is important to extend the Monte Carlo studies in ref. [52] to much longer time. For instance, it would be interesting to see whether the inflation and the Big Bang occurs in this model as is generally believed in modern cosmology. Moreover, if we can go further and reach the time region in which the classical approximation is valid, we should be able to determine the solution which actually describes the late-time behaviors in the dominant configurations. Once this has been done, we should be able to derive the effective field theory below the scale where gravity decouples by considering the fluctuations around the classical solution. In particular, it would be interesting to see whether the Standard Model particles appear at low energy, for instance, in a way speculated in refs. [56][57][58].
To this end, we develop a new method that enables us to investigate a long timeevolution of the universe in the Lorentzian IIB matrix model by Monte Carlo simulation. Note, in particular, that the time scale that one would hope to achieve is, at least, a few orders of magnitude larger than the typical time scale of the model. If one attempts to study it directly, one would need a huge matrix size, which makes the calculation impractical. In this paper we show that there exists a "renormalized theory", which corresponds to a theory obtained from the original model by integrating out the dynamical degrees of freedom at earlier times. The renormalized theory has two extra parameters compared with the original model, which can be used to optimize the length of the time region that one can actually probe. By simulating the renormalized theory with optimized parameters, we can investigate the late-time behaviors with much less dynamical degrees of freedom than the direct approach would require.
In order to show how the method works, we consider a simplified model, which can be obtained from the original model by neglecting the coupling of fermionic matrices to the spatial bosonic matrices. This approximation is expected to be valid at early times, where the space is not so large. Then the Pfaffian that arises from integrating out fermionic matrices can be expressed by some power of the van der Monde determinant, which is written explicitly in terms of the eigenvalues of the temporal matrix only, and the simulation becomes as fast as the bosonic model. The simplified model indeed retains the important properties of the original model such as the spontaneous breaking of the rotational symmetry at some critical time. Moreover, we find that the size of the universe grows exponentially after the critical time. We apply the renormalization group method to the simplified model and confirm the exponential expansion more clearly with smaller matrices.
The rest of the paper is organized as follows. In section 2 we review some important properties of the Lorentzian IIB matrix model. In section 3 we define the simplified model, and present results obtained by direct Monte Carlo studies. In particular, we show that the exponential expansion is realized in this model after the spontaneous breaking of rotational symmetry. In section 4 we describe the renormalization group method. In section 5 we apply the method to the simplified model and show that it allows us to study the timeevolution much more efficiently. Section 6 is devoted to a summary and discussions. In appendix A we give the details of our Monte Carlo simulation.

Brief review of the Lorentzian IIB matrix model
The IIB matrix model has an action [5]

1)
2) where the bosonic N × N matrices A µ (µ = 0, · · · , 9) and the fermionic ones Ψ α (α = 1, · · · , 16) are both traceless and Hermitian. Γ µ are 10-dimensional gamma-matrices after the Weyl projection and C is the charge conjugation matrix. Since the coupling constant g can be absorbed by rescaling A µ and Ψ appropriately, it is merely a scale parameter. The IIB matrix model is conjectured to be a nonperturbative definition of superstring theory [5]. There are various pieces of evidence for this conjecture. First of all, the action (2.1) can be regarded as a matrix regularization of the worldsheet action of type IIB superstring theory in the Schild gauge [5]. 2 Secondly, D-branes in type IIB superstring theory can be described in the matrix model, and the interaction between them can be correctly reproduced [5]. Thirdly, under a few reasonable assumptions, the string field Hamiltonian for type IIB superstring theory can be derived from Schwinger-Dyson equations for the Wilson loop operators, which are identified as creation and annihilation operators of strings [59].
In all these connections to type IIB superstring theory, the target space coordinates are identified with the eigenvalues of the matrices A µ . In particular, this identification is consistent with the supersymmetry algebra of the model, in which the translation that appears from the anti-commutator of supersymmetry generators is identified with the shift symmetry A µ → A µ + α µ 1 of the model, where α µ ∈ R. Also the fact that the model has extended N = 2 supersymmetry in ten dimensions is consistent with the fact that the model actually includes gravity since it is known in field theory that N = 1 supersymmetry is the maximal one that can be achieved in ten dimensions without including gravity.
The partition function for the Lorentzian version of the IIB matrix model is proposed as [52] Z = dA dΨ e iS . (2.4) Integrating out the fermionic matrices, we obtain the Pfaffian which is real. Note that the bosonic action (2.2) can be written as where we have defined Hermitian matrices F µν = i [A µ , A ν ]. Hence the bosonic action is not positive semi-definite.
The partition function (2.4) is not finite as it stands, but it can be made finite by introducing infrared cutoffs in both temporal and spatial directions as [52] It turned out that these two cutoffs can be removed in the large-N limit in such a way that physical quantities scale. The resulting theory thus obtained has no parameter except one scale parameter. After some manipulation, the partition function can be rewritten into the form [52] where θ(x) is the Heaviside step function. Since the Pfaffian PfM(A) is real unlike in the Euclidean case [39][40][41]48], the model (2.10) can be studied by Monte Carlo simulation without the sign problem. 3 Let us then discuss how we can extract the time-evolution from a configuration generated by simulating the system (2.10). First we choose the SU(N ) basis in such a way that the temporal matrix A 0 is diagonalized as (2.11) In that basis, it turned out that the spatial matrices A i has a band-diagonal structure with off-diagonal elements (A i ) ab for |a − b| ≥ n being small for some integer n. Therefore, we may naturally consider n × n block matrices where I, J = 1, · · · , n and ν = 0, 1, · · · , N − n, as representing a state of the universe at the time For instance, the extent of space at time t is defined by where the trace here is taken over the n × n block. In order to see the spontaneous breaking of SO(9) symmetry, we define the "moment of inertia tensor" which is represented by a real symmetric 9 × 9 matrix. We denote the real positive semidefinite eigenvalues of T ij (t) as λ j (t) with the ordering If the SO(9) symmetry is not spontaneously broken, the expectation values λ i (t) become equal in the large-N (and large-n) limit. We find that this is indeed the case at early times, while at sufficiently late times, three of the eigenvalues become considerably larger than the others, suggesting that the SO(9) symmetry is spontaneously broken down to SO(3) after a critical time.
The necessity for introducing the cutoff (2.8) in the temporal direction can be understood as follows. Let us consider a situation in which the eigenvalues of A 0 are well separated from each other and estimate the effective action for the eigenvalues perturbatively. By fixing the gauge to (2.11), we rewrite the integration over A µ as where ∆(α) is the van der Monde (VDM) determinant. The action can be expanded as where the omitted terms are subleading for large |α a − α b |. Integrating out A i at the one-loop level, one obtains ∆(α) −18 neglecting the zero modes corresponding to diagonal elements. Integrating out Ψ α at the one-loop level, one obtains ∆(α) 16 neglecting the zero modes. Thus one finds that the ∆(α) 2 in (2.17) is canceled exactly at the one-loop level, which is actually a consequence of the supersymmetry [5] of the model (2.4). Due to this property, the eigenvalue distribution of A 0 extends to infinity even for finite N if it were not for the cutoff (2.8).

A simplified model and its properties
The argument given above motivates us to generalize the Lorentzian IIB matrix model to (d + 1)-dimensional versions (d = 9, 5, 3), which can be obtained by dimensional reduction of (d + 1)-dimensional N = 1 super Yang-Mills theory. (The d = 9 case corresponds to the Lorentzian IIB matrix model.) In general, integration over A i gives ∆(α) −2d , while integration over the fermionic matrices gives ∆(α) 2(d−1) . Hence the VDM determinant that appears as in (2.17) is canceled exactly at the one-loop level and there is no interaction among the eigenvalues of A 0 at the one-loop level. In Monte Carlo simulations, the eigenvalue distribution indeed extends to infinity if one does not introduce the temporal cutoff κ as in (2.8). If one omits fermions, one obtains an attractive force between the eigenvalues of A 0 , and the eigenvalue distribution of A 0 has a finite extent without any cutoff. In fact, the d = 5 supersymmetric model turns out to have very similar properties as the original model. 4 In particular, the SO(5) rotational symmetry is broken spontaneously down to SO(3) after a critical time. The (5+1)-dimensional model contains bosonic matrices A µ (µ = 0, · · · , 5) and fermionic matrices Ψ α andΨ α (α = 1, · · · , 4). The action for the fermionic matrices is given by where Γ µ are 6-dimensional gamma-matrices after the Weyl projection. Integrating out the fermionic matrices, we obtain the determinant which is real. This model can be studied by Monte Carlo simulation using the partition function (2.10), where PfM(A) should be replaced by det M(A).
In performing Monte Carlo simulation of the model (2.10) or its (5+1)-dimensional version, the most time-consuming part comes from calculating the contribution from the fermions. Here we consider a simplified model, which can be obtained by replacing the Pfaffian (2.5) or the determinant (3.2) by ∆(α) 2(d−1) , which we obtained as the leading contribution when the separation |α a − α b | of the eigenvalues of A 0 is large. Note that this amounts to neglecting the terms proportional to the spatial matrices A i (i = 1, . . . , d) in the fermionic action (2.3) or (3.1). Therefore we expect that the simplified model captures the qualitative behaviors of the original models at early times before the expansion of space proceeds much. Thus we arrive at the model where A 0 is given by (2.11). This model, which we call the VDM model in what follows, can be simulated as easily as the bosonic model. Moreover, it shares important properties with the original supersymmetric models such as the spontaneous symmetry breaking of SO(d) and expanding behavior of the three-dimensional space after the critical time.
In this paper we study the VDM model in the d = 5 case 5 for simplicity by Monte Carlo simulation. (See appendix A for the details of our simulation.) We set L = 1 in (3.3) 4 Similarity of (5+1)-dimensional version and the (9+1)-dimensional version is seen also in the Euclidean IIB matrix model. It was found that SO(D) symmetry of the D-dimensional model is broken down to SO(3) symmetry for D = 10 [49] and D = 6 [60], and various properties associated with the SSB turned out to be common to both models [49]. 5 The properties of the d = 5 VDM model observed here are confirmed also in the d = 9 case. In particular, we observe the SSB from SO(9) to SO(3) at some critical time.  without loss of generality since it only fixes the scale of the model. In figure 1 (Left) we plot the eigenvalues α a of A 0 with the ordering (2.11) against its label a for N = 64 and κ = 4. Figure 1 (Right) shows the magnitude of the off-diagonal elements of A i against the time separation α a − α b . We find that the off-diagonal element decreases rapidly as one goes away from the diagonal element. Moreover, we observe a nice scaling behavior for sufficiently large |α a − α b |. The region with small |α a − α b | that does not scale includes roughly 8 points. Based on this observation, we choose the block size to be n = 8 in this paper.
In figure 2 (Left) we plot R 2 (t) against t for N = 64 and κ = 4.0. This plot shows that the space starts to expand at a critical time. In figure 2 (Right) we plot the expectation value λ i (t) of the five eigenvalues of T ij (t), which shows that the SSB from SO(5) to SO(3) occurs at the critical time.
The definition of the critical time t c is ambiguous at finite N . As a convenient choice we define it as follows. Note first that the appearance of a gap between λ 3 (t) and λ 4 (t) signals the SSB of SO(5) 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 , while in the broken phase we find d 2 (t) < d 3 (t). Therefore we define the critical time t c by the largest value of t ′ such that For instance, the critical time t c obtained in this way from figure 2 (Right) is t c = −0.8813(2) and the extent of space at the critical time is R 2 (t c ) = 0.139(1).
In figure 3 (Left) we plot R 2 (t)/R 2 (t c ) against (t − t c )/R(t c ) for various values of κ and N , which reveals a nice scaling property 6 of the function R 2 (t). The shift in time is necessary since only the difference (t − t c ) is meaningful. We also normalize all dimensionful quantities by R(t c ), which represents the size of the universe when it was born. Interestingly, we observe that our data can be fitted well to (We fix the second coefficient to beC = 1 − C using the constraint f (0) = 1, which follows from the definition of f (x).) This is demonstrated in figure 3 (Right), where we plot R 2 (t)/R 2 (t c ) − C against (t − t c )/R(t c ) in the logarithmic scale. This exponential growth is reminiscent of the inflation, which is expected to have taken place in the early universe. A similar behavior at early times can be seen also in our preliminary results for the Lorentzian IIB matrix model [61].
Let us discuss how we should take the large-N limit. For that we define the "lattice spacing" ǫ and the "time-extent" ∆ by where δt is the mean separation of the eigenvalues of A 0 and t p represents the time t at which the extent of space R(t) becomes maximum. (In fact, t p = 0 as one can see from figure 2 due to the time reflection symmetry.) Results for different κ and N correspond to different ǫ and ∆. As κ is increased for a fixed N , the time-extent ∆ increases and one can see late-time behaviors more. However, the lattice spacing ǫ increases at the same time, which results in deviations from the scaling behavior due to "lattice artifacts". Therefore, one needs to increase N as one increases κ to see the scaling behavior at later times. The fact that the scaling behavior extends with increasing N implies that the two cut-offs (2.8) and (2.9) can be removed in the large-N limit. Whether the time-extent ∆ diverges in the large-N limit or not is an interesting dynamical question. If it diverges, the t > 0 region in figure 2, for instance, becomes invisible in the large-N limit, hence there will be no Big Crunch. In order to study the late-time behaviors, we need to increase the matrix size further. However, we notice from figure 3 (Left) that the symmetric phase extends more than the broken phase with increasing N . Due to this property of the model, it is not efficient to study the late-time behaviors by just increasing the matrix size.

Renormalization group method
In this section we propose a new method based on the idea of renormalization group, which enables us to study the late-time behaviors much more efficiently than in a direct approach. Note first that the late-time behaviors are described by the inner part of the matrices A µ (See figure 4.) if we fix the gauge to (2.11). The corresponding degrees of freedom are given byÑ ×Ñ Hermitian matricesÃ µ , which are defined by where the indices a and b run from 1 toÑ. In principle, we can derive the renormalized theory forÃ µ by integrating out the other degrees of freedom in the original matrices A µ . Once we know the form of the renormalized theory, we can study the late-time behaviors efficiently by simulating the renormalized theory, which has much less degrees of freedom than the original model. In fact, the properties of the renormalized theory can be investigated by simulating the original model written in terms of A µ and measuring quantities written in terms ofÃ µ only. In what follows, we put tildes on all the variables and parameters of the renormalized theory. For instance, corresponding to the cutoffs (2.8), (2.9), we defineκ andL for the renormalized theory bỹ where the symbol · refers to the VEV with respect to the original model for the whole matrices. Let us also define the quantities Let us then consider an effective theory for theN ×N Hermitian matricesÂ µ . (Here and henceforth, we put hats on all the variables and parameters of the effective theory.) We propose  parametersκ,L,Ê andB chosen to beκ,L,B andẼ obtained forÑ =N in figure 5. The results forR 2 (t) obtained in this way are plotted in figure 6. We find that they reproduce the late-time behaviors of the original model very well except the region neart = 0, which is subject to finite "volume" effects anyway. This demonstrates that the effective theory (4.4) indeed captures the late-time behaviors of the original model with much smaller matrix size. Note, in particular, that the symmetric phase, which is not interesting to us, is reduced considerably compared with the original model. In fact, the results forN = 16 do not have a symmetric phase at all. We find it remarkable that the data points agree with those for the original model even in this case.
In application of this method, we actually do not know in advance how to choose the parametersκ,L,Ê andB. This does not spoil the usefulness of the approach at all. In order to show it, we need to answer the question: Can the effective theory (4.4) for arbitrary values ofN ,κ,L,B andÊ be regarded as a renormalized theory forÃ µ in the above sense? By counting the number of parameters in the theory, one finds that the answer is generically "yes". Let us consider the original VDM model with N and κ. Then we define (4.2) and (4.3) for submatrices of sizeÑ , which we denote asκ(Ñ ; N, κ),L(Ñ ; N, κ),B(Ñ ; N, κ), E(Ñ ; N, κ). This specifies a renormalized theory. We try to match it with the effective theory (4.4) by settingÑ =N . In order to matchκ andL withκ(Ñ ; N, κ) andL(Ñ ; N, κ), we can always make a rescaling ofÂ 0 andÂ i as B andÊ should be rescaled accordingly, and we require that they should match (after the Since we have two arbitrary parameters N and κ at our disposal, we can always choose them to satisfy the two conditions in (4.6). (Strictly speaking, since N can take only integer values, the above statement holds as good approximation for sufficiently large N .) Thus we have shown that the effective theory (4.4) for arbitrary values ofN ,κ,L,B andÊ can be regarded as a renormalized theory of the original model for the submatrices corresponding to late times. For this to work, it was necessary to make the rescaling (4.5). This implies that when one makes a plot like the one in figure 3 (Left) for the effective theory (4.4), one should note that the quantity for the x-axis is related to the corresponding quantity in the renormalized theory through where the time-rescaling parameter z is given by . (4.8) Therefore we need to plot our results against the right-hand side of (4.7). Since we do not knowκ(Ñ ; N, κ) in (4.8) a priori, we determine the parameter z in such a way that the results for the model (4.4) scale with the results for the original model at earlier times. In the next section we will show that this is indeed possible, and the method enables us to study the late-time behaviors of the original model in a much more efficient way.

Scaling behaviors in the effective theory
In this section we show how the renormalization group method works by simulating the model (4.4). From now on, we omit the hats on all the variables and the parameters of the model (4.4). We study various values of B and E with N = 32, κ = 4 and L = 1 fixed. From a simulation of the original model with κ = 4 and N = 32, we get B = 7.5 and E = 7.6. The incomplete cancellation E − B ∼ 0.1 is due to numerical artifacts from finiteness of γ C in (A.1). In figure 7 (Left) we show our results for the model (4.4) with various B (including B = 7.5, which corresponds to the original model) for fixed E = 7.6. In figure 7 (Right) we plot the same quantity in physical units. We have introduced the time-rescaling parameter z, which is set to z = 1 for B = 7.5, which corresponds to the original model, and otherwise it is chosen in such a way that the results scale with the results for B = 7.5 at earlier times. Indeed we observe good scaling behavior as anticipated from our arguments in the previous section. We also find that the number of data points in the symmetric phase (t < t c ) decreases as B increases. This means that we can use the matrix degrees of freedom more efficiently for the more interesting broken phase. On the other hand, we find that the lattice spacing ǫ increases slowly as we increase B. We have to make sure that the lattice artifacts are kept under control when we increase B. E = 7.6, which corresponds to the original model, and otherwise it is chosen in such a way that the results scale with the results for E = 7.6 at earlier times. We observe good scaling behavior as anticipated from our arguments in the previous section. We also find that the number of data points in the symmetric phase (t < t c ) decreases as E increases. On the other hand, we find that the lattice spacing ǫ and hence the time-extent ∆ decrease rapidly as we increase E.
The above results suggest a simple strategy for optimizing B and E. First we increase B from the value for the original model until the lattice spacing ǫ becomes a bit too large. Then we can increase E a little in order to make the lattice spacing ǫ sufficiently small. If necessary, we repeat this procedure a few times until the number of data points in the symmetric phase becomes sufficiently small. This way we can increase the time-extent ∆ keeping the lattice spacing sufficiently small with the same matrix size. In figure 9 the region between the two curves in the B-E plane correspond to the case with only one data point in the symmetric phase. Within this region, the lattice spacing ǫ decreases as one increases E.  figure 9. In order to determine the time-rescaling parameter z, we use the results for the original model with N = 64 and κ = 4 as a reference. The lattice spacing ǫ is larger for larger B and smaller E. We find that it becomes too large for B = 35 and E = 7 judging from the deviation from the scaling behavior. The scaling region extends until one reaches B = 29 and E = 8. This gives the maximum time-extent ∆ that one can probe using the renormalization group method with N = 32.
In figure 10 (Right) we plot R 2 (t)/R 2 (t c ) − C against (t − t c )/(zR(t c )), where C is determined by fitting the data in the left panel to an exponential behavior R 2 (t)/R 2 (t c ) = C + (1 − C) exp(bx) with x = (t − t c )/(zR(t c )). We observe a clear straight line behavior    providing strong evidence for the exponential expansion of the early universe in the VDM model.

Summary and discussions
In this paper we have developed a new method for studying the Lorentzian IIB matrix model for a long time period based on the idea of the renormalization group. The method is tested in a simplified model, which captures the behaviors of the supersymmetric model at early times. We were able to confirm the exponential expansion of the space observed in the simplified model with much smaller matrix size using the new method.
On the conceptual side, we consider it interesting that the idea of the renormalization group works in the Lorentzian matrix model. The renormalization group has been applied to matrix models some time ago by refs. [62][63][64][65][66] and more recently by refs. [67,68]. In particular, ref. [68] studies a Yang-Mills two-matrix model as a simplified model of the Euclidean IIB matrix model. In these papers some elements of the matrices were integrated out explicitly to obtain a renormalized theory for matrices of smaller size. A crucial difference from these works is that we have a notion of time in the Lorentzian matrix model. This allows us to consider a renormalized theory for the submatrices representing the degrees of freedom at later times.
The effective theory for the submatrices representing the later time behaviors contains two extra parameters. We have shown how one can optimize them to probe the late-time behaviors most efficiently. The time-rescaling parameter denoted by z has to be determined by requiring that the results should scale at earlier times with the results obtained for the original model. This procedure becomes more complicated if one applies the present method to the original supersymmetric model since the fermionic action (2.3) or (3.1) contains two terms; one of them being proportional to A 0 and the other being proportional to A i . The necessity to rescale A 0 and A i differently as in (4.5) requires us to make the coefficient of the term proportional to A i in the fermionic action of the effective theory, a new unknown parameter. This new parameter can be fixed by probing the scaling behavior. Despite this complication, we consider that the renormalization group method is useful in extracting the late-time behaviors in the Lorentzian IIB matrix model. In particular, from the viewpoint of cosmology, we consider it important to confirm the exponential expansion observed in our preliminary results for the Lorentzian IIB matrix model reported in ref. [61], and to see whether it turns into a power-law expansion at later times as suggested there. We hope to report on these issues in future publications.

A. Details of Monte Carlo simulation
In this appendix we explain how we actually deal with the simplified model (3.3) in Monte Carlo simulation. Generalization to the effective theory (4.4) is straightforward.
First we replace the delta functions and the step function in (3.3) by Gaussian potentials as where the coefficients γ C , γ L , γ κ should be taken large enough to fix each observable to the specified value. In actual simulation we have used γ C ∼ 1 × N 2 and γ L = γ κ ∼ 100 × N 2 . Another important issue concerns the spontaneous breaking of the shift symmetry A 0 → A 0 + α1. For instance, when we try to calculate the expectation value R 2 (t) defined in (2.14), the peak of the quantity measured for each configuration fluctuates considerably. This simply reflects the ambiguity in choosing the origin of the time coordinate, and we should fix it somehow before taking the ensemble average. Here we fix it by introducing a potential 7 where the coefficient is typically taken to be γ sym ∼ 100. We have checked that the results do not alter within error bars for larger values of γ sym .
To summarize, the model we study by Monte Carlo simulation is given by Here p a are real variables, whereas X i are traceless Hermitian matrices. We update all the variables in the model (A.6) as follows. 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 dα a dτ = ∂H ∂p a = p a , dp a dτ = − ∂H ∂α a = − ∂S VDM ∂α a , for some fictitious time τ . This part of the algorithm is called the Molecular Dynamics. In solving the Hamilton equations (A.7) numerically, we discretize them using the socalled 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), based on the idea of the Metropolis algorithm to satisfy the detailed balance. The crucial 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 is non-zero due to the discretization, but it is a small quantity of order (∆τ ) 2 . Therefore, one can move around efficiently in the configuration space.
Since the auxiliary variables p a and (X i ) ab appear only as the Gaussian terms in (A.6), 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. Thus the HMC algorithm as applied to our system can be described as follows.
2. Evolve the fields p a (τ ), X i (τ ), α a (τ ) and A i (τ ) for fictitious time τ f according to the discretized Molecular Dynamics.
The HMC algorithm involves two parameters ∆τ and τ f , which can be optimized. (See, for instance, appendix B of ref. [38] for more details.) For fixed τ f , we have to choose ∆τ so that ∆τ × (acceptance rate) is maximized. Typically this is achieved for acceptance rate of 50∼60%. Then τ f can be optimized to minimize the autocorrelation time in units of one step in the Molecular Dynamics.