Phases of a matrix model with non-pairwise index contractions

Recently a matrix model with non-pairwise index contractions has been studied in the context of the canonical tensor model, a tensor model for quantum gravity in the canonical formalism. This matrix model also appears in the same form with different ranges of parameters and variables, when the replica trick is applied to the spherical $p$-spin model ($p=3$) in spin glass theory. Previous studies of this matrix model suggested the presence of a continuous phase transition around $R\sim N^2/2$, where $N$ and $R$ designate its matrix size $N\times R$. This relation between $N$ and $R$ intriguingly agrees with a consistency condition of the tensor model in the leading order of $N$, suggesting that the tensor model is located near or on the continuous phase transition point and therefore its continuum limit is automatically taken in the $N\rightarrow \infty$ limit. In the previous work, however, the evidence for the phase transition was not satisfactory due to the slowdown of the Monte Carlo simulations. In this work, we provide a new setup for Monte Carlo simulations by integrating out the radial direction of the matrix. This new strategy considerably improves the efficiency, and allows us to clearly show the existence of the phase transition. We also present various characteristics of the phases, such as dynamically generated dimensions of configurations, cascade symmetry breaking, and a parameter zero limit, to discuss some implications to the canonical tensor model.


Introduction
Quantization of gravity is one of the most challenging fundamental problems in physics, and various approaches to this problem have been proposed so far. These include sophisticated applications of the renormalization group procedure to general relativity [1] as well as approaches that use discretization of spacetime in the definition of theory, for instance the approaches in [2,3,4,5], matrix models [6,7,8,9,10], and tensor models [11,12,13,14]. One of the goals of these discretized approaches is to show the emergence of macroscopic spacetime as a continuous manifold, with general relativity emerging as the effective description of dynamics. This is still a challenging goal for any of these approaches.
In this paper, we study the dynamics of a matrix model which contains non-pairwise index contractions [15,16]. This matrix model has only the O(N ) × S R symmetry for the index spaces of the matrix variable φ i a (a = 1, 2, . . . , N, i = 1, 2, . . . , R), where S R denotes the symmetric group and O(N ) the orthogonal group. The symmetry is not enough to diagonalize an arbitrary matrix, and therefore this matrix model is not solvable by the methods usually employed to solve matrix models [6,7,8,9,10] or rectangular matrix models [17,18,19]. Our matrix model can also be regarded as a vector model of R vector variables, but our setup is different from the exactly solved ones in [20,21].
The background motivation for our matrix model comes from the fact that this model has an intimate connection [15,16] to an exact wave function [22] of a tensor model in the Hamilton formalism, which we call the canonical tensor model [23,24]. Previously it has been found that this wave function peaks around Lie-group symmetric configurations of the tensor-variable of the model [25,26]. This is encouraging towards potential emergence of a spacetime as mentioned above, because Lie-group symmetries, such as Lorentz, deSitter, gauge, and so on, are ubiquitous in the universe. However, it is still difficult to show whether the peaks contain configurations which can be interpreted as some sort of spacetime, for instance, in the manner described for a classical treatment in [27]. Understanding the properties of the dynamics of our matrix model will potentially provide useful insights about the relation between the wave function of the tensor model and spacetime emergence.
It is an intriguing coincidence that a matrix model with the same form has previously appeared in the context of spin glasses. It is obtained, when the replica trick is applied to the spherical p-spin model (p = 3) [28,29] for spin glasses, where R designates the replica number. However, the physics of the spin glass and that of our model will be largely different, because the parameter and variable regions of interests are different from each other. In the spin glass case, the replica number R is taken to the limit R → 0 as part of its process, while our interest is rather in the limit R ∼ N 2 /2 → ∞, the reason of which comes from the consistency of the tensor model, as explained more in Section 7. In addition, the coupling parameter (called λ in later sections) of the models has opposite signs, and the spin glass case has a spherical constraint, φ i a φ i a = 1, for each i. Considering these differences, it seems necessary to analyze our model independently from the spin glass case.
In the previous paper [16], Monte Carlo simulations of the model were performed with the usual Metropolis update method. This has revealed various interesting characteristics of the model. However, there was an issue which affects the reliability of the Monte Carlo simulations: For some values of the parameters important to study its properties, the iterative updates in the radial direction of the matrix variable were too slow to reach thermodynamic equilibriums in a reasonable amount of time. For instance, it could not be determined with confidence whether the transition is a phase transition or just a crossover, since the parameters could not be tuned to make the transition more evident. The major improvement of the present paper is that we integrate out the troublesome radial direction before doing the numerical calculation, obtaining a model essentially defined on a compact manifold (the hypersphere). In addition, the more efficient Hamiltonian Monte Carlo method is employed instead of the more straightforward Metropolis algorithm. This replacement of the model drastically improves the efficiency of the simulations, and we have successfully obtained much more evident results than the previous ones.
We summarize below the properties of the transition and the phases derived from the numerical results: • The transition becomes sharper as N is taken larger. This implies the transition is a phase transition in the thermodynamic limit. We have not observed any discrete behavior of observables around the transition point, implying that the transition is continuous.
• The value of R at the transition point, which we call the critical value R c , is a little smaller than (N +1)(N +2)/2, that was previously obtained by the perturbative analytic computations in [15,16]. The critical value R c is better approximated by R c ∼ (N + 1)(N + 2)/2 − N + 2 in our numerical results, where the parameters of the model are taken in the range k/λ O(10 −10 ) and N 12.
• It has been shown that the Monte Carlo results and the results of the perturbative analytical computations in [15,16] do not agree with each other near R ∼ R c . The ratios between them on the peaks increase with the decrease of k/λ, and increase or converge 1 to some k/λ-dependent limiting values with the increase of N . Away from R ∼ R c , the two approach each other, quickly for R > R c and gradually for R < R c .
• It has been reported in [16] that the dimensions of the configurations change under the change of R in the vicinity of the transition point R c . This behavior is more precisely investigated in this paper. For small k/λ, the dimensions take the smallest values at the transition point, and take larger values as R is taken further away from R c .
• Though we worked with the improved setup explained above, we still encountered rapid slowdown of iterative updates of Monte Carlo simulations in the parameter region R R c and k/λ O(10 −8 ). However, the slowdown seemed to be smoothly improved by taking smaller step sizes and performing longer simulations. This implies that there is no transition associated to the slowdown. Thus, for R R c , the model behaves like a fluid with a viscosity which continuously grows as k/λ decreases. 1 We could not conclude which one is the right behavior from the simulation datas, as we will see later.
• For R R c , it has been observed that the dynamics of the model converges in the k/λ → +0 limit, in which expectation values of observables and the free energy converge to finite values.
• For R R c , it has been observed that the free energy diverges in the limit k/λ → +0.
• SO(N ) symmetry breaking occurs at R R c due to large φ i a . This occurs in a cascade manner as R increases: the breaking of SO(N ) occurs first, then SO(N − 1), . . . , and finally SO(2) breaks down.
We also discuss some implications of the numerical results to the tensor model. The most important is the coincidence between the location of the transition point and a consistency condition of the tensor model in the leading order of N . Combining this with the result that the phase transition is continuous, this suggests the possibility that a continuum theory can be associated to the tensor model. Moreover, the fact that the transition point is where the dimensions of the configurations quickly decrease towards low values suggests the possibility of emergent spacetimes with sensible dimensions in the tensor model. This paper is organized as follows. In Section 2, we explain the matrix model and derive the new setup for the numerical simulations, which is obtained by integrating out the radial direction of the matrix variable. In Section 3, observables are introduced. There are roughly two classes of observables, one directly related to the matrix model, and the other directly related to our setup. The two classes are connected by a formula. In Section 4, we derive some formulas which compute the expectation values of some observables by using the analytic results obtained previously in [15,16]. In Section 5, we comment on our actual Hamiltonian Monte Carlo method for the angular variables in our setup. In Section 6, we summarize our results of the Monte Carlo simulations. In Section 6.1, we present several pieces of evidence of the phase transition. In Section 6.2, we compare the results of the numerical simulations and the analytic perturbative computations, and show that there are differences in the vicinity of the transition point, which grow or converge to some k/λ-dependent values as N increases. In Section 6.3, the k/λ → +0 limit is discussed. Its behaviour severely differs in the two phases. In Section 6.4, the geometry of dominant configurations is discussed. In particular, the dimensions take minimum values at the transition point. In Section 6.5, symmetry breaking in a cascade manner for R ≥ R c is shown. In Section 6.6, the slowdown of iterative updates in our simulations is discussed. This appears to occur quickly as k/λ becomes smaller at k/λ O(10 −8 ) in our simulations, but a quantitative investigation shows that this is a smooth change, implying that there is no transition to another phase with slow dynamics. In Section 7, the implications of the numerical results to the tensor model are discussed. In Section 7.1, the coincidence of the transition point with a consistency condition of the tensor model is discussed. In Section 7.2, the behavior of the dimensions is explained from the symmetrypeak relation argued in [25,26]. In Section 7.3, the normalizability of the wave function of the tensor model is discussed. The last section is devoted to a summary and future prospects.

The matrix model and the setup for simulations
The matrix model we consider in this paper is defined by the partition function, where φ i a (a = 1, 2, . . . , N, i = 1, 2, . . . , R) denote the matrix variable, the integration is over the whole N R-dimensional real space, and the coupling parameters, k and λ, are assumed to be positive real for the convergence of the integral as will be explained in more detail below.
where the repeated lower indices are assumed to be summed over. Throughout this paper, repeated lower indices always appear pairwise, and we assume the common convention they are summed over, unless otherwise stated. On the other hand, the upper indices are triply or sixfold contracted in (1), and summation over them will always be written explicitly.
The matrix model (1) has the O(N ) × S R symmetry, where O(N ) denotes the orthogonal group transformation in the N -dimensional vector space of the lower index, and S R denotes the permutation symmetry for the upper index values {1, 2, . . . , R}. The O(N )×S R symmetry is generally not enough to diagonalize the matrix φ i a , and therefore the model cannot exactly be solved by the well-known methods often applied to the usual matrix models [6,7,8,9,10] or the rectangular matrix models [17,18,19]. Because of the O(N ) symmetry, the model can also be regarded as a vector model [20,21] with the multiplicity of vectors labeled by the upper index. In fact, our model can be solved in the N → ∞ limit with finite R [15], as in the vector models [20,21] and in the spherical p-spin model [28,29]. However, this solution is not so useful, because our major interest is the vicinity of the phase transition point with R ∼ N 2 /2, as will be explained later.
The first term of the exponent of the matrix model (1) is positive semi-definite, since The equality on the rightmost is actually satisfied by various configurations, including straightforward ones like φ 1 a = −φ 2 a , . . .. Moreover, when R is larger than a certain value, there will be a continuous infinite number of solutions. 2 Therefore, if k = 0, it is not obvious whether the integral (1) is convergent or not. On the other hand, if k > 0, one can immediately see 2 A simple counting of degrees of freedom implies that the dimension of the solution space of R i=1 φ i a φ i b φ i c = 0 will be given by N R − N (N + 1)(N + 2)/6, where the former counts the degrees of freedom of φ i a and the latter the number of independent conditions. Therefore, in general for R > (N + 1)(N + 2)/6, the solutions to the equality will exist continuously.
Thus k > 0 assures the convergence of the integral (1) for general cases, while it will be shown later that the k/λ → +0 limit can be taken if R < R c .
As explained above, the second term in the exponent of (1) acts as a regularization of the integral. A term with the same role existed in the previous studies of the model [15,16], but had a different, namely quadratic, form, The main reason for this choice of quadratic form was that then the action (the exponent) had the standard form used in perturbative computations. It is however not necessary to take a quadratic term as a regularization term for the perturbative computations 3 , as we will review in Section 4. The present choice R i=1 U ii (φ), which has the same order as the first term, is more convenient in the current analysis, because then the radial direction of φ i a can be integrated out in a straightforward way, as we will perform below. Since the radial direction was the main source of the difficulties in the previous simulations [16], the present choice will ease the deadlock of the simulations. Now let us divide φ i a into the radial and the angular coordinates, φ i a = rφ i a , where r denotes the radial coordinate, andφ i a denote the angular coordinates. Putting this reparameterization into (1) and integrating over r, one obtains where Γ(·) is the gamma function, and S N R−1 denotes the unit N R − 1 dimensional sphere. We will use (4) as the weight of our simulations, where the variables are only the angular ones. Our implementation of the Hamiltonian Monte Carlo method for this system will briefly be explained in Section 5.
Finally, let us comment about the relation between the matrix model (1) and the spherical p-spin model for spin glasses, leaving the relation to the canonical tensor model for Section 7. The partition function of the spherical p-spin model for p = 3 is given by with a random real coupling P abc to simulate a spin glass system. Considering R replicas of the same system in the replica trick and simulating the random coupling by a Gaussian distribution e −αP abc P abc with positive α, one obtains where N is an overall coefficient. The righthand side has a similar form as (1), but there are two major differences. There are the restrictions, φ i a φ i a = 1 for each i, which assure the finiteness of the integration, taking the role of the second term in the exponent of (1). The other difference is that the coefficient of the exponent has the inverse sign compared to (1). This physically means that the dominant configurations will be largely different between the matrix model (1) and that of the spherical p-spin model. In addition, the R → 0 limit is finally taken as part of the replica trick, while this is not necessary in the matrix model (1) itself. We are rather interested in the dependence on R of the system, especially in the regime R ∼ N 2 /2, as we will see later. In particular, the last relation requires R → ∞ in the thermodynamic limit N → ∞, which is opposite to the spin glass case.

Expectation values of observables
For convenience, let us first slightly generalize the definition of the partition function (1) of the matrix model to We assume the symmetric matrix coupling Λ is taken so that the integral is convergent. This includes the original case (1) with Λ = Λ λ,k , where with positive λ, k.

Let us introduce
where U ij (φ) := (φ i aφ j a ) 3 , and O(φ) is an arbitrary observable expressed as a function ofφ i a . As derived in Section 3, by integrating over r, the partition function (7) can be expressed by where ∆ N R := 1 6 N R. From (10), the expectation value of an observable O(φ) is given by These are the observables which are directly obtained in our Monte Carlo simulations for the angular variables.
There is another kind of observable that is expressed as a function of φ i a . The difference is just the normalization of φ i a . Let us introduce a weight [·] which counts the multiplicity of φ i a contained in an observable that is assumed to be a homogeneous function of φ i a . For example, the weight of φ i a φ i a is given by [φ i a φ i a ] = 2. Let us consider an observable O(φ) with weight w. This can be rewritten as O(φ) = O(φ) r w by the reparameterization φ i a = rφ i a with the radial and angular variables. Then the expectation value is given by where we have used the obvious property, z N, This formula relates the expectation values of the observables of the matrix model (1) expressed by φ i a with those in our Monte Carlo simulations for the angular variables. This formula is used for deriving some results in Section 6.

Analytic computations by a perturbative method
In the previous papers [15,16], the authors introduced a function defined by where V S N R−1 designates the volume of the unit sphere, S N R−1 dφ, for the normalization f N,R (Λ = 0) = 1. This function is related to the partition function (7) of the matrix model by A merit of introducing the function (13) is that it can obviously be defined for arbitrary complex values of Λ ij since the integral region in (13) is compact, meaning that it is an entire function of Λ ij [15]. Therefore, if the whole perturbative series expansion in Λ ij of this function is obtained, this is convergent for any Λ ij = ∞, and hence determines the function completely in the whole complex region of Λ ij . This means that, in principle, the dynamics of the matrix model (1) can be determined as precisely as one can by improving the perturbative series expansion of f N,R (Λ). This is in contrast with the partition function Z N,R (λ, k), which is singular at λ = 0 or k = 0 and merely an asymptotic perturbative series expansion of it in λ or k can be obtained.
The perturbative computations of f N,R (Λ) using Feynman diagrams have been performed in [15,16]. This can be done by mapping the integrals S N R−1 dφφ i 1 a 1φ i 2 a 2 · · ·φ in an , which appear in the expansion of the integrand in (13), to the standard computations using Wick contractions. The final result derived in the leading order is given by where the product is over the eigenvalues of the matrix Λ ij with degeneracies taken into account, and with In the case with Λ λ,k ij = λ + k δ ij , the eigenvalues are k + λR for the eigenvector (1, 1, . . . , 1) and k for all the other vectors transverse to that. Therefore To obtain formulas for expectation values of observables, let us introduce From (13), one finds Therefore, it has a relation with (9) as Then, by comparing with the results in Section 3, the correlation functions of U ij (φ) := (φ i aφ j a ) 3 for the angular variables can be expressed as Therefore, by combining with (15) (or (18)) and (19) and numerically integrating over t, one can compute the correlation functions of U ij (φ) in the leading order of the analytic perturbative computation.
Let us next consider the correlation functions of U ij (φ) for the variable φ i a . We can use the formula (12), where each of U ij (φ) has weight w = 6. We obtain This also gives the correlation functions of U ij (φ) in the leading order from the analytic perturbative method.
Later we consider an observable given by 3 . In our actual case with Λ = Λ λ,k given in (8), the derivatives in (22) and (23) for the observable can be performed by ∂ ∂k . Therefore the correlation functions are given by This formula is used when we compare the numerical results with the analytical ones in Section 6.2.

Hamiltonian Monte Carlo method for angular variables
In this paper, we use Hamiltonian Monte Carlo method [30] for the numerical simulations. This method upgrades the configuration space of some integral to a phase space by introducing conjugate variables, and creates a Hamilton system and (locally) solves the equations of motion in order to find new, more remote, candidates for the Metropolis update. This process is called leapfrog, which consists of a sequence of discrete jumps from one phase space location to another. While it is enough for presenting update candidates to approximately solve classical equation of motion, the time reversal symmetry and the conservation of phase space volume must be exactly satisfied under the discrete jumps for correct sampling of configurations. For a flat configuration space, these conditions are easily satisfied by alternately sequencing the following two processes: where (q i , p i ) designate phase space variables indexed by i, is the size of one jump, and V (q) is the Gibbs potential for a weight exp(−V (q)). Observe that (i) is a free motion in a flat space, and only (ii) takes effects from V (q). Each of the two jumps obviously satisfies the conservation of the phase space volume, det |∂(q i + δq i ), ∂(p i + δp i )/∂q j , ∂p j | = 1, due to the fact that q i and p i do not jump simultaneously. The time reversal symmetry is also satisfied, When the configuration space q i is constrained to a non-flat sub-manifold embedded in a flat space, a free motion corresponding to (i) is generally a simultaneous jump of p i and q i , since the tangent space of the sub-manifold containing p i changes along q i . In such a case, finding an appropriate jump corresponding to (i) satisfying the two necessary conditions above is generally a difficult problem. An obvious solution to an appropriate jump is to exactly solve the classical equation of the free (geodesic) motion on the sub-manifold [31]. This is possible when a sub-manifold is simple enough to allow us to obtain such exact solutions. In our case, the embedded manifold is a unit hypersphere, which gives the constraints, i q 2 i = 1 and i q i p i = 0, and the jump describing the exact free (geodesic) motion on the sphere is given by where |p| = i p 2 i and θ = |p|. The second jump (ii) does not contain a jump in q i , therefore there are no difficult issues, and it can just be replaced by where the additional term takes into account the constraint i q i p i = 0.
In our present case (4), the coordinatesφ i a are constrained on a unit sphere R i=1φ i aφ i a = 1, and we employ these jumps (i') and (ii'). The potential energy can be read from (4) as with Λ = Λ λ,k .

Results of Monte Carlo simulations
In this section, we summarize the results of our Hamiltonian Monte Carlo simulations from several view points. Since the overall factor of the exponent of (1) can be absorbed in the rescaling of φ i a , we set λ = 1 in all the simulations, leaving N , R, and k as variable parameters. Errors were estimated by the Jackknife method described for example in [32]. We took the leapfrog numbers to be about 1000-10000, depending on the hardness of the simulations explained in Section 6.6, and the step sizes were tuned so that the acceptance rates were about 80-99 percent, which were a little higher than the commonly taken ones because of the reason explained in Section 6.6. Parallel tempering [33] was also used in some of the computations to take some datas which systematically study k-dependencies. However, as will be explained more in Section 6.6, parallel tempering did not seem to essentially affect the expectation values computed.
In the following subsections, we show the results of the simulations of the expectation values of various observables depending on the purposes. The observables are taken to be invariant under the O(N ) × S R symmetry.

Phase transition point
There are various observables which can be used to study the location of the phase transition. We will present one example forφ i a and another for φ i a . The observable we first consider is An important reason for considering this observable is that this has the natural normalization factor N . This factor is determined by the uncorrelated case, in which each ofφ i a is regarded as an equally independent variable. More precisely, the uncorrelated case corresponds to φ i aφ j b uncorrelated ∼ δ ab δ ij /(RN ) up to sub-leading corrections in N and R by taking into account the constraint, R where we have ignored sub-leading corrections in N and R. Figure 1 shows the results of the Monte Carlo simulations for O 1 . The normalization factor R c for the horizontal axes is chosen as R c = (N + 1)(N + 2)/2 − N + 2. The perturbative computations in the leading order predict the transition point to be at R c = (N + 1)(N + 2)/2 [15,16]. However, for the datas shown in the left figure for k = 10 −8 , it is better to take R c = (N + 1)(N + 2)/2 − N + 2 to locate all the peaks near R/R c = 1. The values of O 1 approach 1 as R takes more distant values from R c , implying that the correlations become more independent there. On the other hand, the values of O 1 at the peaks become larger for larger N . This can be checked more clearly in Figure 2. This means that the correlation becomes larger at the transition point for larger N , which is a typical signature of a continuous phase transition.
The right picture of Fig. 1 shows the dependence of O 1 on k for N = 10. The dependence on k seems little for R R c , as we will discuss more of this aspect in Section 6.3. On the other hand, at R R c , O 1 seems to become larger as k becomes smaller, slightly shifting the locations of the peaks to the right. This implies that the correlations become larger for smaller k and the critical value R c depends not only on N but also on k as well. The last statement implies that what we have taken as R c above cannot be considered to be a correct expression valid for the general values of the parameters, but can at most be considered to be an approximate expression valid for our parameter range N 12 and 10 −10 k 10 −8 .
Let us next turn to the observable, From the formula (12), by setting the weight w = 2 and noting O 2 (φ) = 1 identically, we obtain with Λ = Λ λ,k . The results of the simulations for r 2 are plotted in Figure 3. The figures clearly show that the two phases are characterized by r 2 ∼ 0 for R < R c and r 2 > 0 for R > R c , respectively. The transition becomes sharper as N becomes larger or k becomes smaller. r 2 changes continuously at R ∼ R c , supporting the claim that the transition is continuous.

Comparison with the perturbative computation
In this section, we compare the results of the simulations with the analytic perturbative computation in the leading order, which was reviewed in Section 4. In particular, we see that the analytic computation does not explain the peaks of the correlations ofφ i a , which was shown in Section 6.1. We find clear deviations between them around the phase transition point R ∼ R c , while they converge as R takes distant values from R c .
To see this we consider the observables, , whose formulas of the analytic computation are given in (24). The explicit values are obtained by performing the numerical integration of (19) with (18) contained in (24) for M = 1. On the other hand, we compare these with U d (φ) and from the simulations, where we have used (12) for w = 6. Figure 4 shows the comparison between the Monte Carlo results and the analytic computations. There exist systematic deviations in the vicinity of R = R c , as was previously reported in [16]. For R > R c , they quickly converge as R leaves R c . For R < R c , they slowly converge as R becomes smaller.
One can see similar deviations for the U d (φ) . Figure 5 plots the ratio U d (φ) / U d (φ) pert between the Monte Carlo results and the perturbative analytic computations. Indeed the ratio deviates from 1 in the vicinity of the transition point. The deviations at the peaks become larger as k becomes smaller. On the other hand, as shown in Figure 6, it seems that the deviations increase with N for k < 10 −8 , but this is not clear for k ≥ 10 −8 . We cannot rule out the possibility that they actually converge in the large N limit to some values which increase with the decrease of k.
From the comparisons above, we conclude that the perturbative analytic computation in the leading order does not correctly reproduce the behavior of the matrix model in the vicinity of the transition point. As was previously performed in [16], the situation does not   essentially change, even if we take into account the next leading order corrections to the analytic computation.

k/λ → +0 limit
In this subsection, we focus on the k/λ → +0 limit of the matrix model (1). There are a few reasons to study this. One is the characterization of the phases separated at R = R c . We find different limits for each phase at R > R c and R < R c . Another is its relevance to the tensor model. The behavior determines whether the wave function is normalizable or not. This will be discussed in Section 7.
Firstly, let us show that the limit k/λ → +0 converges in the phase R < R c . This can be seen by looking at the behavior of expectation values of observables. Figure 7 shows the result of the simulation about the behavior of U d (φ) in (33) against k for a case with R < R c . As can be seen in the figure, the expectation value approaches a constant value in the k → +0 limit. In fact, similar convergence can be observed also for other observables in other cases with R < R c .
Let us discuss the consequence of this behavior to the free energy defined by F N,R (λ, k) := − log Z N,R (λ, k). By taking the derivative of (1) with respect to k, we obtain Therefore, as a function of k, F N,R (λ, k) can be determined by studying the k-dependence of U d (φ) and performing integration: In particular, lim k→+0 U d (φ) will determine lim k→+0 F N,R (λ, k). Below let us discuss the behavior of the free energy in k/λ → +0. By performing the rescaling of the variable φ i a → λ −1/6 φ i a in (1), one obtains In the region R < R c , there is a finite limit of lim k→+0 U d (φ) as shown above. Considering (35) and (36), the behavior of the free energy is obtained as where U 0 d := lim k→+0 U d (φ) λ=1 , and p N,R (k/λ) is smaller than k/λ in order and has a finite limit p N,R (+0). We comment that this finiteness was proven analytically for R = 2 and any N previously in [16] 4 . This finiteness of F N,R (λ, k) in the k → +0 limit is non-trivial, as discussed in Section 2.
On the other hand, for R > R c , the simulations show that U d (φ) diverges in the k → +0 limit. An interesting matter is that, instead, k U d (φ) converges in the k → +0 limit, as can be seen from Figure 8. This implies that, from (35), F N,R (λ, k) logarithmically diverges in the limit k/λ → +0.
Let us discuss this divergence of the free energy in more detail. As we have seen in Figure 8, if we take k small enough, k U d (φ) can be regarded as its limiting value, lim k→+0 k U d (φ) . By assuming this for the N = 10 data in the large R region in the left figure of Figure 4, and fitting a linear function of R for the data in the region R > 1.4 · R c , one obtains, This curiously agrees very well with what can be obtained by putting N = 10 to a hypothetical expression for the righthand side, where #P := N (N + 1)(N + 2)/6 is the number of independent components of a symmetric three-index tensor P abc . We have performed similar analyses for N = 5, 7 cases and have found good matches with the hypothesis (39). Assuming the hypothesis and reminding the form (36), we obtain whereŨ 0 d := lim k→+0 k U d (φ) ,p N,R (k/λ) is smaller than log(k/λ) in order, and with δŨ 0 d sub-leading in large R. Note thatŨ 0 d ≥ 0 due to U d (φ) > 0, and therefore δŨ 0 d takes positive values in the range R c < R < (N + 1)(N + 2)/2.
Here it is a non-trivial question whetherp N,R (k/λ) has a finite limitp N,R (+0). For example, a slow correction of order ∼ 1/ log(k) to k U d (φ) for k ∼ +0 leads to a double logarithmic divergence ofp N,R (+0). However, as shown in Figure 8, the data points of k U d (φ) can be fitted very well with a correction of order √ k, and there is no good motivation for introducing such slow corrections. Therefore it would be reasonable to assumep N,R (+0) to exist as a finite value. We also comment that the hypothesis (39) is nothing but what can be obtained from the perturbative computation in the leading order [15], as the coincidence in the left figure of Figure 4 shows. Therefore δŨ 0 d is a correction beyond the leading order perturbative computation.
The difference between the behavior of the free energy (37) and (40) characterizes the two phases separated by R = R c . These formulas will be used in Section 7.3, where we will discuss the normalizability of the wave function of the tensor model.

Geometric properties
In Section 6.2, we have found the deviation between the results of the simulations and the perturbative analytic results in the vicinity of the phase transition point. This suggests that some non-perturbative configurations are important in the vicinity of the phase transition point. In this subsection, to discuss the characteristics of the configurations around the phase transition point, we study the distributions of the vectors φ i a (i = 1, 2, . . . , R) in the Ndimensional vector space associated to the lower index. These vectors define a point cloud with R points, where φ i a for each i determines the location of each point in the N -dimensional vector space. We study the dimensions of such point clouds. It turns out that the dimensions depend on the parameters of the matrix model. In particular, the dimensions take the smallest values at the transition point as functions of R.
To study the dimension of a point cloud we use angle distributions among the vectors. The angle between two vectors, say φ i a and φ j a (i = j), in the N -dimensional vector space is given Assuming that the vectors approximately form a rotationally symmetric d-dimensional point cloud, the distribution of the angles should be approximately given by where θ designates the angle, and N is a normalization factor. This formula can easily be obtained by radially projecting points to the unit sphere S d−1 , and computing infinitesimal areas associated with given mutual angles. The dimensions can be computed by fitting the formula (43) to the angle distributions obtained from the datas.
In Figure 9, we show two examples of the fitting. As shown in the figures, the fitting is generally quite good for high dimensions but not so much for lower dimensions. The reason behind this is that the point clouds cannot be characterised as a single dimensional object but are a mixture of objects with different dimensions, as we discuss in Section 7.2. Yet, to characterize the configurations in terms of dimensions, we perform the fitting restricted to a portion around the center, namely θ ∼ π/2, because there exist dominant numbers of cases in this region. In this sense, the dimension is merely a qualitative characterization, but it still gives a fairly interesting observable: The dimension takes the lowest value at the phase transition point as a function of R. For instance, this can be observed for N = 10, k = 10 −8 in Figure 10.
It is instructive to directly see a point cloud itself. A point cloud exists in an N -dimensional space, but if its dynamical dimension is lower than three, one can project it into a threedimensional space by extracting the main three extending directions through principal component analysis (PCA). Figure 11 shows a collection of a number of such projected point  Left: The collection of the point clouds obtained from the simulation with N = 10, R = 57, and k = 10 −8 . The point cloud from each data of φ i a is projected into the three-dimensional space and the collection through all the datas are plotted. For the projection, PCA is used to take three major directions out of N dimensions. Right: The corresponding density plot. The shape is like a squashed rugby ball, which may be regarded as an object with a dimension between 2 and 3. clouds, which have been sampled from the simulation with N = 10, R = 57, k = 10 −8 . According to Figure 10, the point cloud has a dimension nearly two in this case, and we indeed find an approximately two-dimensional object which has the shape of a squashed rugby ball as shown in the right figure of Figure 11.
Finally, let us discuss the k-dependence of the dimension. The general behavior is that the dimensions decrease with the decrease of k and converge to limiting values, as is shown in Figure 12.

Symmetry breaking
As shown in Section 6.1, the phase at R > R c is characterized by large values of r 2 . Since a non-vanishing value of φ i a breaks the O(N ) symmetry associated to the lower index vector space, the phase at R > R c will be characterized by symmetry breaking. In this subsection, we will study this aspect.
Let us consider one of the generators T ab of SO(N ). The size of the breaking of T ab by a vector φ i a will be characterized by the size of the vector T ab φ i b . By considering its square and summing over all the vectors, the breaking by a configuration can be characterized by Thus the natural quantity to study is where T where e φ a (a = 1, 2, . . . , N ) are the eigenvalues of the matrix m. Figure 13 gives two examples of the eigenvalues eg(M ). In the figures, the eigenvalues are plotted in ascending order along the horizontal direction. In the case of the left figure, one can find an interesting stair-like pattern of the eigenvalues. This pattern means that the original SO(N ) symmetry is hierarchically broken to SO(N −1), SO(N −2), . . .. In fact, the horizontal locations of the steps agree with the numbers of the generators of these symmetries. On the other hand, in the case of the right figure, all the symmetries are broken with no obvious hierarchal structure.
Since the pattern above generally fluctuates over the samples of φ i a in a simulation, we consider an average, eg(M ) . The precise definition of this quantity is as follows: we run a simulation with a certain choice of parameters; for each sample of φ i a in a simulation, we compute eigenvalues eg(M ) and order them in ascending order; then we take mean values of each entry over all the datas of the simulation. Figure 14 shows eg(M ) computed from the simulations respectively for R = 60 (left) and for R = 80 (right) with N = 10, k = 10 −8 . Figure 15 shows the dependence of eg(M ) over the change of R for N = 10, k = 10 −8 . The eigenvalues start to increase from R ∼ 0.9R c with the increase of R. The symmetry breaking occurs one by one: first SO(10), then SO(9), and so on, until finally all the symmetries are broken at R ∼ 1.3R c . In the figure, one can find some gaps between the eigenvalues  in the vicinity of R ∼ R c . They correspond to the differences of the step heights, which for example exist in left figure of Figure 14. As R becomes larger, the gaps gradually disappear, approaching the situation in the right figure of Figure 14.
The above symmetry breaking in a cascade manner is consistent with the results in the previous subsections. When R < R c , since φ 2 is small, there is no symmetry breaking. As R increases from R ∼ R c , the vectors φ i a (i = 1, 2, . . . , R) start to take larger values and fill a subspace, the dimension of which increases with the increase of R. Since the subspace breaks part of the SO(N ) symmetry, depending on its dimensions, more symmetries are broken with the increase of R.
Let us comment about the fate of the discrete symmetry φ i a → −φ i a ∀i. For N = odd, this corresponds to the Z 2 subgroup of the O(N ) symmetry. The quantity 6 , R i=1 φ i a /R, is not invariant under the discrete symmetry, and therefore the expectation value of its square, R i,j=1 φ i a φ j a /R 2 , will provide a good quantity to measure its breaking. It has turned out that the expectation values computed from the simulation datas stay small in the order O(1) over the range, and we have not observed any signatures of its breaking.

Slowdown of Monte Carlo updates
In the previous work [16], we encountered a rather serious difficulty of the Monte Carlo simulation: For R R c and small k, the step sizes of the simulations had to be tuned very small for reasonable acceptance rates of Metropolis updates, but then the updates of configurations were too slow for the system to reach thermodynamic equilibriums within our runtimes. Therefore, in this paper, we have improved the strategy: Integrating out the radial direction of the model and using the so-called Hamiltonian Monte Carlo method for simulations. Indeed the new strategy drastically improves the efficiency of the simulations, but we still encounter the slowdown for smaller k, which is however several orders of magnitude smaller than that in the previous work. This implies that this slowdown is an intrinsic property of the model, which is independent from methods of simulations, and would even suggest a possibility of the presence of a transition to a new phase characterized by slow dynamics. However, in this subsection, we will show that the last possibility is unlikely, and the system in the phase at R > R c is rather like a fluid with a viscosity which continuously grows for smaller k.
The speed of updates can be quantified by the mean value of distances between neighboring configurations in a sequence of updates,φ i a (1),φ i a (2), . . . ,φ i a (M + 1): where In the ideal maximum situation that each entry of the sequence is independent from the others, (δφ) 2 = 2 because of the normalization |φ(m)| 2 = 1. Figure 16 shows the dependence of (δφ) 2 against the value of k for R = 45 (left) and R = 80 (right), respectively, with N = 10. In the simulations, the step sizes, namely the value  in Section 5, are properly chosen for reasonable acceptance rates 7 for each k, while the other parameters of simulations, such as leapfrog numbers, are fixed. The R = 45 case keeps the ideal values around log 10 2 ∼ 0.3 throughout the shown range of k. On the other hand, the R = 80 case has a rapid decrease of the speed with the decrease of k at k 10 −6 .
The speed of updates defined above is dependent on the parameters of the simulation such as step size, leapfrog number, and even the frequency at which the data is saved, and is therefore not a quantity intrinsic to the model. For instance, the starting point of decreasing, k ∼ 10 −6 , has no physical meaning, since this can easily be changed by taking different simulation parameters. However, in the data above, the only parameter which is varied is the step size among different values of k, and it is therefore meaningful to compare the data for different values of k by rescaling (δφ) 2 / 2 to cancel the obvious dependence on . The left figure of Figure 17 plots the values of taken for the simulation of N = 10, R = 80, and the right is for the corrected values, (δφ) 2 / 2 . In the right figure, leaving aside the irrelevant ideal region k 10 −6 , one can see that the values are almost flat in the region k 10 −6 with no essential change. This implies that the system is basically similar up to the obvious rescaling among different values of k. We also used parallel tempering [33] in addition to the Hamiltonian Monte Carlo method for taking some datas which systematically study k-dependencies. The exchanges of configurations were performed among different values of k, typically taken k = 10 −n (n = 2, 3, . . . , 11), with common values of the other parameters. In the region R > R c , as R increases from R c , the exchange rate quickly reduces for the above choices of k's. Therefore, parallel tempering does not seem effective to solve the slow update problem, which exists at R R c for small k. On the other hand, the exchange rate is high for small k at R R c , which can easily be understood by the presence of the well-defined k/λ → +0 limit at R < R c , as discussed in Section 6.3: the sets of configurations are similar among different values of k, when k is small enough. However, we did not observe any major differences between the datas with or without parallel tempering. This would imply that there are no major isolated dominant configurations which can only be reached by employing parallel tempering. All in all, we have not observed any essential improvement by employing parallel tempering in addition to the Hamiltonian Monte Carlo method.
Another interesting aspect of our actual Hamiltonian Monte Carlo simulation is that a relatively smaller choice of the step size seems to give better sampling, and we even took such small values that acceptance rates were nearly 1. This seems to be in contradiction with the more common situation that larger with a reasonable acceptance rate like several 10% would give better sampling. However, this apparent contradictory aspect could be explained in the following manner in our case. The positive semi-definite (2) of the first term in the exponent of the matrix model (1) implies that the dominant configurations for small k are around R i=1 φ i a φ i b φ i c ∼ 0, and this condition becomes tighter as φ i a can take larger values when k is taken smaller. Therefore, the space of dominant configurations can be illustrated as in Figure 18: the dominant configuration space is broad in the small φ i a region, but it becomes narrower as φ i a becomes larger. Here, we also assume that dominant configurations are connected, as suggested in the previous paragraph. Assuming the dominant configuration space as shown in the figure, the updates with relatively smaller will smoothly visit the narrow region as well as the broad region. On the other hand, sampling with relatively larger mainly moves within the broad region, occasionally jumps to the narrow region, and is trapped for a while to compensate the low possibility to visit the narrow region. We have actually observed such trapping to occur more frequently for relatively larger values of . This occasional trapping damages quality of sampling and it generally takes longer time to obtain a dataset with lower margins of error.
Let us summarize this subsection. We encountered the slowdown of the Monte Carlo updates in the region R R c with small k. The speed of updates becomes slower as k becomes smaller, but the dependence is continuous and is subject to the explanation with obvious rescaling. Therefore we have not observed any qualitative changes of the system under the change of the value of k, and it is unlikely that there is a phase transition to a new phase with characteristics of slow dynamics. Rather it seems that the system continues to behave like a fluid with a viscosity which continuously grows for smaller k in the region R R c .

Implications to the tensor model
In this section, we discuss the implications of the results of the simulations to the tensor model in the canonical formalism, the canonical tensor model [23,24].

Phase transition point and the consistency of the tensor model
The wave function of the canonical tensor model, that is obtained by solving a number of first-class constraints to the wave function, 8 has the following form [22], where P abc , a real symmetric tensor, is the configuration variable of the tensor model, λ H = (N + 2)(N + 3)/2, I is the imaginary unit (so I 2 = −1), Ai(·) designates the Airy Ai function, and κ is a real constant in the tensor model. It is particularly important that λ H is determined by the hermiticity condition for the Hamiltonian constraint of the tensor model, and therefore must have this particular form depending on N . Physically, the sign of the parameter κ is supposed to be opposite to that of the cosmological constant, based on the argument relating the mini-superspace approximation of GR and the tensor model with N = 1 [34] (See an appendix of [16] for more details).
The simplest observable for the physical state represented by the wave function (47) would be given by where R = λ H , we have introduced replicas φ i a (i = 1, 2, . . . , λ H ) to replace the power coming from that of (47) in the first line, and have performed the Gaussian integration over P abc . Here α is an arbitrary positive number, N is an unimportant factor independent from α, #P = N (N + 1)(N + 2)/6, i.e. the number of independent components of P abc .
The system (48) is complicated due to the presence of the Airy functions. However, when κ is taken to be positive, which physically corresponds to a negative cosmological constant, the Airy function Ai (κ φ i a φ i a ) is a function that rapidly decays with the increase of φ i a φ i a . Therefore, as an interesting simplification, we could replace the Airy function by a rapidly damping function with a simpler form. In particular, to make the correspondence to the matrix model (1), we consider a simplified wave function, with R = λ H and a positive k by performing the replacement Ai (47). With the observable mentioned above, this leads to where R = λ H .
One important matter in the relation (50) between the tensor and matrix models is that the parameter R of the matrix model (1) is related to N by R = λ H = (N +2)(N +3)/2. What is striking is that this value agrees with the critical value R c ∼ (N + 1)(N + 2)/2 − N + 2 in the leading order of N . Considering the ambiguity of the approximate relation (50), we could say that the tensor model is exactly on or at least in the vicinity of (or a little above of) the continuous phase transition point of the matrix model. This is quite intriguing, because our common knowledge tells that continuum theories can often be obtained by taking continuum limits around continuous phase transition points in discretized theories. We could say that the consistency of the tensor model automatically puts the tensor model at the location where a continuum limit may be feasible, though it is currently difficult to conclude this because of the ambiguity contained in the simplification above.

Dimensions and symmetries of the configurations
In this subsection we will discuss the results of the simulations concerning dimensions and symmetries obtained in Section 6. For this purpose we refer to a property of the wave function (47) that the peaks (ridges) of the wave function are located on the values of P abc which are invariant under Lie-group transformations. This symmetry highlighting phenomenon has been found in [25,26], where the qualitative argument was given as follows. The integration (47) is of an integrand which oscillates rather widely due to the pure imaginary cubic function in the exponent. Therefore, for a "generic" value of P abc , the contributions from different integration spots generally have different phases and mutually cancel among themselves so that the total amount of integration does not take a large value. However, at the location where P abc is invariant under a representation H of a Lie group, P abc = h a a h b b h c c P a b c (∀h ∈ H), the integration along the gauge orbit h a a φ a (∀h ∈ H) contributes coherently in (47), and the wave function has the chance to take a large value compared to that at a "generic" location. This is indeed realized and has concretely been shown for some tractable cases in [25,26].
The above qualitative argument will hold at least partially after the simplification (49), since we can expect a similar coherence phenomenon in this case, too. Then, from the relation (50), the symmetry highlighting phenomenon of the wave function explained above will have a corresponding phenomenon in the matrix model (1) [16]. Note that this will be valid for general values of R, since the constraint R = λ H = (N + 2)(N + 3)/2 coming from the consistency of the tensor model has nothing to do with the equalities in (50). In the relation (50), the contribution of a peak with a Lie group representation H in the second line will correspond on the matrix model side to the contributions of N -dimensional vectors φ i a (i = 1, 2, . . . , R) being distributed along a gauge orbit h a a φ a (∀h ∈ H). In the simulation data, such distributed vectors will appear as a point cloud discussed in Section 6.4. This point cloud will have the dimension of the Lie group representation, and will break part of the SO(N ) symmetry which is not commutative with H. Generally, the wave function contains a number of peaks with various Lie group representations, and therefore the point cloud will be that of a mixture of various gauge orbits. This mixed structure will induce a non-obvious pattern of symmetry breaking, which would be consistent with the hierarchical symmetry breaking in Section 6.5.
To get more information from the behavior obtained in Section 6, let us rewrite the second line in (50) as R #P N a,b,c=1 a≤b≤c dP abc e −αP abc P abc where N , N are some unimportant coefficients. Here, to the second line, we have separated P abc into the radial and angular variables, P abc = PP abc , where P = √ P abc P abc , and have introduced a rescaled variable, ϕ a = P 1/3 φ a . Then, to the last line, we have rescaled P 2 → kP 2 , have divided ϕ a into the radial and angular variables, ϕ a = ϕφ a with ϕ = √ ϕ a ϕ a , and have introducedΨ P ϕ 3 := The functionΨ(P ϕ 3 ) will have a number of peaks at Lie-group symmetricP abc . On such a peak, the value ofΨ(P ϕ 3 ) will generally become smaller as ϕ increases, because the oscillation of the integrand in (52) will become wilder. In the following paragraphs we will further argue that the ϕ-dependence ofΨ(P ϕ 3 ) qualitatively depends on the symmetry ofP abc as in Figure 19 to explain the dimensional behavior in Figure 10: Namely, forP abc symmetric under higher dimensional Lie-groups,Ψ(P ϕ 3 ) takes larger values at small ϕ but quickly decays with ϕ, while it takes smaller values at small ϕ but slowly decays with ϕ forP abc symmetric under lower dimensional Lie-groups.
To see how the dimensional behavior in Figure 10 can be explained by the profile in Figure 19, let us first consider R < (N + 1)(N + 2)/2. In this case, the power of P in the integrand of the last line of (51) is positive, and therefore the integral over P will be over the range 0 ≤ P 1/ √ αk with some preference to larger P . As k is taken smaller, the larger region of ϕ in the integral of (51) becomes more dominant, making the peaks associated with lower dimensional Lie-groups more dominant than higher dimensional ones. Then the increase of the power R in (51) will enhance the peaks of lower dimensional Lie groups. This explains the decrease of the dimensions with the increase of R in the region R < R c in Figure 10.
Let us next consider R > (N + 1)(N + 2)/2. In this case the power of P in the integrand of (51) is negative, and P will have the preference to smaller values as R increases. Then, the last term in (51) will bound the range of ϕ in the integration, as R increases. Because of the profile in Figure 19, increase of R will enhance the peaks with higher dimensional Lie groups, explaining the increase of dimensions in the region R > R c in Figure 10.

Normalizability of the wave function of the tensor model
From the physical point of view, it would be interesting to discuss the norm of the wave function of the tensor model. If the wave function of the tensor model successfully represents a spacetime in some manner, the norm of the wave function will linearly diverge in the time direction, which is supposed to form a trajectory in the space of P abc . Thus, the normalizability of the wave function has a connection to the question concerning time in the tensor model.
As an approximation or as an example case study similar to the actual case, we discuss the norm of the simplified wave function (49). More precisely, we study the α → +0 limit of the relation (50): where N is an unimportant factor independent of α.
When R < R c , by putting (37) with λ = 1/(4α) on the righthand side of (53), the dominant α dependence in the α → +0 limit is obtained as This concludes that the norm diverges for this case by assuming that the critical value satisfies R c < (N + 1)(N + 2)/2.
On the other hand, when R > R c , by putting (40) with λ = 1/(4α) on the righthand side, we obtain under the assumption thatp(+0) is finite, which has been supported from the data in Section 6.3. This is divergent in the limit α → +0 in the range R c < R < (N + 1)(N + 2)/2, since δŨ 0 d > 0, as discussed in Section 6.3. On the other hand, since R = λ H = (N + 2)(N + 3)/2 of the tensor model is in the region R > (N +1)(N +2)/2, we cannot currently determine whether the simplified wave function of the tensor model is normalizable or not or how rapidly this diverges if it does. As explained in Section 6.3, this is beyond the leading order perturbative computation.
The simplification (49) of the real wave function (47) is to approximate the case with a positive κ, which corresponds to a negative cosmological constant. Therefore, it is an interesting future study to determine δŨ 0 d to answer the physical question concerning the emergence of time in the tensor model for a negative cosmological constant. Note also that the above discussion deals with finite N , and therefore taking N → ∞ would also require more study on this matter.

Summary and future prospects
In this paper, we have numerically studied a matrix model with non-pairwise index contractions by Monte Carlo simulations. The matrix model has an intimate connection to the canonical tensor model, a tensor model for quantum gravity in the Hamilton formalism [23,24], and also has a similar structure as a matrix model that appears in the replica trick of the spherical p-spin model (p = 3) for spin glasses [28,29]. The matrix model had previously been analyzed by a few analytic methods and Monte Carlo simulations in [15,16], which had suggested the presence of a continuous phase transition around R ∼ N 2 /2. This relation between N and R is particularly interesting, because this agrees with a consistency condition of the tensor model in the leading order of N , implying that the tensor model is automatically located exactly on or near a continuous phase transition point. However, in the previous works the evidence for the phase transition was not very clear. In this paper we have presented a new set up for Monte Carlo simulations by first integrating out the radial direction, and have studied the model by employing the more efficient Hamiltonian Monte Carlo method. We have obtained considerable improvement of the efficiency of the simulations, and have found a rather sharp continuous phase transition around R = R c ∼ (N + 1)(N + 2)/2 − N + 2. We have also studied various properties of the phase transition and the two phases: the dimensions of the configurations take the smallest values at the transition point; the phase at R > R c is characterized by cascade symmetry breaking; and the k/λ → +0 limit is convergent in one phase and diverges in the other.
We have also discussed some implications to the tensor model. In particular, the most striking is the coincidence above between the location of the continuous phase transition point and the consistency condition of the tensor model in the leading order of N . A well known fact is that continuum theories can often be obtained by taking a continuum limit near a continuous phase transition point. This means that the tensor model seems to automatically put itself at the location where it is possible to find a sensible continuum limit. We have also discussed the wave function of the tensor model by using the connection between the matrix model and an approximation of a known exact wave function of the tensor model. In particular, we have provided a qualitative argument for the dependence on R of both the dimension of the preferred class of configurations and the observed symmetry breaking patterns, using the intimate connection between Lie-group symmetries and peaks of the wave function as has been investigated before [25,26].
While we have numerically obtained a rather clear picture of the phase structure of the matrix model, we are still seriously lacking analytic understanding. As shown in Section 6.2, there seem to exist essential differences between the numerical results and the analytic perturbative results performed in [15,16]. Moreover, other than the qualitative argument given in Section 7.2, a more rigorous understanding of the behavior of the dimensions and the symmetry breaking would be desirable. Analytic understanding is also necessary to discuss the continuum limit discussed in the previous paragraph, since taking a large N limit while simultaneously tuning R and k is difficult to do exclusively through numerical methods. Developing an analytical non-perturbative understanding is an important future direction to understand the dynamics of the matrix model.
Though this paper has given several clear pieces of evidence for the phase transition in the matrix model, explaining its interesting connection to the canonical tensor model, various things still need to be explored in order to understand more about the canonical tensor model through matrix models of the similar sort. Most importantly, the simplification of the wave function discussed in Section 7.1 by approximating the Airy function for κ > 0 by a conveniently chosen damping function does not explain whether the obtained results are universal under a different choice of a damping function. Moreover, from a physical point of view we would like to explore the κ < 0 case corresponding to the positive cosmological constant, rather than the κ > 0 case corresponding to the negative cosmological constant. In the case of κ < 0, the Airy function becomes oscillatory, and the dynamics will most likely be different from the κ > 0 case. Since this case suffers from the notorious sign problem, it is technically very challenging. It would also be an interesting future direction to apply new Monte Carlo methods developed to analyze various other systems suffering from sign problems to our case, as well as to develop analytical treatment.