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

We study a matrix model that has ϕai(a=1,2,…,N,i=1,2,…,R)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _a^i\ (a=1,2,\ldots ,N,\ i=1,2,\ldots ,R)$$\end{document} 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∼N2/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\sim N^2/2$$\end{document}, 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 S1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S^1$$\end{document} 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 thea e-mail: sasakura@yukawa.kyoto-u.ac.jp b e-mail: shingo.takeuchi@phenikaa-uni.edu.vn ory 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 Ref. [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 Ref. [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 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: (1) 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 → ∞; (2) a coupling constant 1 takes the opposite sign in our model compared to the spin glass model; (3) 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 1 More specifically, λ in (1). This sign difference originates from the pure imaginary value of the tensor coupling in (59), which is real in the spherical p-spin model. 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.) [15,18,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 next-leading 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 Sect. 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 equillibriums 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 Sect. 2, we define the model and summarize the previous results [21] that are relevant to the present paper. In Sect. 3, some observables are introduced and the analytical expressions of their expectation values are obtained in a leading order. In Sect. 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 Sect. 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 Sect. 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 Sect. 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 (1) 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φ : 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, representing the angular coordinates, we obtain 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 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, 2 One can avoid this unusual convention by introducing some external variables. An example is 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 Sect. 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 Then the expectation values can respectively be expressed 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 Ref. [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 Sect. 6. 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 Sect. 2: It is an entire function of t; For λ > 0, λ d ≥ 0, it is positive and decreasing in t for real , the observables (9) can be expressed by where the normalization is given by From the leading order computation, which is detailed in Sect. 4, we obtain (14) with 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 Sect. 4 is different from the one taken previously in Ref. [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.

Computations of f N,R,λ,λ d (t) in the leading and the next-leading 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 Ref. [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 Ref. [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 Ref. [21] in the leading order, as commented below (17).
Let us define where with a real symmetric matrix i j . The eigenvalues of the matrix i j are assumed to be non-negative for the convergence of the corresponding partition function. The f N ,R, i j (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, i j (0) = 1. If we take i j = λ + λ d δ i j , (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˜ i j is a symmetric matrix satisfying, 5 The constant prefactor in (21) can be determined by 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 .
with arbitrary s. Then (21) can be rewritten as with where S e f 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: 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).
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 selfcontained.
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 Fig. 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. Now let us apply the Wick contractions to what is obtained by replacing (28) with (29). We find the two diagrams shown in the right figure of Fig. 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˜ i j˜ i j = ii as a factor (see 22). Thus, we obtain 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

S
(2) 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 i j 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 (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 P abc = P T abc + P V abc , 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 7 This decomposition can be understood as follows. First P abb represents the O(N )-vector part of P abc . Its embedding to P abc is given by an expression proportional to (37). Then the coefficient can be determined by the condition that P T abc does not contain the vector part: P T abb = 0.
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, i j (t) from the quadratic order S (2) where we have determined the overall factor by requiring f (2) N ,R, i j (0) = 1, the product is over all the eigenvalues of the matrix i j , and h N ,R (x) is defined in (15).
For the computation of the observables discussed in Sect. 3, we consider i j = λ + λ d δ i j . In this case, the matrix i j 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 Sect. 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 is the sum of the leading and the nextleading orders, and with the definitions of G i , x i , y i given by from (D.21) to (D.31).

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 Sect. 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 (43) 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 Eq. (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 Fig. 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 Fig. 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 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 9 dr e −a2(r −r * ) 2 generates a contribution proportional to log a 2 O(log N ) to the free energy, and this is subdominant.
insertion of the observables O as in (18)

Comparison with Monte Carlo simulations
In Sect. 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 e f 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.
The results of the comparison are summarized in Figs   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 Sect. 5 (see Fig. 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 e f 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) e f 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) e f 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 Ref. [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 Ref. [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 Ref. [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 [15,18,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 Refs. [18,20], 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 with-out 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 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 Ref. [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 simula- 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 threedimensional homologies from the left to the right). The results of 100 samples are plotted on the same persistent diagrams. Blue dots repre-sent 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 longlife 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 tions 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 nonpairwise 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 Ref. [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") [15,18,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.
Acknowledgements The Monte Carlo simulations in this study were performed using KEKCC, the cluster system of KEK. The work of N.S. is supported in part by JSPS KAKENHI Grant no. 19K03825.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: We can provide the data and the programs on the basis of each request from interested researchers.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendices Appendix A: The motivation for the model (1) from the viewpoint of the canonical tensor model
In this appendix, we explain the motivation for the model (1) by summarizing the developments of the canonical tensor model so far.
An interesting direction in the attempt to formulate quantum gravity is to pursue the existence of the spacetime as an emergent phenomenon. In this context, the matrix models have achieved great success in the description of the twodimensional quantum gravity [31]. To extend the success to higher dimensions, the tensor models [11][12][13][14] were proposed as generalization of the matrix models. Though various interesting results have been obtained, emergence of spacetimes has not been achieved so far by these tensor models, suffering from the dominance of singular configurations instead of macroscopic space-like ones. Moreover, in 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 three-index 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, (A.1) 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 Refs. [15,18] (P) phys = ϕ(P) λ H /2 , where λ H = (N + 2)(N + 3)/2, and 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 d P ≡ N a≤b≤c d P 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 (A.3) so thatφ(P) in (59) is used as ϕ(P) in (A.2). 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 (A.2) and performing the Gaussian integration over P in (A.4), 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, φ 2 → ∞ with U (φ) = 0. 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 For R = 2, one can see that this is the only solution as follows. By contracting the indices b and c in (B.2), we obtain This implies that φ 1 and φ 2 are linearly dependent, and putting this back to (B.2), 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 (θ ) (B.5) 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 (B.4) 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 .

(B.7)
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 Ref. [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 = 2 p 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 (C.1), 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 (C.2), we obtain the relation, By solving for A(β) and putting it into (C.1), we obtain (29).
where f (2) N ,R, i j (t) is given in (39), the allover constant has been determined by requiring f (4) N ,R, i j (0) = 1, and · P is defined by with the quadratic action S (2) e f f (P) given in (33). The computation of f (4) N ,R, i j (t) in (D. 3) has now been reduced to that of S (4) e f f (P) P . This can be computed by the Wick theorem, using the Wick contraction for P i abc determined from S (2) e f f (P). The Wick contraction of P i abc can be obtained by taking the inverse of the coefficient matrix in the quadratic action S (2) e f f (P) of P i abc . Since S (2) e f f (P) has the form of the direct product with respect to the upper and lower indices, they can be treated separately.
Let us first treat the lower indices by starting with (36). Let us introduce the following matrices (see Fig. 9): 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 V,λ ev = 1 + 6(N + 4)γ 3 λ ev t. (D.12) The inverse of the matrix in (D.11) is given by (D.13) For later convenience, let us define (D.14) 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 i j . 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 i j λ ev denotes the projector to the eigenspace of the eigenvalue λ ev of i j .
Let us next start computing S (4) e f f (P) P by using the Wick contraction (D.15). Let us first compute the factors coming from the projectors. Let us restrict ourselves to the case i j = λ + λ d δ i j , which is of our interest as explained in Sect. 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 i j λR+λ d = 1/R for λ ev = λR + λ d , and M i j λ d = δ i j − 1/R for λ ev = λ d , respectively. For later usage, let us compute the projectors sandwiched between˜ i j : where we have used the following explicit solution to (22) for i j = λ + λ d δ i j : As for the factor coming from the˜ i j in (D.2), the Wick contractions (D.15) of P i abc will generate the factor j (˜ M λ ev˜ ) j j (˜ M λ ev˜ ) j j . Thus, using (D.16), we obtain the following results for each case: for λ ev = λ ev = λR + λ d , , for λ ev = λ ev = λ d , , otherwise.

(D.18)
It is easy to see that all the other diagrams in Fig. 8 have the same factor. As for the Wick contractions (D.15) of (Pφ 3 ) 2 c φ 2 in (D.1), 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 (D.15), the factors coming from the i j 's are obtained as  for the disconnected and the connected cases, respectively. The last ingredient for the computation of S (4) e f f (P) is to take into account the lower index part of the Wick contraction (D.15). As can be seen in (D.15), in general, this is to insert x A + y AB A with some x, y at the location of the Wick contraction (see Fig. 10 for an example). Diagrammatically, this is to insert the diagrams in Fig. 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 where   A graph made by connecting the points with relative distances smaller than 1.5. The data is the same one used for Fig. 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 Fig. 11 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 Fig. 12, we provide the graph constructed by connecting the points with relative distances u < 1.5, using the same φ i a for Fig. 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 two-dimensional 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 Fig. 13, which corresponds to the right barcode diagram in Fig. 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 Fig. 13 The persistent diagram corresponding to the right barcode diagram in Fig. 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 Sect. 7 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.
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.