Numerical and analytical analyses of a matrix model with non-pairwise contracted indices

We study a matrix model that has $\phi_a^i\ (a=1,2,\ldots,N,\ i=1,2,\ldots,R)$ as its dynamical variable, whose lower indices are pairwise contracted, but upper ones are not always done so. This matrix model has a motivation from a tensor model for quantum gravity, and is also related to the physics of glasses, because it has the same form as what appears in the replica trick of the spherical $p$-spin model for spin glasses, though the parameter range of our interest is different. To study the dynamics, which in general depends on $N$ and $R$, we perform Monte Carlo simulations and compare with some analytical computations in the leading and the next-leading orders. A transition region has been found around $R\sim N^2/2$, which matches a relation required by the consistency of the tensor model. The simulation and the analytical computations agree well outside the transition region, but not in this region, implying that some relevant configurations are not properly included by the analytical computations. With a motivation coming from the tensor model, we also study the persistent homology of the configurations generated in the simulations, and have observed its gradual change from $S^1$ to higher dimensional cycles with the increase of $R$ around the transition region.


Introduction
Quantization of gravity is one of the major fundamental problems in theoretical physics. The quantization of general relativity by the standard perturbative methods of quantum field theory fails due to non-renormalizable divergences. Various approaches have been proposed and being studied to solve the fundamental problem, depending on views of authors. In one approach, general relativity (with higher derivative terms) is directly quantized as quantum field theory with the modern technique of the functional renormalization group [1]. In other approaches, fundamental discreteness is introduced to represent spacetimes, which include (causal) dynamical triangulations [2], loop quantum gravity [3], causal sets [4], quantum graphity [5], matrix models [6,7,8,9,10], tensor models [11,12,13,14], and so on. In these discretized approaches, an important criterion for success is whether macroscopic spacetimes are generated, or in other words, whether there exist appropriate continuum limits that recover the usual continuum picture of spacetime with dynamics described by general relativity as low-energy effective theory.
The criterion above can in principle be checked by studying the properties of a wave function of each theory. If the wave function has a peak at a configuration that can well be described by a macroscopic spacetime picture, then the model can be considered to be potentially successful. An indirect motivation for the present paper is to understand the properties of the wave function [15] that is an exact solution to a tensor model in the Hamilton formalism [16,17] (See Appendix A for the tensor model.). It has been argued and shown for some simple cases that the wave function has peaks at the tensor values that are invariant under Lie groups [18]. By using the correspondence developed in [19] between tensor values and spaces with geometries, this would imply that the spacetimes symmetric under Lie groups are favored quantum mechanically. However, the main difficulty in arguing this is that only little part of the peak structure (often called landscape in such contexts) of the wave function is known, not enough to discuss "probabilities of spacetimes." To simplify the problem keeping the main structure from the tensor model as much as possible, one of the authors of the present paper and his collaborators considered the following two simplifications in the former papers. One is that they considered a toy wave function rather than the actual wave function of the tensor model [20]. The actual wave function is expressed by a certain power, say R-th, of a function expressed by an integration over N + 1 variables, but, in the toy wave function, the function is simplified to the one expressed by an integration over N variables by fixing a certain integration variable to a constant. While this substantially simplifies the analysis, the toy wave function keeps the most crucial property mentioned above that there appear peaks at the tensor values that are symmetric under Lie groups as the actual wave function of the tensor model does [20].
Though this toy wave function is simpler than the actual one, it is still difficult to perform thorough analyses, because the dimension of the argument (a symmetric tensor with three indices) of the wave function is very large with the order of ∼ N 3 /6. Therefore, as an additional simplification, the authors in [21] considered a model that can be obtained by integrating over the argument of the toy wave function. This gives a dynamical system of a matrix, say φ i a (a = 1, 2, . . . , N, i = 1, 2, · · · , R), rectangular in general, where the lower indices are 1 pairwise contracted, but the upper ones are not always done so. While this model does not fall into the known solvable models such as the rectangular matrix models [22,23,24] or the vector models [25,26], it has the same form as what appears when the replica trick is applied to the spherical p-spin model for spin glasses [27,28], where R plays the role of the replica number.
Here, though the form is exactly so, the concerned ranges of the dynamical variables and the parameters are different between our model and the spin glass model, and it seems reasonable to reanalyze the matrix model with fresh eyes: (i) While the replica number R is taken to the vanishing limit in the replica trick, it takes a finite value R ∼ N 2 /2 in the correspondence to the tensor model, and should rather be taken to infinity in the thermodynamic limit N → ∞; (ii) A coupling constant 1 takes the opposite sign in our model compared to the spin glass model; (iii) There is a spherical constraint on the dynamical variable in the spherical p-spin model, but there is none in our model.
In this paper, we numerically study the matrix model by the Monte Carlo simulation with the standard Metropolis update method. This contrasts with the perturbative analytical computations performed in the previous paper [21]. We also perform some additional analytical computations to compare with the numerical results. We have obtained the following main results: • The expectation values of some observables are computed by the numerical simulations, and it is observed that there exists a transition region around R ∼ N 2 /2. Intriguingly, the location is in good coincidence with R = (N + 2)(N + 3)/2 that is required by the consistency of the tensor model (Namely, the hermiticity of the Hamiltonian constraints. See Appendix A for more details.) [18,15,29]. Presently, this coincidence is mysterious, since there are no apparent connections between the transition and the consistency. The observables seem to continuously but substantially change their behavior at the transition region, but it has not been determined whether this transition is a phase transition or a crossover in the thermodynamic limit N → ∞. The method for the Monte Carlo simulations performed in this paper is not powerful enough for the determination because of an issue explained below.
• The expectation values of some observables are analytically computed in the leading order, mostly based on the treatment in the previous paper [21], and are compared with the numerical results. Good agreement between them is obtained outside the vicinity of the transition region, while there exist deviations in the transition region. The deviations are such that they soften the transition to make it look more like a crossover. A nextleading order computation has also been performed, but this does not well correct the deviations.
• The tensor model suggests the presence of topological characteristics for the dominant configurations of the matrix φ i a (See Section 7 and Appendix A.). Therefore, we have studied topological characters of the configurations that are generated in the simulations by using the modern technique called persistent homology [30] in the topological data analysis. This technique extracts homology groups possessed by a data, which is a value of the matrix φ i a in our case. The dominant topology gradually changes from S 1 to higher-dimensional cycles with the increase of R in the vicinity of the transition region.
• The Monte Carlo simulation becomes substantially difficult in the region with R N 2 /2 and large λ/k 3 , where λ and k are the parameters of our model (1). In the region, the step sizes of the Metropolis updates chosen for reasonable acceptance rates become too small to reach thermal equilliburiums in ∼ 10 10 Metropolis updates.
• A characterization of the transition can be done by the sizes of the matrix elements, which take relatively large values at the region with R N 2 /2 and large λ/k 3 , but otherwise fluctuate around small values. In the former case, our model may behave like the spherical p-spin model, since the matrix elements are effectively constrained to non-zero sizes, which would approximately realize the spherical constraint in the spherical p-spin model. This may partly explain the bad performance of the Monte Carlo simulation in the region, seemingly reflecting the high viscous nature of glassy fluids.
This paper is organized as follows. In Section 2, we define the model and summarize the previous results [21] that are relevant to the present paper. In Section 3, some observables are introduced and the analytical expressions of their expectation values are obtained in a leading order. In Section 4, the details of the computation in the leading order are given. The result of the next-leading order is also presented, while the details of the computation are given in Appendix D. In Section 5, we perform a saddle point analysis of the expectation values of the observables in the leading order. This describes the transition as a continuous phase transition in the large N limit, where the first derivatives of the expectation values of the observables with respect to R are discontinuous. In Section 6, we compare the Monte Carlo and the analytical results. They agree well outside the transition region. In the transition region, however, there exist deviations, which smoothen the transition to make it look more like a crossover. In Section 7, we analyze the homology structure of the configurations generated by the simulations. The preference changes from S 1 to higher-dimensional cycles with the increase of R in the vicinity of the transition region. The last section is devoted to a summary and future prospects. In Appendix A, the motivation for the matrix model is explained from the viewpoint of the tensor model. In Appendix B, an instructive computation of the partition function for R = 2 is given. In Appendix C, Appendix D, and Appendix E, some equations used in the main text are explicitly derived. In Appendix F, a brief introduction to persistent homology is given.

The model
The partition function of our matrix model is given by where λ and k are the coupling constants assumed to be positive, φ is a (generally rectangular) real matrix, φ i a (a = 1, 2, . . . , N, i = 1, 2, . . . , R), and dφ := N,R a,i=1 dφ i a . The integration is over the N R-dimensional real space denoted by R N R . The coupling terms are defined by where the repeated lower indices are assumed to be summed over. We use this standard convention for the lower indices throughout this paper, unless otherwise stated. On the other hand, we do not use this convention for the upper indices: A sum over them must always be written explicitly 2 . The background motivation for the matrix model is explained in Appendix A from the viewpoint of the tensor model.
In (2), the lower indices are contracted pairwise, while the upper indices are not necessarily so. Therefore the model has the symmetry of the O(N ) transformation on the lower indices, but only the permutation symmetry S R of relabeling {1, 2, . . . , R} on the upper indices. These symmetries are not enough to diagonalize φ i a in general, and therefore this model cannot be solved in a similar way as the usual matrix model.
In the previous paper [21], we considered an expression which can just be obtained by separating the radial and angular part of the integration in (1): By the change of variable, φ i a = rφ i a , with r 2 = R i=1 φ i a φ i a andφ representing the angular coordinates, we obtain where This rather trivial change of expression is actually very useful, because f N,R (t) can be shown to be an entire function of t and therefore has a Taylor expansion in t with the infinite convergence radius around t = 0 (actually around arbitrary t = ∞). Therefore, in principle, the dynamics can be solved by obtaining the entire perturbative series of f N,R (t). Note that the corresponding perturbative expansion of Z N,R (λ, k) in λ around λ = 0, often obtained by perturbative methods, is merely an asymptotic series, because Z N,R (λ, k) has an essential singularity at λ = 0. The f N,R (t) has also the property that it is a decreasing positive function of t with f N,R (0) = 1 for real t. This property provides a good criterion for assessing the 2 One can avoid this unusual convention by introducing some external variables. An example is validity of an approximation of f N,R (t). In the previous paper [21], f N,R (t) in the leading order of 1/R has been determined by a Feynman diagrammatic method with the result, In particular, (5) indeed satisfies the properties above: It is a decreasing function for real t with f 1/R,leading N,R (0) = 1, and is almost an entire function, since the locations of the singularities are far away from the relevant region t ≥ 0 for large N, R.
Since there are two parameters N, R, which can be taken large, the range of validity of (5), which was derived in the leading order of 1/R, is not obvious. However, in later sections, we will find that (5) 3 will give results which agree well 4 with those of the numerical simulations except in the transition region around R ∼ N 2 /2.

Observables
The purpose of this section is to introduce some observables, say O(φ), and discuss the expectation values defined by The observables must respect the symmetry O(N ) × S R mentioned in Section 2. Among various possibilities, we consider Tr(φ t φ) and U (φ) in (2), and also The last one is the diagonal part of the sum in U (φ) in (2). Since these observables are some parts contained in the exponent of (1), they can be implemented in the numerical simulations with little additional computational costs.
To compute the expectation values of these observables, it is convenient to extend (1) by introducing the coupling constant λ d conjugate to U d (φ) as 3 More precisely, because of the difference of our strategy of computations taken in this paper, the expression (17) is slightly different from (5) obtained in [21]. However, the difference is negligible for large N, R, and is not essential. 4 In fact, the expression (5) cannot be correct for small R such as R = 2. This is explicitly shown in Appendix B. However, the difference shown there in the asymptotic region t ∼ ∞ is not relevant for the thermodynamic properties, since for small R, the dominant contributions come from t ∼ 0, as can be explicitly observed in the Monte Carlo simulations in Section 6.
Then the expectation values can respectively be expressed as where , which is the free energy of the model. Here we have put λ d = 0 at last, since our interest is in (1) corresponding to λ d = 0 of (8).
To compute the partition function (8), it is convenient to first separate the radial and the angular part as in (3): where The f N,R,λ,λ d (t) has the similar properties as f N,R (t) explained in Section 2: It is an entire function of t; For λ > 0, λ d ≥ 0, it is positive and decreasing in t for real t; f N,R,λ,λ d (0) = 1. With f N,R,λ,λ d (t), the observables (9) can be expressed by where the normalization is given by From the leading order computation, which is detailed in Section 4, we obtain with h N,R (x) := (1 + 12γ 3 x) − N (N −1)(N +4)

12
(1 + 6(N + 4) where When λ d = 0 is taken, (14) becomes This could be expected to agree with (5), but there is a slight difference coming from (16) with n = 3. This difference originates from the fact that the strategy of the computation we take in Section 4 is different from the one taken previously in [21], and therefore the expressions of the leading orders are slightly different from each other. However, they agree with each other in the leading order of 1/R, since γ n ∼ (N R) −n for N R n, as expected.
By putting these expressions into (12), we obtain where N f is given by (13) with f leading N,R,λ,0 (t). The actual values of these integrations can be obtained numerically.
4 Computations of f N,R,λ,λ d (t) in the leading and the nextleading orders In this section, we will compute f N,R,λ,λ d (t) in the leading order of t, and will also show the result in the next-leading order, whose detailed computations are given in Appendix D. In [21], the computation in the leading order of 1/R has been performed by using the Feynman diagrams for the φ i a variables. In this paper, however, we will take a different strategy. This is because the new strategy makes more transparent the rather complicated counting of combinatorics performed in [21], and make it straightforward to include the extra coupling λ d U d (φ) and also to consider the next order in t. For λ d = 0, the new strategy gives essentially the same result as [21] in the leading order, as commented below (17).
Let us define where with a real symmetric matrix Λ ij . The eigenvalues of the matrix Λ ij are assumed to be nonnegative for the convergence of the corresponding partition function. The f N,R,Λ ij (t) also has the same nice properties as f N,R (t) that it is an entire function, which has a Taylor series expansion in t with the infinite convergence radius around t = 0, and is a decreasing positive function of t for real t and Λ = 0 with f N,R,Λ ij (0) = 1. If we take Λ ij = λ + λ d δ ij , (19) gives f N,R,λ,λ d (t) in (11). By introducing a new variable P i abc (a, b, c = 1, 2, . . . , N, i = 1, 2, . . . , R), which is symmetric for the lower indices, one can rewrite (19) as where I denotes the imaginary unit andΛ ij is a symmetric matrix satisfying 5 , The constant prefactor in (21) can be determined by f N,R,Λ ij (0) = 1.
To compute (21), let us first integrate overφ. This change of the order of the integrations can be done, because the integration over P i abc with the infinite integration region converges uniformly for anyφ ∈ S N R−1 . Then our task is to compute where we have used a short-hand notation, and · φ denotes the expectation value for the uniform probability distribution on the unit sphere S N R−1 .
For further computations, let us introduce the cumulants O n c defined by 6 with arbitrary s. Then (21) can be rewritten as with where S ef f (P ) can be regarded as an effective action of P i abc after integrating outφ, and it is given in terms of the perturbative expansion in t. Due to the form (24), the n-th order cumulant gives the n-th order interaction term of P i abc , and all the terms with odd n vanish because of the obvious invariance of the integration overφ underφ → −φ.
Let us compute the quadratic term with n = 2 in (27). Since Pφ 3 φ = 0, we obtain The integral overφ on S N R−1 , which is a curved compact space, is not easy to handle, so we use a formula which maps this integration to the Gaussian integral: where γ n is defined in (16), and Here β is a sort of dummy variable, which can be chosen freely with β > 0, and does not appear in the final expressions. In fact, as shown below, the factor (2β) m 2 in (29) is exactly canceled by the same factor from the Wick contraction (31). The formula (29) was previously used in [21], and is proven in Appendix C so that the present paper be self-contained. 6 Cumulants are more familiar as connected correlation functions in particle physics, because they can be computed by summing over connected Feynman diagrams. However, we do not use this terminology here, because we will rewrite the cumulants forφ in terms of φ as in (29), and one can explicitly see that these cumulants contain some disconnected diagrams in terms of the Feynman diagrams of φ in general, because of the extra factor γ m/2 in (29).
Through the replacement (29), (28) can be computed by the standard procedure using the Wick theorem and Feynman diagrams. A Wick contraction is performed by which can be derived by an explicit computation of (30). The Feynman diagram for the vertex Figure 1. Each leg is supposed to bring the two indices of φ i a , and a Wick contraction connects two legs with the identification of their indices as in (31). A caution is that each leg on one vertex brings independent lower indices, but a common upper index. (28) with (29). We find the two diagrams shown in the right figure of Figure 1. Their degeneracy factors are 6 and 9, respectively, by counting the numbers of the ways to connect the legs of the two vertices. Since j and j in (28) are identified by the Wick contractions, we also get R j=1Λ ijΛi j = Λ ii as a factor (See (22).). Thus, we obtain

Now let us apply the Wick contractions to what is obtained by replacing
where one notices that the factor (2β) 3 coming from the replacement (29) is exactly canceled by the factors of the three Wick contractions (31) performed for the evaluation. Putting the result (32) into (27), one obtains the effective action in the second order of P i abc as The computation of (26) has now been reduced to diagonalizing the quadratic expression (33). The upper and lower indices can independently be diagonalized, because (33) has the form of the direct product with respect to these indices. More explicitly, since Λ ij is real and symmetric, we can consider the following decomposition into the eigenspaces: where v λev are the orthonormal eigenvectors, and the sum is over all the eigenvalues (with their degeneracies). By putting this and λev v λev i v λev j = δ ij into (33), we obtain a decomposition, where P λev abc := R i=1 v λev i P i abc . Next let us diagonalize the lower index part in (35), where for brevity we omit λ ev from P λev abc . Let us separate P abc into the tensor part P T abc and the vector part P V abc , which are defined by 7 It is easy to check that P T abc P V abc = 0 and P V abc P V abc = 3 N +2 P abb P acc . In particular, the former identity implies that P T abc and P V abc are independent degrees of freedom. Then, by using (37) and the identities above, (36) can be expressed as The number of independent components contained in P T abc and P V abc are #P T = N (N + 1)(N + 2)/6 − N = N (N + 4)(N − 1)/6 and #P V = N , respectively. Therefore, by putting this diagonal form into (26) and integrating over P V and P T , we finally obtain the expression of f N,R,Λ ij (t) from the quadratic order S where we have determined the overall factor by requiring f (2) N,R,Λ ij (0) = 1, the product is over all the eigenvalues of the matrix Λ ij , and h N,R (x) is defined in (15).
For the computation of the observables discussed in Section 3, we consider Λ ij = λ + λ d δ ij . In this case, the matrix Λ ij has one eigenvalue λR + λ d with the eigenvector (1, 1, . . . , 1), and the eigenvalue λ d with degeneracy R − 1 with any of the vectors transverse to (1, 1, . . . , 1) as the eigenvectors. Therefore, from (39), we obtain This is the leading-order result shown in (14).
As we will see later in Section 6, there are some deviations between the leading-order result above and the numerical simulations. To see how the situation changes by adding some corrections, we have also computed the next-leading order. The details of the computation are given in Appendix D. The final result is where f next−leading is the sum of the leading and the next-leading orders, and with the definitions of G i , x i , y i given by from (99) to (109).

Saddle point analysis in the leading order
The integral expressions of the observables (18) in the leading order do not seem to have explicit expressions with known functions. Therefore, a way to obtain their explicit values is to numerically perform the integrations. This will be performed in Section 6 to compare with the Monte Carlo results. In this section, on the other hand, we will apply the saddle point approximation to the integrals to obtain a qualitative global picture of the phase structure of the model.
To discuss the saddle point approximation of the partition function (3), let us consider the minus of the logarithm of the integrand, which is given by 8 with the second derivative being positive at the point. As f N,R , we take (5), which is the expression in the leading order of 1/R: In the saddle point approximation, the free energy of the model is approximated by where the unimportant terms contain log V N R−1 , which does not depend on λ or k, and also some lower order terms in N , which will be discussed in the last paragraph of this section.
Let us first show that there exists a unique solution to the saddle point equation (44), for the leading order expression (45), in the integration region r ≥ 0. To see this, it is convenient to use a new parametrization of R in terms of α as R = R c (1 + α) with R c = (N + 1)(N + 2)/2 and −1 < α. Then, by noting N R c = 6(A 0 + B 0 ), the saddle point equation (44) with (45) can be written as The lefthand side is obviously a decreasing function of r * with a maximum at r * = 0, while the righthand side is an increasing function from zero to the infinity. Since the maximum on the lefthand side is N R c α − 1 + 6A 0 + 6B 0 = N R − 1, there always exists a unique solution of r * for N, R ≥ 1. Moreover, the solution is smooth: r * does not jump in a discrete manner, when the parameters are continuously changed, because the r * -dependence of each side continuously changes. This turns down the possibility that the model has a discontinuous phase transition in this treatment.
To discuss the solution more quantitatively with approximations, let us restrict ourselves to the parameter range of our interest: λ ∼ O(1), N O(10), and k O(1). In addition, for the simplicity of the following discussions, let us avoid the region around α ∼ 0. By noting that N R c , and A 0 are of order O(10 3 ) or larger, one can find that, for each case of α < 0 and α > 0, there are only two relevant terms among all in (48). For α < 0, the first and third terms on the lefthand side are relevant, and for α > 0, the first term on the lefthand side and the one on the righthand side are relevant. By solving the equations under taking these relevant terms only, the solutions are respectively given by The first case shows divergent behavior for α → −0. However, this should not be taken as it is, since the transition should not have such an abrupt behavior as discussed above. In fact, the simplification taken above brakes down in the vicinity of α ∼ 0, and the real behavior is such that r 2 * smoothly interpolates between the two parameter regions in the vicinity of α ∼ 0. The r 2 * in (49) has different large-N behavior in the two regions: N 7 3 for α < 0 and N 3 for α > 0. By normalizing r 2 * with the common factor (N R c ) −1 for both the regions, we obtaiñ in the N → ∞ limit, wherer 2 * := (N R c ) −1 r 2 * . See the leftmost figure in Figure 2. This characterizes the transition at α = 0 as a continuous phase transition with Tr(φ t φ) /(N R c ) = 0 for α ≤ 0 and Tr(φ t φ) /(N R c ) = α 2k for α > 0 in the thermodynamic limit. The other observables can be treated in similar manners. By taking (18), putting r = r * , and taking the leading order in large N , one obtains The divergence of U d (φ) leading in α → −0 should not be taken as it is, because of the same reason mentioned above for r 2 * . By normalizingŨ (φ) : See the middle and rightmost figures in Figure 2. The results support the same conclusion that there is a continuous phase transition at α = 0.
Let us finally comment on the consistency of the above saddle point approximation in the large N case. From (45) and (49), it is straightforward by some explicit computations to find the expansion of F where a n := 1 , and Therefore the integral is dominated by the saddle point value a 0 , and the Gaussian integration around the saddle point 9 and the higher orders in r − r * are subdominant. In addition, the insertion of the observables O as in (18)

Comparison with Monte Carlo simulations
In Section 4, we have computed the leading order approximation of f N,R,λ,λ d (t) defined in (11) in a perturbative method, and have obtained the result (14) along with (15) and so on. We have also derived the next-leading order correction (41) with (42), the details of which are given in Appendix D. With those results, we can numerically calculate the expectation values (e.v.) of the observables given in (12) by the expressions on the righthand sides. However, note that the above approximations of f N,R,λ,λ d (t) are based on taking the perturbative expansion of S ef f given in (27) up to the second order in t. Therefore they seem to require the implicit assumption of small values of t, and may not generally be trusted for the computation of (12), because t is finally assigned with r 6 , and r is integrated over zero to infinity.
In view of the question above, it would be interesting to compute our model without any adoption of approximation methods. More specifically, in this section we compute the e.v. of the observables, (6) , by the Monte Carlo simulations, and compare them with the analytical results obtained by numerically integrating the righthand sides of (12), where we put our perturbative results for f N,R,λ,λ d (t). Note that, in our strategy of the approximations, R is not an expansion parameter, only t is, and therefore it is meaningful to compare the results in the full range of R.      An important thing that can be observed in the figures is that, for each N , there exists a region of R that separates the smaller and larger regions of R with different qualitative behaviors of the observables. This is more clearly seen for larger N and smaller k. The transition region we observe indeed exists around the value R c = (N + 1)(N + 2)/2, which was obtained as the critical point from the saddle point analysis in Section 5 (See Figure 2.).
It is an important physical question whether this transition of behavior is a phase-transition or just a crossover in the thermodynamic limit N → ∞. However, we cannot currently answer this question for certain with the Monte Carlo results presently available, and this would require larger scale Monte Carlo simulations. It seems also difficult to answer this question by our perturbative analytical methods because of the following reason. In the figures, we can find good agreement between the perturbative computations and the Monte Carlo results in the regions away from the transition region. This would support the validity of our perturbative calculations in those outside regions. On the other hand, we can observe that there exist some deviations between the perturbative computations and the Monte Carlo results in the transition region. The deviations appear in such a way that the Monte Carlo results smoothen the transition to make it more like a crossover. Therefore, it seems that the analytical expressions we have obtained as approximations do not seem to be reliable in the transition region.
We can further discuss this complication from another view point as follows. Let us look at the figures more closely. Then we can find that, as to the numerical relations among Trφ t φ , U , and U d in the leading and the next-leading orders, the following hold: Tr φ t φ : including next-leading > leading, U : including next-leading < leading, U d : including next-leading < leading, for all R. Therefore, while the next-leading order corrections indeed improve the approximations so that they approach the Monte Carlo results in the outside regions and this is also so for Tr φ t φ and U in the transition region, the last inequality about U d is in the opposite direction. This suggests that our perturbative treatment seems to have some difficulties in correctly taking into account some configurations that mainly contribute to U d in the transition region. It would be an interesting future problem to identify these configurations.
Let us briefly explain our actual Monte Carlo simulations. We have performed Monte Carlo simulations with the standard Metropolis update method for the model (1) by using KEKCC, the cluster system of KEK. For each calculation shown in the figures, we performed 2 billion sweeps, where the time taken for this was generally about 7 and 23 hours with R = 10 and 130, respectively. We stored the data of the observables once per 400 sweeps, and computed their mean values and the 1σ-errors by the Jackknife resampling method. In each calculation we always set the acceptance rate to be around 60%. However, to realize this 60%, we had to tune the step sizes in our Metropolis method to quite small values, especially in the region R N 2 /2.
Let us further comment on the last peculiar nature we encountered in the simulations. We have performed the simulations for k = 1, 0.1 and 0.05 with N = 10 and k = 1.0 and 0.1 with N = 5, respectively, as shown in the figures. As suggested by the results in the figures, the transition could be sharpened, if we performed simulations with smaller values of k than those in the figures. However, when we tried to do so, we encountered a serious difficulty in particular in the region R N 2 /2. It was that the Metropolis step sizes must be tuned to very small values to keep the reasonable acceptance rate like 60%. Then, the performance of the simulations became so slow that we could not find the timing when the system had reached thermodynamic equilibriums: The system always looked like being in the middle stage of moving very slowly toward thermodynamic equilibriums, at least during one week of continuous running or so. Therefore we took relatively large values of k as those in the figures to avoid the serious difficulty that makes the simulations unreliable.
Finally, let us qualitatively explain why the analytical results computed by the perturbative method and the Monte Carlo results agree with each other outside the transition region of R. Let us start with a qualitative estimation of the effective action (27) with respect to the orders in R and t. One can obtain where we have ignored all the index structures of P i abc for notational simplicity, and the dominant R-dependencies of the coefficients are given by Here the dominant R-dependence of b 0 is obtained by putting the eigenvalue λ ev ∼ R into (36) and recalling γ 3 ∼ 1/R 3 , as defined in (16); The estimation of b 1 is also straightforward and is given in Appendix E. Then, normalizing the quadratic part by rescaling P → P/ √ 1 + b 0 t in (56), one finds that the actual quartic coupling of S ef f is estimated as which is different from the naive value b 1 t 2 . When R is small, r 2 = Tr(φ t φ) is dominated by small values as shown in the Monte Carlo computations above. This means that, since t ∼ r 6 in the usage (3), (58) is also dominated by small values, and therefor the quadratic term S (2) ef f will give good analytic estimations. On the other hand, when R is large, t ∼ r 6 is dominated by large values as shown in the Monte Carlo results. Nonetheless what is remarkable is that (58) is always suppressed by 1/R irrespective of the values of t. Therefore the leading order term S (2) ef f will again give good estimations. In the middle region, however, (58) could generally become large, and higher order and/or non-perturbative corrections may substantially contribute, as suggested by the deviations between the Monte Carlo and analytical results.

Topological structure of configurations
In this section, we will explain our observation on the topological structure of the configurations generated by the Monte Carlo simulations. Topology of a value of the matrix φ i a can be analyzed by persistent homology [30], which is a modern technique of the topological data analysis (See Appendix F for a brief introduction of persistent homology.). More specifically, we performed the Monte Carlo simulations for N = 4 and R = 10, 15, 20, 25 with λ = 1, k = 0.01, and, for each case, uniformly took 100 samples of the values of φ i a during a large number of sequential updates of order 10 8 after thermodynamic equilibriums were seen to be reached. Then, the samples were analyzed in terms of persistent homology. The analysis shows that the favored topology of the configurations is S 1 for R = 10, 15, but gradually changes to higher dimensional cycles, when R is increased. We will first explain the background motivation for this analysis, and will then show the results.
One of the present authors and his collaborators have been studying a tensor model in the canonical formalism, which we call canonical tensor model [16,17], as a model of quantum gravity. In [18], it has been shown that the exact wave function of the tensor model has peaks at the configurations that are invariant under Lie groups. This phenomenon, which we call symmetry highlighting phenomenon, potentially has an important physical significance, since this phenomenon would imply the dominance of spacetimes symmetric under Lie groups through the correspondence between tensors and spaces developed in [19]. This symmetry highlighting phenomenon has first been shown for a toy wave function [20], which slightly simplifies the wave function of the canonical tensor model. The toy wave function is given by 11φ where I denotes the imaginary unit, and is a small positive regularization parameter to assure the convergence of the integral. The symmetry highlighting phenomenon is that the wave function has large peaks at P abc that are invariant under Lie groups: where H is a representation of a Lie group. The phenomenon can qualitatively be understood by the following rough argument: If P abc is invariant under a Lie group, the integration over φ a in (59) will contribute coherently along the gauge orbit h a a φ a ( ∀ h ∈ H), but, if this is not so, the contributions tend to cancel among themselves due to the phase oscillations of the integrand and the wave function takes relatively small values. In [20], some tractable simple cases have explicitly been studied, and the presence of the phenomenon has indeed be shown.
Other than the simple case studies, the peak structure of the toy wave function and that of the tensor model are largely unknown. One reason is that the number of independent components of P abc , which is about ∼ N 3 /6, is so large that it is practically not possible to go over the whole configuration space of P abc . Rather, we will be able to obtain rough knowledge by integrating over P abc : where we have denoted k = −Iκ + , and we have considered an arbitrary power R of the wave function, because the actual wave function of the tensor model has the corresponding power 12 which specifically takes R = (N + 2)(N + 3)/2 [18,15,29] (See Appendix A.). In (60), we see that the integration of a power of the wave function over P abc with a Gaussian weight produces the matrix model (1). If the integration is dominated by such peaks associated to Lie groups, we can expect that the N -dimensional vectors, φ i a (i = 1, 2, . . . , R), in the matrix model tend to exist along some gauge orbits in the vector space.
In general, there exist a number of peaks (or rather ridges) associated with various Lie group symmetries and gauge orbits for the wave function (59). As argued and shown explicitly in [20,18], peaks associated with lower dimensional Lie group symmetries generally exist more abundantly than those with higher dimensional symmetries, because the number of symmetry conditions that must be satisfied by P abc is smaller for the former than for the latter. On the contrary, the peaks of the latter are generally higher than those of the former, because the gauge orbits of the latter have larger dimensions and provide more coherent contributions than the former. Therefore, there are competitions between height and abundancy, and it is generally a subtle question which of lower or higher dimensional Lie group symmetries is probabilistically favored in a given case.
Let us discuss this question in view of the correspondence to the matrix model (60), especially by considering changing the value of R. As R is the power of the wave function as in (60), larger R will enhance higher peaks compared to the lower ones. Therefore we expect that, for larger R, the contributions from peaks with higher dimensional symmetries will be more dominant. Otherwise, peaks with lower dimensional Lie group symmetries will dominate because of their abundant existence. In addition, when R is very small, non-symmetric configurations will dominate, since they exist most abundantly.
In the Monte Carlo simulation of our model, if a symmetric peak dominates as explained above, the N -dimensional vectors φ i a (i = 1, 2, . . . , R) will be randomly distributed along the associated gauge orbit. A caution here is that the gauge orbit can take any O(N )-transformed location in the N -dimensional vector space due to the O(N ) symmetry of the model. Therefore we should not simply plot all the samples of φ i a generated by the Monte Carlo simulations in the N -dimensional vector space, since this will merely provide a trivial O(N )-invariant spherical distribution of points without any characteristics of the Lie-group symmetry associated to the dominant peak. Rather, we have to analyze each sample of φ i a generated by the Monte Carlo simulations to extract its characteristics, and then pile up all the extractions over all the samples to find allover characteristic properties.
Topological structure of each sample of φ i a can be analyzed by using the technique of persistent homology [30]. This is a modern applied mathematical technique of the topological data analysis, and can extract homology groups of a data (See Appendix F.). Here an input data should be a set of points with relative distances. We used an open-source c++ program that is called Ripser 13 for the analysis and plotted the output with Mathematica. For a configuration φ i a , we consider the replica number i = 1, 2, . . . , R to represent the label of "points" of a data set, and define the distances between two points i and j as To avoid the dependence of the initial values, 10 independent Monte Carlo sequences were run, and the sampling were performed uniformly from the sequence of updates of ∼ 10 9 after thermodynamic equilibriums were seen. 100 configurations of φ i a were uniformly sampled and the persistent homologies were analyzed (one-to three-dimensional homologies from the left to the right). The results of 100 samples are plotted on the same persistent diagrams. Blue dots represent the longest-life elements in each dimensional persistent homology group of each data, the yellow ones the second, and the green ones all after the second. The dots away from the diagonal line represent long-life persistent homology group elements, which are considered to be characteristics of a data, while those near the diagonal line are regarded as "noises". The highest blue dots, namely those with the largest u end that represent the largest structure, move from H 1 to H 2 and then H 3 with the increase of R. 23 The gist of this definition is that the N -dimensional vectors φ i a (i = 1, 2, . . . , R) are projected onto the unit sphere S N −1 , and the geodesic distances along the sphere are regarded as the distances. In particular, this definition is suited for detecting a gauge orbit h b a φ b (h ∈ H), since it is projected on the sphere irrespective of the size of the vector φ a .
We want to see the phenomenon explained above in the actual Monte Carlo simulations. For this initial study, choosing small N would be preferred for simper analysis, because then there exist a small number of possibilities of gauge orbits with small dimensions, and also thermodynamic equilibriums can easily be reached due to the small number of degrees of freedom. Note however that this trades off the clarity of the homology structure being detected by persistent homology. This is because, for small N , the values of R at which the phenomenon appears are also rather small in the order of R ∼ N 2 /2, as we will see in the analysis below. The homology cycles formed by small numbers (namely R) of points necessarily become obscure, especially higher dimensional cycles are difficult to be clearly detected.
For the actual simulation, we considered N = 4. In N = 4, as explicitly solved in [20], there exist only two possibilities of Lie group symmetries, SO(2) and SO (3), and the gauge orbits are S 1 and S 2 , respectively. In fact, the ridges of the P abc with these symmetries reach the origin P abc = 0, and therefore we should also add the trivial possibility of S 3 with the SO(4) symmetry, which is the symmetry of P abc = 0 and is maximal. Figure 7 shows the persistent diagrams obtained from the Monte Carlo simulations with R = 10, 15, 20, 25. Statistically speaking, one can observe that, starting from S 1 at R = 10, 15, higher dimensional cycles gradually appear and become the largest structure when R is increased, while lower dimensional cycles gradually take smaller values of u.

Summary and future prospects
In this paper, we studied a matrix model containing non-pairwise index contractions [21], which has a motivation from a tensor model of quantum gravity [16,17]. More specifically, it has φ i a (a = 1, 2, . . . , N, i = 1, 2, . . . , R) as its degrees of freedom, where the lower indices are pairwise contracted, but the latter are not always done so. This matrix model has the same form as what appears in the replica trick of the spherical p-spin model for spin glasses [27,28], though the variable and parameter ranges of our interest are different. We performed Monte Carlo simulations with the Metropolis update method, and compared the results with some analytical computations in the leading order, mostly based on the previous treatment in [21]. They are in good agreement outside the transition region, that is located around R ∼ N 2 /2. In the transition region, however, there exist deviations between the simulations and the analytical results, and the deviations cannot be corrected well, even if the next-leading order contributions are included. It has not been determined whether the transition is a phase transition or a crossover in the thermodynamic limit N → ∞, because of the limited range of the parameters like N 10 available in our Monte Carlo simulation. Our Monte Carlo simulation tended to slow down especially at R N 2 /2 with large λ/k 3 , suspecting that the system gets glassy nature in the region, but no conclusive argument has been made for this aspect. We also studied the topological characteristics of the configurations φ i a generated in the Monte Carlo simulations by using the modern technique called persistent homology [30] in topological data analysis. This technique extracts the homology structure of a data, which is a configuration of φ i a in our case. We observed that, in the vicinity of the transition region, the homology structure of the configurations gradually changes from S 1 to higher-dimensional cycles with the increase of R.
A particularly interesting result of this paper is that there seems to exist a transition region around R ∼ N 2 /2. Intriguingly, this value of R coincides with what is required by the consistency of the tensor model (namely, the hermiticity of the hamiltonian constraint. See Appendix A.) [18,15,29]. Moreover, our model seems to have the most interesting properties in this region, but they are not well understood: There are some deviations between the simulation and the analytical results in this region, but the reason is not clear; The transition of the homological structure of the dominant configurations in this vicinity is peculiar but not well understood; Whether the transition is a phase transition or a crossover in the thermodynamic limit N → ∞ is not determined. For the better understanding in the future, it seems necessary to treat larger N cases with large λ/k 3 by employing more efficient methods of Monte Carlo simulations and finding more powerful analytical methods. these tensor models under the interpretation employed in these original papers [11,12,13,14], the numbers of the indices of the tensorial dynamical variables have direct connections to the dimensions of the spaces supposed to be emergent. This is a drawback from the viewpoint of quantum gravity, since the number of dimensions should also be dynamically determined rather than predetermined as an input.
This drawback may be overcome by introducing a new interpretation of the tensor models. In fact, tensors themselves have rich structures enough to describe spaces. This can be seen by supposing that a three-index tensor, say C ab c , define the structure constants of an algebra of functions, say f a , over a space by f a f b = C ab c f c . It is known that spaces can be defined by the algebras of functions on them satisfying certain properties (see for example [32].). Though such algebras for usual spaces are commutative and associative, we would be able to suppose that more general C ab c , which defines noncommutative or/and nonassociative algebras, would define "fuzzy spaces". This idea is explored to propose a new interpretation of three-index tensor models as models of dynamical fuzzy spaces [33]. Under this new interpretation, a threeindex tensor model can in principle generate various dimensional spaces, not being restricted to certain dimensions [33]. One can also see the emergence of Euclidean general relativity on them [34,35,36,37].
However, these favorable results were obtained only by doing some fine-tunings to the models [34,35,36,37]. Moreover, time is missing in these tensor models. In fact, time might play essentially important roles in the mechanism of emergence of spacetime, as the success of the causal dynamical triangulation [2] shows over the dynamical triangulation. Generally speaking, introducing time requires delicate treatment in quantum gravity, since such models must respect the spacetime diffeomorphism invariance on emergent spacetimes. To introduce time with this delicate requirement, a tensor model, which we call the canonical tensor model (CTM), has been formulated as a constrained system in the Hamilton formalism [16,17]. The model is composed of a number of first-class constraints analogous to those in the ADM formalism of general relativity. Therefore, CTM embodies an analogue of the spacetime diffeomorphism as its gauge symmetry. CTM is supposed to describe time-evolutions of fuzzy spaces.
The classical aspects of CTM have been compared with GR in a formal continuum limit with N → ∞: The constraint algebra of CTM has been shown to agree with that of the ADM formalism of GR [38]; the classical equation of motion of CTM has been shown to agree with that of GR in the Hamilton-Jacobi formalism with a certain Hamilton's principal function [39]. Simple cases with large but finite N have also been studied, and mutual agreement has been obtained [19].
The canonical quantization of CTM is straightforward [29]. The classical constraints are transformed to the quantized ones with one additional term coming from operator ordering. The algebra essentially remains the same form, especially implying the closure of the commutation algebra among the quantized constraints. This assures the consistency of the physical state conditions given byĤ a |Ψ phys =Ĵ ab |Ψ phys = 0, whereĤ a andĴ ab denote the quantized constraints of CTM, respectively corresponding to the Hamiltonian and momentum constraints of ADM. By taking a certain representation, these equations become a system of partial differential equations for physical wave functions. At first sight they looked too complicated to solve, but it has turned out that there exist various explicit exact solutions to these equations [15]. In particular, a systematic solution is given in the P -representation by [15,18] where λ H = (N + 2)(N + 3)/2, and with P φ 3 ≡ P abc φ a φ b φ c , φ 2 ≡ φ a φ a , and λ is an arbitrary real parameter 14 , which can be normalized as 0, ±1. The above value of λ H is uniquely determined by the hermiticity of the "Hamilotonian constraints"Ĥ a of CTM.
In principle, emergence of spacetime can be analyzed by investigating the properties of the wave function Ψ(P ) phys . The ideal situation would be that, in the configuration space, the quantum probability density |Ψ(P ) phys | 2 form some ridges that can be considered to be trajectories of classical time-evolutions of spaces. In fact, it has been qualitatively argued and actually been found in some tractable simple cases that such ridges exist at P abc 's that are invariant under Lie-groups 15 . Since symmetries are ubiquitous in the actual spacetimes 16 , this result is encouraging in the pursuit for spacetime emergence. However, further analysis of Ψ(P ) phys is not easy, and most of the properties of the wave function are still unknown. Therefore, it is presently not possible to discuss emergent spacetimes in CTM.
Yet, aiming for better understanding of Ψ(P ) phys , we consider in this paper a model which can be obtained after two simplifications. One is to integrate the quantum probability density over the configuration space: where dP ≡ N a≤b≤c dP abc , P 2 ≡ P abc P abc and α is an arbitrary positive number. The other is to fix (or ignore) the integration variableφ of ϕ(P ) in (64) so thatφ(P ) in (59) is used as ϕ(P ) in (63). The wave function after this replacement is still interesting enough, becauseφ(P ) has the similar connection between Lie-group symmetries and peaks [20]. Then, by puttingφ(P ) into (63) and performing the Gaussian integration over P in (65), one obtains where we have used the reality of the wave functionφ(P ) for the real value of k = −Iκ + .
Here we have actually obtained the partition function Z in (1) of the matrix model with the particular choice R = λ H (= (N + 2)(N + 3)/2). Since λ H ∼ R c (= (N + 1)(N + 2)/2), CTM seems to be located in the transition region of the matrix model.

Appendix B R = 2 case
In this appendix, we consider the partition function for R = 2, and see that the partition function is finite even for k = 0. This is not trivial, because, for general R > 1, the solution to U (φ) = 0 is non-empty (see below), and therefore there is a potential risk of runaway behavior, Let us first see that U (φ) = 0 is non-empty for general R > 1. Since An obvious set of solutions are given by φ 2i−1 a = −φ 2i a , (i = 1, 2, R/2 ) with φ R a = 0 if R=odd. For R = 2, one can see that this is the only solution as follows. By contracting the indices b and c in (68), we obtain This implies that φ 1 and φ 2 are linearly dependent, and putting this back to (68), we obtain φ 1 a = −φ 2 a . In the R = 2 case, since U (φ) depends only on the relative directions and the sizes of φ 1 a and φ 2 a , the partition function for k = 0 can be written as where U (r 1 , r 2 , θ) = r 6 1 + r 6 2 + 2r 3 1 r 3 2 cos 3 (θ) with r 1 and r 2 being the sizes, and θ being the relative angle. As shown above, the only case with U = 0 is given by r 1 = r 2 , θ = π. Therefore let us perform the reparameterization, Then the integral of (70) at large r can be approximated by expanding the integrand for small x, y, and we obtain ∼ drdxdy r 2N −1 y N −2 e −3λr 6 (3x 2 +y 2 ) ∼ dr r −N −1 .
This shows that the partition function is convergent for R = 2 and k = 0.
There are two things which can be learned from this simplest case. One is that, if R is small enough, the partition function is convergent even for k = 0. Another thing is that the large-r asymptotic behavior of the integrand is much slower than the leading order result, Therefore the asymptotic behavior derived from the leading order result cannot be correct down to R = 2.
Appendix C Derivation of (29) In this section, we derive (29). This was previously derived in [21], but this is repeated here to make the present paper self-contained.
For m =odd, the equation holds, because the both sides vanish.
Let us assume m = 2p with a positive integer p. Let us start with the following equation: where β is an arbitrary positive constant. This can easily be proven by reparameterizing φ i a with the radial and the angular variables as φ i a = rφ i a on the righthand side, and observing that the integrations over r decouples from the angular part and cancels between the numerator and the denominator.
Let us consider the numerator on the righthand side of (75), Taking the p-th derivative of A(β) with respect to β cancels the (Trφ t φ) p in the denominator of the inetegrand. On the other hand, by performing the rescaling φ i a → β −1/2 φ i a , it is obvious that A(β) has the dependence β −N R/2 on β. Therefore, by performing the p-th derivative of the both sides of (76), we obtain the relation, By solving for A(β) and putting it into (75), we obtain (29). Appendix D Computation of f N,R,Λ ij (t) in the next-leading order In this appendix, we will compute the fourth-order term (Pφ 3 ) 4 c φ in (27). From the definition (25) of cumulants and the formula (29), we obtain The (P φ 3 ) 2 c φ in the last term has already been computed in Section 4. As for (P φ 3 ) 4 c φ , Wick contractions (31) give the five connected diagrams in Figure 8. By counting the number of ways to connect the legs, one can find that the degeneracies are given by 2 3 · 3 5 , 2 4 · 3 5 , 2 3 · 3 4 , 2 4 · 3 4 , and 2 3 · 3 5 , respectively, from the left to the right diagrams 17 in Figure 8. These Feynman diagrams represent the ways of the index contractions of P i abc 's in the fourth order interaction term. For example, the leftmost diagram gives where the numerator of the numerical factor is the degeneracy, and the denominator comes from the factor of the Wick contraction (31). It is also straightforward to write down the explicit expressions for all the other diagrams in Figure 8. Now let us suppose we have obtained the explicit expressions of S (4) ef f (P ) by the above procedure. An immediate difficulty of this forth order term is that S (4) ef f (P ) has the negative all over sign due to I 4 as in (27), and therefore the system with S ef f (P ) = S (2) ef f (P )+S (4) ef f (P ) is where the summation over σ denotes the sum over all the permutations of d, e, f , and the product of two matrices, say X and Y , is defined by Note that A acts as the identity on a symmetric tensor, namely, (AP ) abc = P abc . One can easily check the following properties: Furthermore, one can check that A − B and B give the projectors to the tensor and vector parts of P abc , respectively, as Therefore, by using these matrices, (38) can be rewritten in the form, where c T,λev and c V,λev are the coefficients associated to the tensor and vector parts of P abc , namely, c T,λev = 1 + 12γ 3 λ ev t, c V,λev = 1 + 6(N + 4)γ 3 λ ev t.
The inverse of the matrix in (88) is given by For later convenience, let us define Let us next take into account the upper indices. To derive the above result we started with the expression (36) with an eigenvalue λ ev of the matrix Λ ij . Therefore it is the result for the corresponding eigenspace. Considering the projection to each eigenspace, the final form of the Wick contraction for P i abc is obtained as where M ij λev denotes the projector to the eigenspace of the eigenvalue λ ev of Λ ij . Let us next start computing S (4) ef f (P ) P by using the Wick contraction (92). Let us first compute the factors coming from the projectors. Let us restrict ourselves to the case Λ ij = λ + λ d δ ij , which is of our interest as explained in Section 4. As given there, the eigenvector for λ ev = λR + λ d is (1, 1, · · · , 1), and those for λ ev = λ d are the vectors transverse to that. Therefore the projectors are M ij λR+λ d = 1/R for λ ev = λR + λ d , and M ij λ d = δ ij − 1/R for λ ev = λ d , respectively. For later usage, let us compute the projectors sandwiched betweenΛ ij : where we have used the following explicit solution to (22) for Λ ij = λ + λ d δ ij : As for the factor coming from theΛ ij in (79), the Wick contractions (92) of P i abc will generate the factor j (ΛM λevΛ ) jj (ΛM λ evΛ ) jj . Thus, using (93), we obtain the following results for each case: for λ ev = λ ev = λR + λ d , It is easy to see that all the other diagrams in Figure 8 have the same factor.
As for the Wick contractions (92) of (P φ 3 ) 2 c φ 2 in (78), there exist two cases. One is to take the contractions within each (P φ 3 ) 2 c φ , or otherwise. The former case can be called the disconnected case, and the latter the connected case, based on their diagrammatic characters. From (32) and (92), the factors coming from the Λ ij 's are obtained as for the disconnected and the connected cases, respectively.
The last ingredient for the computation of S ef f (P ) is to take into account the lower index part of the Wick contraction (92). As can be seen in (92), in general, this is to insert xA + yABA with some x, y at the location of the Wick contraction (See Figure 10 for an example). Diagrammatically, this is to insert the diagrams in Figure 9 at the location of the Wick contraction. This insertion generates many diagrams with a number of loops. The number of loops gives the degeneracy of each diagram in powers of N . The summation over all the diagrams is too many to do so by hand, so we performed this task by using Mathematica. We have obtained for (P φ 3 ) 4 c φ . As for (P φ 3 ) 2 c φ 2 , we have obtained G 3 (x, y) =2N (2 + N )(4 + N )(15x 2 + 24xy + 6N xy + 8y 2 + 6N y 2 + N 2 y 2 ), respectively, for the disconnected and connected cases.
Appendix E Estimation of the R-dependence of b 1 As for b 1 , there are two kinds of contributions as given in (78). The dominant R-dependence of the former term of (78) can be estimated by considering the expression (79), because the contractions of the upper indices among P i abc are the same for all the diagrams in Figure 8. It is sufficient to consider λ d = 0, since this is of our interest, and the mode of P i abc contributing in this case with the appropriate normalization has the form 1 √ R (1, 1, . . . , 1)P abc , which has the form of the eigenvector of the matrix Λ ij = λ for the eigenvalue λR. SinceΛ ij ∼ 1/ √ R (see (94)) and there are five free summations over the upper indices in (79), the former term of (78) is estimated as As for the latter term of (78), (P φ 3 ) 2 c φ ∼ RP 2 by a similar argument as above. Therefore, b latter where it is important to note that the contents of the parentheses give R −7 due to the cancellation of the leading R −6 terms of both.
The above two estimations conclude b 1 ∼ R −5 .

Appendix F Brief introduction of persistent homology
In this appendix, we give a brief introduction of persistent homology for this paper to be self-contained. More details can be found for instance in [30].
Persistent homology is a notion that characterizes the topological aspects of a data in terms of homology. A data to be analyzed should be a set of points with relative distances. From a data, a stream (or a filtration) of simplicial complexes parameterized by a scale parameter, say u, is constructed. The details of the construction of a stream will be given at the end of this appendix. Once a stream is constructed, the homology groups of the simplicial complexes at each value of u are computed. By increasing the value of u from zero, a homology group element will appear at a certain value of u, say u start , and will continue to exist until it disappears at another value, say u end . If the life period u end − u start is large, one may regard the element as a persistent homology group element, which has a long life. The collection of persistent homology group elements characterizes the topological property of a data set. There will also be a number of short-life elements, but they are often regarded as "noises", which are not robust against small perturbations of the data. Roughly speaking, the scale u parameterizes the sizes of topological structure of interest, and persistent homology characterizes a data with multi-scale homology groups.
There are two kinds of diagrams that are convenient for visualizing the persistent homology of a data. One is called barcode diagram, where each horizontal line segment [u start , u end ] represents a homology group element that exists during the period. An example of barcode diagrams constructed from a sample of φ i a generated in the Monte Carlo simulation is shown in Figure 11. How to construct a set of points with relative distances from φ i a is given in Section 7. The left figure shows the barcode diagram for the 0-dimensional homology, and the right that for the 1-dimensional homology. The left diagram indicates that the initially separated points form one connected component over the scale u ∼ 1.4. The right figure shows that there exists a one-dimensional cycle which has the size of u ∼ 1.8, while there is a small "noise" around u ∼ 0.9. In Figure 12, we provide the graph constructed by connecting the points with relative distances u < 1.5, using the same φ i a for Figure 11. One can actually see the presence of a one-dimensional cycle consistent with the barcode diagram.
The other kind of diagram is called persistent diagram. An element that is represented by a line segment [u start , u end ] in a barcode diagram is represented by a dot located at the twodimensional coordinate (u start , u end ) in a persistent diagram. Since there are multiple elements in general, and u start < u end , a persistent diagram consists of a number of dots in the region over the diagonal line. The long-life elements are represented by the dots that exist away from the diagonal line, and those in the vicinity of the diagonal line are regarded as "noises". An example of a persistent diagram is given in Figure 13, which corresponds to the right barcode diagram in Figure 11. What is convenient in a persistent diagram is that one can easily superimpose persistent diagrams from multiple data. If there is a common characteristics through multiple data, one can find it as a characteristic pattern in a superimposed persistent diagram. Therefore, we use persistent diagrams to find statistically favored structure common in the configurations generated by the Monte Carlo simulation.  The data is the same one used for Figure 11. This choice of the distance cut-off is made because a one-dimensional cycle is expected to exist at this scale, seeing the right of Figure 11. Figure 13: The persistent diagram corresponding to the right barcode diagram in Figure 11. A line segment in a barcode diagram is represented by a dot in a persistent diagram. The longest one is represented by a blue dot, and the second by a yellow one. The colors are also used in the persistent diagrams in Section 7.
Let us finally explain the actual construction of a stream of simplicial complexes parameterized by u. In fact, there exist various streams depending on purposes in the literature, but let us restrict ourselves to the Vietoris-Rips stream, which is used in the open-source c++ code called Ripser. For a given data that stores distances between points, the Vietoris-Rips stream VR(V, u) is defined as follows: • The vertex set is given by the point set V of a data.
• For vertices i and j with distance d(i, j), the edge [i, j] are included in VR(V, u), if and only if d(i, j) ≤ u.
• A higher-dimensional simplex is included in VR(V, u), if and only if all of its edges are.
From the above definition, there is an obvious property, VR(V, u) ⊂ VR(V, u ) for u < u . An important fact is that this induces a map: H k (VR(V, u)) → H k (VR(V, u )) for u < u . Therefore, the development of each homology group element under the change of the value of u can be followed, and its life is characterized by the two endpoint values of u. This makes barcode and persistent diagrams convenient ways for the description.